#!/usr/bin/env python3
"""RevenueMix assessment v1. Standard library only; CLI never uses the network.

python3 revenue.py data.csv --date-column date --revenue-column revenue --drivers ads sales
python3 revenue.py --request request.json --output assessment.json
"""
import argparse
import csv
import datetime as dt
import io
import json
import math
from pathlib import Path
import sys

VERSION = '1.0.0'
MAX_BYTES = 200_000
MAX_ROWS = 2000
MAX_DRIVERS = 8


def _dot(a, b):
    return math.fsum(x*y for x, y in zip(a, b))


def _solve(r, b):
    out = [0.0]*len(b)
    for i in range(len(b)-1, -1, -1):
        out[i] = (b[i] - _dot(r[i][i+1:], out[i+1:])) / r[i][i]
    return out


def _fit(x, y):
    """Centered/scaled OLS using twice-reorthogonalized modified Gram-Schmidt."""
    n, k = len(x), len(x[0])
    means = [math.fsum(row[j] for row in x)/n for j in range(k)]
    scales = [math.sqrt(math.fsum((row[j]-means[j])**2 for row in x)/n) for j in range(k)]
    if any(s <= 1e-12*max(1, abs(m)) for s, m in zip(scales, means)):
        raise ValueError('A driver is constant or has too little variation. Remove it or supply more varied observations.')
    z = [[1.0]+[(v-m)/s for v, m, s in zip(row, means, scales)] for row in x]
    p = k+1
    q, r = [], [[0.0]*p for _ in range(p)]
    for j in range(p):
        col = [row[j] for row in z]
        for _ in range(2):
            for i in range(j):
                projection = _dot(q[i], col)
                r[i][j] += projection
                col = [v-projection*u for v, u in zip(col, q[i])]
        r[j][j] = math.sqrt(_dot(col, col))
        if r[j][j] < 1e-8*math.sqrt(n):
            raise ValueError('Drivers are linearly dependent or nearly indistinguishable. Remove redundant drivers.')
        q.append([v/r[j][j] for v in col])
    scaled_beta = _solve(r, [_dot(col, y) for col in q])
    beta = [scaled_beta[0]-math.fsum(b*m/s for b, m, s in zip(scaled_beta[1:], means, scales))]
    beta += [b/s for b, s in zip(scaled_beta[1:], scales)]
    fitted = [_dot(row, scaled_beta) for row in z]
    residual = [v-f for v, f in zip(y, fitted)]
    inv_columns = [_solve(r, [float(i==j) for i in range(p)]) for j in range(p)]
    inv_r = list(map(list, zip(*inv_columns)))
    gram_inv = [[_dot(a,b) for b in inv_r] for a in inv_r]
    leverage = [math.fsum(col[i]**2 for col in q) for i in range(n)]
    # HC3 sandwich covariance: sum (G z_i)(G z_i)' e_i²/(1-h_i)².
    cov = [[0.0]*p for _ in range(p)]
    for row, err, h in zip(z, residual, leverage):
        influence = [_dot(g, row)*err/max(1e-12,1-h) for g in gram_inv]
        for a in range(p):
            for b in range(p):
                cov[a][b] += influence[a]*influence[b]
    transform = [[1.0]+[-m/s for m,s in zip(means,scales)]]
    transform += [[0.0]+[1/scales[j] if i==j else 0.0 for j in range(k)] for i in range(k)]
    se = [math.sqrt(max(0.0, math.fsum(t[a]*cov[a][b]*t[b] for a in range(p) for b in range(p)))) for t in transform]
    return beta, residual, se, [n*gram_inv[j][j] for j in range(1,p)], leverage


def _beta_fraction(a, b, x):
    # Continued fraction for the regularized incomplete beta function.
    tiny=1e-300
    c=1.0; d=1-(a+b)*x/(a+1); d=1/(d if abs(d)>tiny else tiny); h=d
    for m in range(1,201):
        for aa in [m*(b-m)*x/((a+2*m-1)*(a+2*m)), -(a+m)*(a+b+m)*x/((a+2*m)*(a+2*m+1))]:
            d=1+aa*d; d=d if abs(d)>tiny else tiny
            c=1+aa/c; c=c if abs(c)>tiny else tiny
            d=1/d; delta=d*c; h*=delta
        if abs(delta-1)<3e-14:
            return h
    raise ValueError('Uncertainty calculation did not converge.')


def _ibeta(a,b,x):
    if x<=0: return 0.0
    if x>=1: return 1.0
    front=math.exp(math.lgamma(a+b)-math.lgamma(a)-math.lgamma(b)+a*math.log(x)+b*math.log1p(-x))
    if x<(a+1)/(a+b+2): return front*_beta_fraction(a,b,x)/a
    return 1-front*_beta_fraction(b,a,1-x)/b


def _critical_t(df):
    lo,hi=0.0,100.0
    for _ in range(70):
        mid=(lo+hi)/2
        tail=_ibeta(df/2,0.5,df/(df+mid*mid))
        if tail>0.05: lo=mid
        else: hi=mid
    return (lo+hi)/2


def _parse(payload):
    if not isinstance(payload,dict): raise ValueError('Expected a JSON object.')
    unknown=set(payload)-{'csv','date_column','revenue_column','driver_columns'}
    if unknown: raise ValueError('Unknown fields: '+', '.join(sorted(unknown)))
    text=payload.get('csv'); date_col=payload.get('date_column'); target=payload.get('revenue_column'); drivers=payload.get('driver_columns')
    if not isinstance(text,str) or not text.strip(): raise ValueError('Provide csv as a nonempty string.')
    if len(text.encode('utf-8'))>MAX_BYTES: raise ValueError(f'CSV must be at most {MAX_BYTES} UTF-8 bytes.')
    if not isinstance(drivers,list) or not 1<=len(drivers)<=MAX_DRIVERS or not all(isinstance(s,str) and s for s in drivers):
        raise ValueError('driver_columns must contain 1–8 column names.')
    if not isinstance(date_col,str) or not isinstance(target,str): raise ValueError('Specify date_column and revenue_column.')
    if 'intercept' in drivers: raise ValueError('Rename the reserved driver column name intercept.')
    selected=[date_col,target]+drivers
    if len(set(selected))!=len(selected): raise ValueError('Date, revenue, and each driver must use different columns.')
    try:
        reader=csv.reader(io.StringIO(text.lstrip('\ufeff')),strict=True)
        header=next(reader)
        if len(header)>40 or len(set(header))!=len(header) or any(not h.strip() or len(h)>100 for h in header):
            raise ValueError('Use up to 40 unique, nonempty column headers of at most 100 characters.')
        missing=[c for c in selected if c not in header]
        if missing: raise ValueError('Missing columns: '+', '.join(missing))
        indices=[header.index(c) for c in selected]
        records=[];dates=set()
        for line,row in enumerate(reader,2):
            if not row: continue
            if len(records)>=MAX_ROWS: raise ValueError(f'Use at most {MAX_ROWS} observations.')
            if len(row)!=len(header): raise ValueError(f'CSV row {line}: expected {len(header)} cells, received {len(row)}.')
            try:
                date_text=row[indices[0]].strip()
                date=dt.date.fromisoformat(date_text)
                if date.isoformat()!=date_text: raise ValueError()
            except ValueError: raise ValueError(f'CSV row {line}: use a valid YYYY-MM-DD date.') from None
            if date in dates: raise ValueError(f'CSV row {line}: duplicate date {date}; aggregate to one observation per period.')
            dates.add(date)
            values=[]
            for column,index in zip(selected[1:],indices[1:]):
                try: value=float(row[index])
                except ValueError: raise ValueError(f'CSV row {line}, {column}: supply a number; missing values are not imputed.') from None
                if not math.isfinite(value) or abs(value)>1e12: raise ValueError(f'CSV row {line}, {column}: use a finite number with magnitude at most 1e12.')
                values.append(value)
            records.append((date,values[0],values[1:]))
    except (csv.Error, StopIteration) as exc:
        raise ValueError('Malformed or empty CSV.') from exc
    return records,target,drivers


def assess(payload):
    records,target,drivers=_parse(payload)
    n=len(records);p=len(drivers)+1;minimum=max(12,5*p)
    result={'schema_version':VERSION,'tool':'revenue_assessment','status':'completed','summary':'',
            'data':{'observations':n,'drivers':drivers,'revenue_column':target,'minimum_observations':minimum},
            'diagnostics':{},'coefficients':[],'warnings':[],
            'limitations':['Associations are not causal effects or incremental ROI. Confounding, reverse causality, trends, and omitted variables can explain estimates.',
                          'The model is additive and contemporaneous; it does not estimate adstock, saturation, seasonality, or delayed effects.',
                          'HC3 intervals address heteroskedasticity, not serial correlation. They use an approximate Student-t reference and may be unreliable in small samples.'],
            'next_steps':[],
            'review':{'url':'https://statwonk.com/contact/?service=revenue-assessment','label':'Request an expert model review','requires_user_decision':True,'pricing':'Consulting is quoted separately. This assessment is free.'}}
    if n<minimum:
        result.update(status='insufficient_data',summary=f'Provide at least {minimum} observations for {len(drivers)} drivers. Received {n}.')
        result['next_steps']=[f'Add at least {minimum-n} complete observations at a consistent frequency, or reduce the number of drivers.','Include independently varying drivers and relevant controls. These minimums are screening rules, not a guarantee of adequate statistical power.']
        return _report(result)
    if records!=sorted(records,key=lambda row:row[0]): result['warnings'].append('Rows were sorted by date before assessment and holdout evaluation.')
    records.sort(key=lambda row:row[0]);dates=[r[0] for r in records];y=[r[1] for r in records];x=[r[2] for r in records]
    result['data'].update(first_date=dates[0].isoformat(),last_date=dates[-1].isoformat())
    gaps={(b-a).days for a,b in zip(dates,dates[1:])}
    monthly=all((b.year-a.year)*12+b.month-a.month==1 and (a.day==b.day or ((a+dt.timedelta(days=1)).month!=a.month and (b+dt.timedelta(days=1)).month!=b.month)) for a,b in zip(dates,dates[1:]))
    if len(gaps)>1 and not monthly: result['warnings'].append('Dates are irregularly spaced or contain missing periods; align the time series before interpreting results.')
    avg=math.fsum(y)/n;tss=math.fsum((v-avg)**2 for v in y)
    if tss<=1e-20*max(1,math.fsum(v*v for v in y)):
        result.update(status='needs_revision',summary='Revenue is constant or nearly constant; there is no useful variation to model.')
        result['next_steps']=['Check the revenue column and provide observations with meaningful variation.']
        return _report(result)
    try: beta,residual,se,vif,leverage=_fit(x,y)
    except ValueError as exc:
        result.update(status='needs_revision',summary=str(exc));result['next_steps']=['Remove constant or redundant drivers and rerun.']
        return _report(result)
    sse=_dot(residual,residual);df=n-p;critical=_critical_t(df)
    dw=math.fsum((b-a)**2 for a,b in zip(residual,residual[1:]))/sse if sse>1e-20 else None
    result['diagnostics']={'r_squared':1-sse/tss,'adjusted_r_squared':1-(sse/df)/(tss/(n-1)),
        'residual_rmse':math.sqrt(sse/df),'residual_degrees_of_freedom':df,'durbin_watson':dw,
        'max_leverage':max(leverage),'variance_inflation_factors':dict(zip(drivers,vif)),
        'interval_method':'HC3 covariance with approximate Student-t 95% intervals', 'holdout':None}
    for name,b,s in zip(['intercept']+drivers,beta,se):
        result['coefficients'].append({'name':name,'estimate':b,'standard_error':s,'interval_95':[b-critical*s,b+critical*s]})
    if n<50: result['warnings'].append('Fewer than 50 observations: treat estimates and uncertainty intervals as preliminary.')
    if max(vif)>5: result['warnings'].append('Some drivers have VIF above 5; correlated drivers make individual estimates unstable.')
    if max(leverage)>min(0.99,3*p/n): result['warnings'].append('High-leverage observations may dominate the fitted model; inspect unusual periods.')
    if dw is not None and not 1.5<=dw<=2.5: result['warnings'].append('Residuals show possible serial dependence. HC3 intervals do not correct this; consider a time-series model and serial-correlation-robust uncertainty.')
    if sse/tss<1e-10: result['warnings'].append('Near-perfect fit: check for target leakage or synthetic/deterministic data.')
    train_n=int(n*0.8)
    if n>=max(30,10*p) and train_n>=minimum:
        try:
            train_beta,*_=_fit(x[:train_n],y[:train_n]);test_y=y[train_n:]
            predicted=[_dot([1.0]+row,train_beta) for row in x[train_n:]]
            baseline=math.fsum(y[:train_n])/train_n
            model_rmse=math.sqrt(math.fsum((a-b)**2 for a,b in zip(test_y,predicted))/len(test_y))
            baseline_rmse=math.sqrt(math.fsum((v-baseline)**2 for v in test_y)/len(test_y))
            result['diagnostics']['holdout']={'training_observations':train_n,'test_observations':len(test_y),'rmse':model_rmse,'training_mean_baseline_rmse':baseline_rmse,'beats_baseline':model_rmse<baseline_rmse,'method':'Final 20% of dates, one fixed-origin split. Uses observed holdout drivers; not an ex-ante revenue forecast.'}
            if model_rmse>=baseline_rmse: result['warnings'].append('On the held-out dates, the model did not beat predicting the training-period mean revenue.')
        except ValueError:
            result['warnings'].append('Holdout fit was not identifiable in the training period; collect more varied history.')
    else: result['warnings'].append('Too few observations for the holdout screen (requires at least max(30, 10 × parameter count)); predictive usefulness is untested.')
    result['summary']='An exploratory revenue model was fitted. Review diagnostics and limitations before making a business decision.'
    result['next_steps']=['Check the source data, unusual periods, and whether each driver was known at the time of the decision.',
                          'Test stability across more time periods and include relevant controls before interpreting driver estimates.',
                          'For spending decisions, establish a causal design or request expert review; do not treat coefficients as guaranteed returns.']
    return _report(result)


def _report(result):
    lines=['# RevenueMix assessment','',result['summary'],'',f"Observations: {result['data']['observations']}"]
    if result['coefficients']:
        lines+=['','| Term | Estimate | Approximate 95% interval |','|---|---:|---|']
        for c in result['coefficients']:
            name=c['name'].replace('|','\\|').replace('\n',' ')
            lines.append(f"| {name} | {c['estimate']:.6g} | [{c['interval_95'][0]:.6g}, {c['interval_95'][1]:.6g}] |")
        lines+=['',f"In-sample R²: {result['diagnostics']['r_squared']:.3f}",f"Residual RMSE: {result['diagnostics']['residual_rmse']:.6g}"]
    for title,key in [('Warnings','warnings'),('Limitations','limitations'),('Next steps','next_steps')]:
        lines+=['',f'## {title}']+['- '+s for s in result[key]]
    lines+=['','Expert review is optional and quoted separately: '+result['review']['url']]
    result['report_markdown']='\n'.join(lines)+'\n'
    return result


def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('csv_path',nargs='?');parser.add_argument('--request');parser.add_argument('--date-column',default='date');parser.add_argument('--revenue-column',default='revenue');parser.add_argument('--drivers',nargs='+');parser.add_argument('--output')
    args=parser.parse_args()
    try:
        if args.request: payload=json.loads(Path(args.request).read_text(encoding='utf-8'))
        elif args.csv_path:
            payload={'csv':Path(args.csv_path).read_text(encoding='utf-8-sig'),'date_column':args.date_column,'revenue_column':args.revenue_column,'driver_columns':args.drivers}
        else: parser.error('Provide a CSV path or --request JSON path.')
        result=assess(payload);text=json.dumps(result,indent=2,allow_nan=False)
        if args.output: Path(args.output).write_text(text+'\n',encoding='utf-8')
        else: print(text)
    except (ValueError,OSError) as exc:
        print(json.dumps({'error':str(exc)}),file=sys.stderr);return 2
    return 0

if __name__=='__main__':
    sys.exit(main())