"""Reproduce selected analyses from the user's course repositories.
Usage: python scripts/analytics/reproduce.py CASE --repo /path/to/original/repo
CASE is s-mobile, fitech, or pentathlon. No customer-level data is exported.
Quarto reads precomputed figures; it does not execute training during a build.
"""
import os
os.environ.setdefault('OPENBLAS_NUM_THREADS', '1')
os.environ.setdefault('OMP_NUM_THREADS', '1')
from pathlib import Path
import argparse, hashlib, json, subprocess, platform
import numpy as np
import pandas as pd
import polars as pl
import scipy
import statsmodels.api as sm
import statsmodels.formula.api as smf
import sklearn
from sklearn.metrics import roc_auc_score, roc_curve
ROOT=Path(__file__).resolve().parents[2]
PROJECTS={'s-mobile':'s-mobile-churn','fitech':'fitech-offer-design','pentathlon':'pentathlon-next-product'}

def save(case,repo,files,result):
    result['provenance']={'commit':subprocess.check_output(['git','-C',str(repo),'rev-parse','HEAD'],text=True).strip(),
      'files':{f:hashlib.sha256((repo/f).read_bytes()).hexdigest() for f in files},
      'python':platform.python_version(),'numpy':np.__version__,'pandas':pd.__version__,
      'polars':pl.__version__,'scipy':scipy.__version__,'statsmodels':sm.__version__,'sklearn':sklearn.__version__}
    out=ROOT/'projects'/PROJECTS[case]/'assets'
    out.mkdir(parents=True,exist_ok=True)
    (out/'results.json').write_text(json.dumps(result,indent=2,allow_nan=False)+'\n')
    print(json.dumps({k:v for k,v in result.items() if k not in ['scores','roc','deciles','importance','provenance']},indent=2),flush=True)


def fitech(repo):
    from itertools import combinations
    hist=pl.read_excel(repo/'data/exhibits.xlsx',sheet_name='exhibit1').to_pandas()
    offers=pl.read_excel(repo/'data/exhibits.xlsx',sheet_name='exhibit2').to_pandas()
    assert (hist.resp+hist.non_resp==hist.nr_emailed).all()
    hist['response_rate']=hist.resp/hist.nr_emailed
    model=smf.glm('response_rate ~ apr + annual_fee + bk_score + C(fixed_var)',data=hist,
                 family=sm.families.Binomial(),freq_weights=hist.nr_emailed).fit()
    scored=[]
    for bk in [150,200,250]:
        d=offers[['offer','apr','fixed_var','annual_fee']].copy();d['bk_score']=bk
        d['p_open']=model.predict(d);d['clv']=offers[f'clv{bk}'].to_numpy()
        d['ev_per_email']=d.p_open*d.clv;d['margin_per_email']=d.ev_per_email-.5
        scored.append(d)
    scores=pd.concat(scored,ignore_index=True)
    margin={(int(r.offer),int(r.bk_score)):r.margin_per_email for r in scores.itertuples()}
    def cost(n):return 800+10000+1000*max(0,n-1)
    solutions=[]
    for k in range(1,13):
        best=None
        for subset in combinations(range(1,13),k):
            pick={bk:max(subset,key=lambda o:margin[o,bk]) for bk in [150,200,250]}
            net=sum(250000*margin[pick[bk],bk] for bk in pick)-cost(k)
            if best is None or net>best['net']:
                best={'offers':list(subset),'pick':pick,'net':float(net),'design_and_round_cost':cost(k)}
        solutions.append(best)
    best=max(solutions,key=lambda r:r['net'])
    grouped=hist.groupby('annual_fee').agg(responses=('resp','sum'),sent=('nr_emailed','sum'))
    grouped['response_rate']=grouped.responses/grouped.sent
    save('fitech',repo,['fitech.ipynb','data/exhibits.xlsx'],{
        'historical_cells':len(hist),'historical_contacts':int(hist.nr_emailed.sum()),'historical_responses':int(hist.resp.sum()),
        'coefficients':model.params.to_dict(),'scores':scores.to_dict('records'),'best':best,
        'best_single':solutions[0],'best_by_offer_count':solutions,
        'observed_fee_groups':grouped.reset_index().to_dict('records'),
        'gross_expected_clv':best['net']+best['design_and_round_cost']+375000,
        'contact_cost':375000,'prospects':750000,
        'notes':['Weighted binomial GLM; grouped count likelihood equivalent to original weighted opened/not-opened rows.',
                 'Enumerates all 4095 nonempty offer subsets. No Round 1 feedback data supplied; static projection only.']})


def smobile(repo):
    from sklearn.ensemble import GradientBoostingClassifier
    from sklearn.preprocessing import StandardScaler,OneHotEncoder,FunctionTransformer
    from sklearn.pipeline import Pipeline,make_pipeline
    from sklearn.compose import ColumnTransformer
    from sklearn.inspection import permutation_importance
    data=pl.read_parquet(repo/'data/s_mobile.parquet').to_pandas()
    numeric=['changer','changem','revenue','mou','overage','roam','conference','months','uniqsubs','custcare','retcalls','dropvce','eqpdays']
    categorical=['refurb','smartphone','highcreditr','mcycle','car','travel','region','occupation']
    logcols=['revenue','mou','roam','conference','uniqsubs','custcare','retcalls','dropvce']
    signed=['changer','changem'];linear=['overage','months','eqpdays']
    for col in numeric:data[col]=data[col].astype(float)
    data['y']=(data.churn=='yes').astype(int)
    train=data[data.training==1].copy();test=data[data.training==0].copy();rep=data[data.representative==1].copy()
    assert set(train.customer).isdisjoint(test.customer) and set(rep.customer).isdisjoint(set(train.customer)|set(test.customer))
    assert not data[numeric+categorical].isna().any().any()
    def correct(p,p_true=.02,p_sample=.5):
        p=np.clip(np.asarray(p),1e-12,1-1e-12)
        odds=p/(1-p)*(p_true/(1-p_true))/(p_sample/(1-p_sample))
        return odds/(1+odds)
    # Preserve the source formula. Drop unsupported Logit freq_weights: correction is post-hoc once.
    formula='y ~ '+' + '.join([f'np.log1p({c})' for c in logcols]+[f'np.sign({c})*np.log1p(np.abs({c}))' for c in signed]+linear+[f'C({c})' for c in categorical])
    lr=smf.logit(formula,train).fit(maxiter=300,disp=False)
    print('S-Mobile logistic fitted',flush=True)
    pre=ColumnTransformer([
       ('log1p',make_pipeline(FunctionTransformer(np.log1p),StandardScaler()),logcols),
       ('signed',make_pipeline(FunctionTransformer(lambda x:np.sign(x)*np.log1p(np.abs(x))),StandardScaler()),signed),
       ('lin',StandardScaler(),linear),('cat',OneHotEncoder(handle_unknown='ignore',sparse_output=False),categorical)])
    gb=Pipeline([('prep',pre),('clf',GradientBoostingClassifier(n_estimators=300,max_depth=4,learning_rate=.05,subsample=.8,random_state=42))])
    gb.fit(train[numeric+categorical],train.y)
    print('S-Mobile gradient boosting fitted',flush=True)
    lr_test=lr.predict(test);gb_test=gb.predict_proba(test[numeric+categorical])[:,1]
    raw=gb.predict_proba(rep[numeric+categorical])[:,1];cal=correct(raw)
    rep['p_raw']=raw;rep['p']=cal
    ordered=rep.sort_values('p',ascending=False).copy()
    ordered['decile']=np.repeat(np.arange(1,11),len(rep)//10)
    deciles=ordered.groupby('decile').agg(n=('y','size'),churns=('y','sum'),observed=('y','mean'),predicted=('p','mean')).reset_index()
    importance=permutation_importance(gb,test[numeric+categorical],test.y,n_repeats=3,random_state=42,scoring='roc_auc',n_jobs=1)
    curves={}
    for name,p in [('Logistic',lr_test),('Gradient boosting',gb_test)]:
        fpr,tpr,_=roc_curve(test.y,p);indices=np.linspace(0,len(fpr)-1,151,dtype=int)
        curves[name]={'fpr':fpr[indices].tolist(),'tpr':tpr[indices].tolist()}
    top=ordered.head(3000)
    save('s-mobile',repo,['s-mobile.ipynb','data/s_mobile.parquet'],{
      'samples':{'training':len(train),'test':len(test),'representative':len(rep)},
      'observed_base_rate':float(rep.y.mean()),'representative_churns':int(rep.y.sum()),
      'auc':{'logistic_test':float(roc_auc_score(test.y,lr_test)),'boosting_test':float(roc_auc_score(test.y,gb_test)),
             'boosting_representative':float(roc_auc_score(rep.y,cal))},
      'mean_probability':{'raw':float(raw.mean()),'corrected':float(cal.mean()),'logistic_corrected':float(correct(lr.predict(rep)).mean())},
      'top_decile':{'n':len(top),'churns':int(top.y.sum()),'rate':float(top.y.mean()),'capture':float(top.y.sum()/rep.y.sum()),'lift':float(top.y.mean()/rep.y.mean())},
      'deciles':deciles.to_dict('records'),'roc':curves,
      'importance':sorted([{'feature':c,'mean':float(m),'std':float(s)} for c,m,s in zip(numeric+categorical,importance.importances_mean,importance.importances_std)],key=lambda x:-x['mean']),
      'notes':['Original 300-tree configuration and train/test split; representative sample never used to fit models.',
               'Single post-hoc prior correction; unsupported statsmodels Logit freq_weights argument removed.',
               'Permutation importance rerun with 3 repeats (original 10); no intervention or CLV causal claims.']})


def pentathlon(repo):
    data=pl.read_parquet(repo/'data/pentathlon_nptb.parquet').to_pandas()
    messages=['control','backcountry','endurance','racquet','strength','team','water']
    features=['age','female','income','education','children','freq_endurance','freq_strength','freq_water','freq_team','freq_backcountry','freq_racquet']
    assert not data[features+['buyer','total_os','message','training']].isna().any().any()
    train=data.training==1;test=data.training==0;y=(data.buyer=='yes').to_numpy().astype(int)
    assert set(data.loc[train,'custid']).isdisjoint(data.loc[test,'custid'])
    # All features interact with message in the source: algebraically equivalent separate message models.
    X=pd.get_dummies(data[features],drop_first=True).astype(float)
    means=X.loc[train].mean();std=X.loc[train].std().replace(0,1)
    X=sm.add_constant((X-means)/std).to_numpy()
    Xt=X[test];yt=y[test];msg=data.message.astype(str).to_numpy()
    probs=[];orders=[];aucs={};best_uniform_training={};metadata=[]
    for m in messages:
        subset=np.asarray(train)&(msg==m);buyers=subset&(y==1)
        model=sm.GLM(y[subset],X[subset],family=sm.families.Binomial()).fit(maxiter=100)
        order_model=sm.OLS(data.total_os.to_numpy()[buyers],X[buyers]).fit()
        p=model.predict(Xt);order=np.maximum(order_model.predict(Xt),0)
        probs.append(p);orders.append(order)
        matched=(msg[test]==m);aucs[m]=float(roc_auc_score(yt[matched],p[matched]))
        best_uniform_training[m]=float(np.mean(model.predict(X[train])*np.maximum(order_model.predict(X[train]),0)*.4))
        metadata.append({'message':m,'training_rows':int(subset.sum()),'training_buyers':int(buyers.sum()),'converged':bool(model.converged)})
        print('Pentathlon fitted',m,flush=True)
    P=np.column_stack(probs);O=np.column_stack(orders);EP=P*O*.4
    best_profit=EP.argmax(axis=1);best_response=P.argmax(axis=1)
    uniform=messages.index(max(best_uniform_training,key=best_uniform_training.get))
    rng=np.random.default_rng(455);random_indices=rng.integers(0,7,len(data))[np.asarray(test)]
    idx=np.arange(len(Xt))
    strategies={'No message':float(EP[:,0].mean()),'Random assignment':float(EP[idx,random_indices].mean()),
                'Best uniform':float(EP[:,uniform].mean()),'Response targeting':float(EP[idx,best_response].mean()),
                'Value targeting':float(EP[idx,best_profit].mean())}
    dist=[{'message':m,'response_share':float(np.mean(best_response==j)),'value_share':float(np.mean(best_profit==j))} for j,m in enumerate(messages)]
    save('pentathlon',repo,['pentathlon_nptb.ipynb','data/pentathlon_nptb.parquet'],{
       'n':len(data),'train':int(train.sum()),'test':int(test.sum()),'buyers':int(y.sum()),'training_buyers':int(y[train].sum()),
       'models':metadata,'holdout_auc':aucs,'uniform_choice':messages[uniform],
       'uniform_expected_contribution':dict(zip(messages,EP.mean(axis=0).tolist())),
       'strategies':strategies,'allocation':dist,'changed_choice_share':float(np.mean(best_response!=best_profit)),
       'margin':.4,'campaign_scale':5000000,
       'notes':['Separate message-specific logistic and conditional OLS fits reproduce the fully interacted source models.',
                'Best uniform message selected using training predictions, not test outcomes.',
                'Conditional order predictions clipped at zero as in source; 40% gross margin; contact costs omitted.',
                'All strategy amounts are model-implied contribution on held-out features, not observed policy effects.']})

if __name__=='__main__':
    parser=argparse.ArgumentParser();parser.add_argument('case',choices=PROJECTS);parser.add_argument('--repo',type=Path,required=True)
    args=parser.parse_args();{'s-mobile':smobile,'fitech':fitech,'pentathlon':pentathlon}[args.case](args.repo.resolve())