#!/usr/bin/env python3
"""Adjacent-decile fusion. Offline, versioned, nested region tuning and full refits.
Run: prepare; fit; bootstrap --workers 4; report. Prior studies remain read-only.
"""
import argparse, csv, gzip, hashlib, importlib.util, json, lzma, math, shutil, time, zipfile
from concurrent.futures import ProcessPoolExecutor
from pathlib import Path
import numpy as np
import scipy
from scipy.optimize import minimize

ROOT=Path(__file__).resolve().parents[1]
OUT=ROOT/'public/research/social-factors/deciles-fusion-v3'
WORK=ROOT/'.fusion-work'
def module(name,path):
    s=importlib.util.spec_from_file_location(name,path);v=importlib.util.module_from_spec(s);s.loader.exec_module(v);return v
d=module('deciles',ROOT/'scripts/social-deciles-v1.py');m=d.m
save=m.save;read=m.read;sha=m.sha
FEATURES=list(d.FEATURES)
MODELS=['B','C-linear','C-smooth','DE','Q','F']
RATIOS=[1.,.3,.1,.03,.01,.003,0.]
SEED=202610103;REPLICATES=199;TOLERANCE=.005;ZERO=1e-7

def compact(path,value):
    path.parent.mkdir(parents=True,exist_ok=True)
    b=json.dumps(value,separators=(',',':'),allow_nan=False).encode()
    if path.suffix=='.xz':b=lzma.compress(b)
    elif path.suffix=='.gz':b=gzip.compress(b,mtime=0)
    tmp=path.with_suffix('.tmp');tmp.write_bytes(b);tmp.replace(path)
def uncompact(path):
    b=path.read_bytes()
    return json.loads(lzma.decompress(b) if path.suffix=='.xz' else gzip.decompress(b) if path.suffix=='.gz' else b)
def base():return read(d.OUT/'frozen-design.json')
def verify():
    p=read(OUT/'frozen-design.json')
    for x in p['inputs']:assert sha(ROOT/x['path'])==x['sha256'],'Frozen input changed: '+x['path']
    return base()

def prepare():
    assert not (OUT/'frozen-design.json').exists(),'Already frozen; amendment required'
    OUT.mkdir(parents=True,exist_ok=True)
    b=d.verify()
    protocol={
      'id':'social-deciles-fusion-v3','registeredAt':d.now(),'status':'Exploratory; 2023–2025 outcomes and earlier model results already inspected',
      'hypothesis':'Training-selected adjacent fusion retains useful adjusted differences with fewer effective social coefficients and improves generalisation.',
      'cohort':'Identical 133 exact-code English authorities, 2018–2025; 399 temporal and 133 geographic assessment rows; original frozen temporal and nine leave-region-out folds.',
      'outcome':'Police-reported personal-injury collision events at road location, including non-motor-only events; log(all-motor miles/100000) offset. Not claims or resident-driver risk.',
      'features':'Income and employment deprivation, and 2011 no-qualifications share. Original historical vintages, baseline controls and training-only history prior. No interactions.',
      'deciles':'Original inverse-CDF, equal unique training area × vintage cutpoints, separately by vintage; 1 most deprived / highest no-qualifications. Preserve all ten source deciles. Fold cuts relearned in every inner/outer/bootstrap fit.',
      'models':{'B':'Existing baseline','C-linear':'Baseline + three continuous linear scores','C-smooth':'Baseline + three natural splines, 3 df each','DE':'Baseline + all ten categorical deciles for each feature','Q':'Baseline + fixed pairs 1–2,3–4,5–6,7–8,9–10 for each feature','F':'Baseline + adjacent fused decile effects; separate groupings within each feature, shared tuned penalty ratio'},
      'estimator':'NB2 jointly estimated dispersion. Minimise mean negative log likelihood plus lambda times sum of absolute adjacent coefficient differences. Reference decile 5; no monotonicity. Original tiny ridge 1e-6 on base and level coefficients for numerical stability.',
      'penalty':'lambda = ratio × maximum absolute mean social score at constrained all-fused baseline NB fit. Recalculate lambda_max in each training fit. Ratios: '+str(RATIOS)+'. Descending warm-start path; ratio one uses the exact constrained baseline endpoint (amendment 1); ratio zero uses the original unpenalised categorical NB fit.',
      'solver':'Positive/negative split of adjacent differences; bounded L-BFGS-B, maxiter 4000, ftol 1e-13, gtol 1e-7, maxls 60. Fixed second polishing attempt from first solution if projected mean-objective gradient >0.0001 or unsuccessful, with ftol 1e-15. Both diagnostics retained. No changed estimator. Adjacent differences <=1e-7 are declared fused; report threshold. Every attempted path fit is archived; only converged candidates can be selected.',
      'innerValidation':'Only outer-training rows. Sort distinct training region names and allocate cyclically into min(3,number of regions) blocks. Hold each block out of ALL inner training years, assess its latest outer-training year; other regions may contribute that year. This is geographic tuning, not a future-year inner test. At least two regions required. All preprocessing refitted inside each inner block. The 2019 deprivation edition is absent before 2020, so early nested temporal fits cannot support its cutpoints; no assessment-derived fallback.',
      'selection':'Sum inner held-out Poisson deviance across all blocks. A ratio is eligible only if every inner fit converges and its outer-training fit converges with adequate merged-band support. Select the strongest penalty whose inner deviance is within 0.5% of the lowest eligible deviance. This is a prespecified simplicity tolerance, not a statistical equivalence claim. No outer outcomes enter selection.',
      'support':'Every feature × vintage × merged outer-training band needs >=10 distinct original authorities and >=100 training collisions, positive motor mileage. Bootstrap replica IDs never inflate distinct support. Inner sparse cells are not a separate eligibility filter; their effects are assessed by inner predictions.',
      'assessment':'Same count outcomes, realised miles and rows for six models. Held-out Poisson deviance primary; O/E and mileage-weighted absolute rate error secondary. Report pooled, each year, all geographic regions, and effective social coefficients.',
      'uncertainty':{'seed':SEED,'replicates':REPLICATES,'method':'Sample nine regions with replacement, retain area histories together; same regional draw for train/test/models. Relearn all cutpoints, inner region splits, tuning and fusion for each temporal fold. Refit all comparators. Percentile 95% intervals across fully paired replicates. Only nine independent region blocks and 199 draws: exploratory tail precision.',
        'failures':'Archive all failures without replacement draws. Incomplete intervals are diagnostic, never pass the uncertainty gate. Bootstrap rank and support issues disclosed.',
        'boundaryStability':'For the 2025 temporal fold, fraction of valid full refits retaining each adjacent decile boundary; cutpoints also refitted. Report selected ratio and band-count distributions. Not a p-value or causal finding.'},
      'candidateGates':'F vs B: >=1% pooled and every annual gain, paired 95% lower >0 with all 199 complete paired draws, pooled/annual O/E within 0.90–1.10, positive pooled geographic gain and >=5/9 improving regions, overall absolute error no worse, every selected point band supported, full-rank converged point designs. Same inherited research gates; no automatic pricing deployment.',
      'replacementGates':'Only after candidate gates pass: vs DE, Q and C-smooth, upper paired 95% excess deviance <0.5%, no worse pooled absolute error or calibration distance. Linear comparison also reported. No claims based on ordinary post-selection coefficient significance.',
      'publication':'Separate v3 with original v1/v2 unchanged. Show 2025 assessment model grouping trained through 2024; coefficients referenced to decile 5, learned-band support, all comparisons, boundary stability, protocol and reproducible files.',
      'limits':['Authority groupings cannot be applied to LSOAs without separate validation. LSOA matched all-road mileage remains unavailable.','Correlated area characteristics, omitted historical joint road-class/rurality exposure and retrospective cohort selection remain limitations.','No significance-based adjacent testing. A merged pair is not proof of equal underlying risk.','No fresh confirmation: retain existing prospective 2027–2029 reserve and access audit.','Bootstrap coefficient intervals are exploratory refit summaries, not post-selection confidence guarantees.'],
      'sources':['https://epub.ub.uni-muenchen.de/10625/1/tr060.pdf','https://scikit-learn.org/stable/auto_examples/model_selection/plot_nested_cross_validation_iris.html']}
    if (OUT/'numerical-amendment.json').exists():protocol['numericalAmendment']=read(OUT/'numerical-amendment.json')
    save(OUT/'protocol.json',protocol)
    paths=[Path(__file__),ROOT/'scripts/social-deciles-v1.py',ROOT/'scripts/social-exploratory-model.py',d.OUT/'frozen-design.json',OUT/'protocol.json']+([OUT/'numerical-amendment.json'] if (OUT/'numerical-amendment.json').exists() else [])
    save(OUT/'frozen-design.json',{'id':protocol['id'],'frozenAt':d.now(),'sourceCommit':'9c870cb2f98e24fd2c864dda0084f8716ab7a857','inputs':[{'path':str(p.relative_to(ROOT)),'sha256':sha(p)} for p in paths],'software':{'numpy':np.__version__,'scipy':scipy.__version__},'temporalFoldIds':[f['id'] for f in b['folds']],'spatialFoldIds':[f['id'] for f in b['futureSpatialFolds']]})
    print('Frozen before new fits',flush=True)

def step_design(train,test,t):
    X,Z,meta=m.design(train,test,d.BASE);p=X.shape[1];names=[]
    for feature in FEATURES:
        a=np.array([d.row_band(r,feature,t) for r in train]);b=np.array([d.row_band(r,feature,t) for r in test])
        for k in range(2,11):
            X=np.column_stack([X,(a>=k).astype(float)-float(5>=k)])
            Z=np.column_stack([Z,(b>=k).astype(float)-float(5>=k)])
            names.append(feature+':boundary'+str(k-1)+'/'+str(k))
    # Level coefficients relative to decile 5: 9 levels × 9 adjacent differences.
    levels=np.array([[float(j>=k)-float(5>=k) for k in range(2,11)] for j in range(1,11) if j!=5])
    T=np.kron(np.eye(3),levels)
    meta['featureNames']+=names;meta['referenceDecile']=5
    return X,Z,p,T,meta

def fused_objective(z,X,y,offset,p,T,lam):
    q=X.shape[1]-p;delta=z[p:p+q]-z[p+q:p+2*q]
    theta=np.r_[z[:p],delta,z[-1]]
    val,grad=m.nb_objective(theta,X,y,offset)
    # Replace the tiny delta ridge with the original categorical-level ridge.
    effects=T@delta
    val+=1e-6*(effects@effects-delta@delta)/2
    grad[p:-1]+=1e-6*(T.T@effects-delta)
    val=val/len(y)+lam*np.sum(z[p:p+2*q])
    grad=grad/len(y)
    return val,np.r_[grad[:p],grad[p:-1]+lam,-grad[p:-1]+lam,grad[-1]]

def projected_gradient(x,g,bounds):
    a=g.copy()
    for i,(lo,hi) in enumerate(bounds):
        if lo is not None and x[i]<=lo+1e-9 and g[i]>0:a[i]=0
        if hi is not None and x[i]>=hi-1e-9 and g[i]<0:a[i]=0
    return float(np.max(np.abs(a)))

def effects_from_delta(delta):
    out={}
    for i,f in enumerate(FEATURES):
        a=np.r_[0.,np.cumsum(delta[i*9:(i+1)*9])];out[f]=(a-a[4]).tolist()
    return out

def grouping(delta):
    groups={}
    for i,f in enumerate(FEATURES):
        bands=[[1]]
        for k,v in enumerate(delta[i*9:(i+1)*9],2):
            if abs(v)>ZERO:bands.append([k])
            else:bands[-1].append(k)
        groups[f]=bands
    return groups

def grouped_support(train,t,groups):
    result=[]
    for f in FEATURES:
        for v in t[f]:
            for group in groups[f]:
                subset=[r for r in train if d.vintage(r,f)==v and d.row_band(r,f,t) in group]
                n=len({r.get('originalAreaId',r['areaId']) for r in subset});c=sum(r['collisions'] for r in subset);miles=sum(r['motorMiles'] for r in subset)
                result.append({'feature':f,'vintage':v,'deciles':group,'areas':n,'collisions':c,'motorMiles':miles,'supported':n>=10 and c>=100 and miles>0})
    return result

def fit_path(train,test,details=True):
    t=d.thresholds(train);X,Z,p,T,meta=step_design(train,test,t);q=27
    y=np.array([r['collisions'] for r in train],float);offset=np.log([r['motorMiles']/1e5 for r in train]);testoff=np.log([r['motorMiles']/1e5 for r in test])
    beta0,diag0=m.fit_nb(X[:,:p],y,offset)
    theta0=np.r_[beta0,np.zeros(q),np.log(diag0['alpha'])]
    _,g0=m.nb_objective(theta0,X,y,offset);lmax=float(np.max(np.abs(g0[p:-1]))/len(train))
    z=np.r_[beta0,np.zeros(2*q),theta0[-1]]
    bounds=[(None,None)]*p+[(0,None)]*(2*q)+[(math.log(1e-6),math.log(10))]
    output=[]
    for ratio in RATIOS:
        try:
            attempts=[]
            if ratio==1:
                # Exact constrained endpoint avoids reporting optimizer-roundoff boundaries.
                delta=np.zeros(q);beta=np.r_[beta0,delta];e=effects_from_delta(delta)
                pred=np.exp(testoff+Z[:,:p]@beta0)
                diag=dict(diag0,endpoint='Exact constrained baseline at lambda_max')
            elif ratio==0:
                A,B,dm=d.design(train,test,'DE',t);bd,diag=m.fit_nb(A,y,offset)
                e={f:[] for f in FEATURES};delta=[]
                for i,f in enumerate(FEATURES):
                    values=list(bd[p+i*9:p+(i+1)*9]);values.insert(4,0.);e[f]=values;delta.extend(np.diff(values))
                delta=np.array(delta);beta=np.r_[bd[:p],delta]
                pred=np.exp(testoff+B@bd)
            else:
                for ftol in [1e-13,1e-15]:
                    opt=minimize(fused_objective,z,args=(X,y,offset,p,T,ratio*lmax),jac=True,method='L-BFGS-B',bounds=bounds,options={'maxiter':4000,'ftol':ftol,'gtol':1e-7,'maxls':60})
                    pg=projected_gradient(opt.x,opt.jac,bounds);attempts.append({'success':bool(opt.success),'iterations':int(opt.nit),'projectedGradient':pg,'message':str(opt.message)})
                    z=opt.x
                    if opt.success and pg<=1e-4:break
                if not opt.success or pg>1e-4 or not np.isfinite(opt.fun):raise ValueError('Fusion optimizer did not converge: '+str(attempts))
                delta=z[p:p+q]-z[p+q:p+2*q];delta[np.abs(delta)<=ZERO]=0
                beta=np.r_[z[:p],delta];e=effects_from_delta(delta)
                pred=np.exp(testoff+Z@beta)
                diag={'alpha':float(np.exp(z[-1])),'objectiveMean':float(opt.fun),'projectedGradient':pg,'attempts':attempts}
            assert np.isfinite(pred).all() and (pred>0).all()
            groups=grouping(delta);support=grouped_support(train,t,groups)
            active=np.r_[np.ones(p,dtype=bool),np.abs(delta)>ZERO]
            # A zero coefficient at lambda=0 is not a selected boundary; same rule as all paths.
            item={'ratio':ratio,'lambda':ratio*lmax,'status':'converged','predictions':pred.tolist(),'groups':groups,'effects':e,'socialParameters':sum(len(v)-1 for v in groups.values()),'support':support,'supported':all(x['supported'] for x in support),'columns':int(active.sum()),'rank':int(np.linalg.matrix_rank(X[:,active])),'diagnostics':diag}
            if details:item.update(coefficients=beta.tolist(),transform=meta,thresholds=t)
            output.append(item)
        except (ValueError,AssertionError,np.linalg.LinAlgError) as exc:
            output.append({'ratio':ratio,'status':'failed','error':str(exc),'attempts':attempts})
    return output

def inner_blocks(train):
    regions=sorted({r['region'] for r in train});assert len(regions)>=2,'At least two regions needed for nested tuning'
    k=min(3,len(regions));latest=max(r['year'] for r in train)
    for i in range(k):
        held=regions[i::k];a=[r for r in train if r['region'] not in held];b=[r for r in train if r['region'] in held and r['year']==latest]
        assert a and b and not {r['region'] for r in a}&{r['region'] for r in b}
        yield a,b,held

def tune(train,test,details=True):
    inner=[];scores={r:0. for r in RATIOS};valid={r:True for r in RATIOS}
    for a,b,regions in inner_blocks(train):
        fits=fit_path(a,b,details=False);yy=np.array([r['collisions'] for r in b])
        summary=[]
        for fit in fits:
            ratio=fit['ratio'];ok=fit['status']=='converged';valid[ratio]&=ok
            loss=float(m.deviance(yy,fit['predictions']).sum()) if ok else None
            if ok:scores[ratio]+=loss
            summary.append({k:fit[k] for k in ['ratio','status','socialParameters','diagnostics','error'] if k in fit}|{'deviance':loss})
        inner.append({'heldRegions':regions,'assessmentYear':max(r['year'] for r in b),'trainingRows':len(a),'assessmentRows':len(b),'trainingIdSha256':hashlib.sha256('|'.join(r['id'] for r in a).encode()).hexdigest(),'fits':summary})
    outer=fit_path(train,test,details)
    eligible=[r for r in outer if r['status']=='converged' and r['supported'] and valid[r['ratio']]]
    if not eligible:raise ValueError('No converged, supported tuning candidate')
    best=min(scores[r['ratio']] for r in eligible)
    chosen=next(r for r in eligible if scores[r['ratio']]<=best*(1+TOLERANCE))
    trace=[{k:r[k] for k in ['ratio','lambda','status','socialParameters','supported','error','diagnostics'] if k in r}|{'innerDeviance':scores[r['ratio']] if valid[r['ratio']] else None,'eligible':r in eligible} for r in outer]
    return chosen,{'inner':inner,'path':trace,'selectedRatio':chosen['ratio'],'minimumEligibleInnerDeviance':best,'tolerance':TOLERANCE}

def q_design(train,test,t):
    X,Z,meta=m.design(train,test,d.BASE)
    for f in FEATURES:
        a=np.array([(d.row_band(r,f,t)-1)//2+1 for r in train]);b=np.array([(d.row_band(r,f,t)-1)//2+1 for r in test])
        for band in [1,2,4,5]:
            X=np.column_stack([X,(a==band).astype(float)]);Z=np.column_stack([Z,(b==band).astype(float)]);meta['featureNames'].append(f+':quintile'+str(band))
    return X,Z,meta

def run_fold(train,test,fold,details=True):
    assert train and test and max(r['year'] for r in train)<min(r['year'] for r in test)
    assert not {r['id'] for r in train}&{r['id'] for r in test}
    chosen,trace=tune(train,test,details)
    t=d.thresholds(train);fits=[];preds=[{k:r[k] for k in ['id','areaId','region','year','collisions','motorMiles']}|{'fold':fold,'predictions':{'F':p}} for r,p in zip(test,chosen['predictions'])]
    y=np.array([r['collisions'] for r in train],float);offset=np.log([r['motorMiles']/1e5 for r in train]);testoff=np.log([r['motorMiles']/1e5 for r in test])
    for model in MODELS[:-1]:
        try:
            X,Z,meta=q_design(train,test,t) if model=='Q' else d.design(train,test,model,t)
            beta,diag=m.fit_nb(X,y,offset);p=np.exp(testoff+Z@beta)
            assert np.isfinite(p).all() and (p>0).all()
            info={'model':model,'status':'converged','columns':X.shape[1],'rank':int(np.linalg.matrix_rank(X)),'diagnostics':diag}
            if details:info.update(transform=meta,coefficients=beta.tolist())
            for row,v in zip(preds,p):row['predictions'][model]=float(v)
        except (ValueError,AssertionError,np.linalg.LinAlgError) as exc:info={'model':model,'status':'failed','error':str(exc)}
        fits.append(info)
    if details:
        for row,r in zip(preds,test):row['deciles']={f:d.row_band(r,f,t) for f in FEATURES}
    fit={'fold':fold,'trainingYears':sorted({r['year'] for r in train}),'trainingAreas':len({r.get('originalAreaId',r['areaId']) for r in train}),'assessmentYear':test[0]['year'],'selected':{k:v for k,v in chosen.items() if k!='predictions'},'tuning':trace,'comparators':fits}
    return preds,fit

def metrics(rows):
    return {model:m.metric(rows,model) if rows and all(model in r['predictions'] for r in rows) else None for model in MODELS}
def gains(ms):
    return {b:1-ms['F']['deviance']/ms[b]['deviance'] if ms['F'] and ms[b] else None for b in MODELS[:-1]}
def describe(rows):
    return {'metrics':metrics(rows),'gain':gains(metrics(rows))}
def fit():
    b=verify();WORK.mkdir(exist_ok=True);lookup={r['id']:r for r in b['rows']}
    preds=[];fits=[]
    for fold in b['folds']+b['futureSpatialFolds']:
        path=WORK/(fold['id'].replace(' ','_')+'.json.gz')
        if path.exists():payload=uncompact(path)
        else:
            train=[lookup[k] for k in fold['train']];test=[lookup[k] for k in fold['test']]
            p,f=run_fold(train,test,fold['id']);payload={'predictions':p,'fit':f};compact(path,payload)
        preds+=payload['predictions'];fits.append(payload['fit'])
        print(fold['id'],'ratio',payload['fit']['selected']['ratio'],'groups',payload['fit']['selected']['groups'],flush=True)
    compact(OUT/'point-fits.json.xz',fits);compact(OUT/'predictions.json.gz',preds)
    temporal=[r for r in preds if r['fold'].startswith('temporal')];spatial=[r for r in preds if r['fold'].startswith('spatial')]
    save(OUT/'point-summary.json',{'temporal':describe(temporal),'annual':[{'year':y,**describe([r for r in temporal if r['year']==y])} for y in [2023,2024,2025]],'spatial':describe(spatial),'regional':[{'region':g,**describe([r for r in spatial if r['region']==g])} for g in sorted({r['region'] for r in spatial})]})

def resample(rows,draws):
    return [dict(r,originalAreaId=r['areaId'],areaId=r['areaId']+'@'+str(i),id=r['id']+'@'+str(i)) for i,g in enumerate(draws) for r in rows if r['region']==g]
def bootstrap_one(i):
    path=WORK/'bootstrap'/f'{i:04d}.json.gz'
    if path.exists():return i,'checkpoint'
    b=base();regions=sorted({r['region'] for r in b['rows']});rng=np.random.default_rng(np.random.SeedSequence([SEED,i]));draws=rng.choice(regions,len(regions),replace=True).tolist();rows=resample(b['rows'],draws)
    preds=[];fits=[];failures=[]
    for fold in b['folds']:
        a=set(fold['train']);z=set(fold['test']);train=[r for r in rows if r['id'].split('@')[0] in a];test=[r for r in rows if r['id'].split('@')[0] in z]
        try:
            p,f=run_fold(train,test,fold['id'],details=False);preds+=p;fits.append(f)
        except (AssertionError,ValueError,np.linalg.LinAlgError) as exc:failures.append({'fold':fold['id'],'error':str(exc)})
    complete=not failures and all(len(r['predictions'])==len(MODELS) for r in preds)
    compact(path,{'replicate':i,'draws':draws,'complete':complete,'failures':failures,'fits':fits,'pooled':describe(preds) if complete else None,'annual':[{'year':y,**describe([r for r in preds if r['year']==y])} for y in [2023,2024,2025]] if complete else []})
    return i,'fitted'
def bootstrap(workers):
    verify();WORK.mkdir(exist_ok=True)
    with ProcessPoolExecutor(max_workers=workers) as pool:
        for n,(i,status) in enumerate(pool.map(bootstrap_one,range(REPLICATES)),1):
            if n%10==0 or n==REPLICATES:print('Full nested regional refits',n,'/',REPLICATES,flush=True)

def interval(values):
    a=[x for x in values if x is not None and math.isfinite(x)]
    return {'lower95':float(np.quantile(a,.025)) if a else None,'upper95':float(np.quantile(a,.975)) if a else None,'validReplicates':len(a),'requiredReplicates':REPLICATES}
def report():
    verify();point=read(OUT/'point-summary.json');fits=uncompact(OUT/'point-fits.json.xz');preds=uncompact(OUT/'predictions.json.gz');boot=[uncompact(WORK/'bootstrap'/f'{i:04d}.json.gz') for i in range(REPLICATES)]
    complete=[b for b in boot if b['complete']];cis={k:interval([b['pooled']['gain'][k] for b in complete]) for k in MODELS[:-1]}
    f=point['temporal']['metrics']['F'];b=point['temporal']['metrics']['B']
    gates={'pooledGain':point['temporal']['gain']['B']>=.01,'pairedLowerBound':len(complete)==REPLICATES and cis['B']['lower95']>0,'everyAnnualGain':all(a['gain']['B']>=.01 for a in point['annual']),'pooledCalibration':.9<=f['observedExpected']<=1.1,'annualCalibration':all(.9<=a['metrics']['F']['observedExpected']<=1.1 for a in point['annual']),'spatialGain':point['spatial']['gain']['B']>0,'fiveRegions':sum(r['gain']['B']>0 for r in point['regional'])>=5,'absoluteRateError':f['exposureWeightedAbsoluteRateError']<=b['exposureWeightedAbsoluteRateError'],'trainingSupport':all(v['selected']['supported'] for v in fits),'convergedFullRankFits':all(v['selected']['rank']==v['selected']['columns'] and all(c['status']=='converged' and c['rank']==c['columns'] for c in v['comparators'] if c['model']=='B') for v in fits)}
    status='supported exploratory candidate' if all(gates.values()) else 'no incremental gain' if point['temporal']['gain']['B']<=0 else 'insufficient evidence' if len(complete)<REPLICATES else 'inconsistent'
    comparisons=[]
    for k in MODELS[:-1]:
        a=point['temporal']['metrics'][k];ci=cis[k];upper=-ci['lower95'] if ci['lower95'] is not None else None
        comparisons.append({'comparator':k,'gain':point['temporal']['gain'][k],'interval':ci,'spatialGain':point['spatial']['gain'][k],'upper95ExcessDeviance':upper,'replacementEligible':all(gates.values()) and len(complete)==REPLICATES and upper<.005 and f['exposureWeightedAbsoluteRateError']<=a['exposureWeightedAbsoluteRateError'] and abs(f['observedExpected']-1)<=abs(a['observedExpected']-1)})
    display=next(v for v in fits if v['fold']=='temporal-2025');selected=[next((v for v in b['fits'] if v['fold']=='temporal-2025'),None) for b in complete];selected=[v['selected'] for v in selected if v]
    stability={}
    for feature in FEATURES:
        stability[feature]=[{'afterDecile':j,'retainedFraction':sum(any(g[0]==j+1 for g in s['groups'][feature]) for s in selected)/len(selected) if selected else None,'validRefits':len(selected)} for j in range(1,10)]
    effect_ci={feature:[interval([s['effects'][feature][j] for s in selected]) for j in range(10)] for feature in FEATURES}
    distributions={feature:[{'bands':k,'refits':sum(len(s['groups'][feature])==k for s in selected)} for k in range(1,11)] for feature in FEATURES}
    rate_rows=[]
    for year in [2023,2024,2025]:
        for feature in FEATURES:
            for band in range(1,11):
                group=[r for r in preds if r['fold']=='temporal-'+str(year) and r['deciles'][feature]==band];miles=sum(r['motorMiles'] for r in group)
                rate_rows.append({'year':year,'feature':feature,'decile':band,'areas':len(group),'motorMiles':miles,'collisions':sum(r['collisions'] for r in group),'observedRate':1e5*sum(r['collisions'] for r in group)/miles if miles else None,**{k:1e5*sum(r['predictions'][k] for r in group)/miles if miles else None for k in MODELS}})
    result={'id':'social-deciles-fusion-v3','computedAt':d.now(),'exploratory':True,'status':status,'gates':gates,'failedGates':[k for k,v in gates.items() if not v],'models':read(OUT/'protocol.json')['models'],'comparisons':comparisons,**point,'grouping':{k:v for k,v in display.items() if k!='tuning'},'foldGroups':[{'fold':v['fold'],'ratio':v['selected']['ratio'],'groups':v['selected']['groups'],'socialParameters':v['selected']['socialParameters']} for v in fits],'stability':stability,'effectIntervals':effect_ci,'bandCountDistribution':distributions,'bootstrap':{'replicates':REPLICATES,'complete':len(complete),'regionBlocks':9,'failedReplicates':[{'replicate':b['replicate'],'failures':b['failures']} for b in boot if not b['complete']],'selectedRatioCounts':[{'ratio':r,'refits':sum(s['ratio']==r for s in selected)} for r in RATIOS],'rankDeficientSelectedFits':sum(v['selected']['rank']<v['selected']['columns'] for b in boot for v in b['fits']),'rankDeficientComparatorFits':sum(c.get('rank',0)<c.get('columns',0) for b in boot for v in b['fits'] for c in v['comparators']),'failedPathFits':sum(s['status']=='failed' for b in boot for v in b['fits'] for s in v['tuning']['path'])},'limitations':read(OUT/'protocol.json')['limits']}
    save(OUT/'findings.json',result);save(ROOT/'db/research/social-deciles-fusion-v3.json',result)
    with (OUT/'predicted-rates.csv').open('w',newline='') as fh:w=csv.DictWriter(fh,fieldnames=list(rate_rows[0]));w.writeheader();w.writerows(rate_rows)
    with (OUT/'predictions.csv').open('w',newline='') as fh:
        flat=[{k:v for k,v in r.items() if k not in ['predictions','deciles']}|r['predictions']|{f+'Decile':x for f,x in r['deciles'].items()} for r in preds];w=csv.DictWriter(fh,fieldnames=list(flat[0]));w.writeheader();w.writerows(flat)
    compact(OUT/'bootstrap-refits.json.xz',boot)
    # Self-contained modelling bundle; only small frozen inputs, no raw national archives.
    paths=[Path(__file__),ROOT/'scripts/social-deciles-v1.py',ROOT/'scripts/social-exploratory-model.py',d.OUT/'frozen-design.json',OUT/'protocol.json',OUT/'frozen-design.json',OUT/'numerical-amendment.json']
    with zipfile.ZipFile(OUT/'reproduce-v3.zip','w',zipfile.ZIP_DEFLATED,compresslevel=9) as z:
        for p in paths:z.write(p,str(p.relative_to(ROOT)))
        z.writestr('README.txt','Requires Python 3.12, numpy '+np.__version__+', scipy '+scipy.__version__+'. No network calls. Frozen v3 design provided; do not rerun prepare. Run with OPENBLAS_NUM_THREADS=1: python scripts/social-deciles-fusion-v3.py fit; python scripts/social-deciles-fusion-v3.py bootstrap --workers 4; python scripts/social-deciles-fusion-v3.py report. The v1 module is imported for formulas only. verify checks the bundled files; original v1 source archives are not required. Compare results excluding computedAt and solver platform rounding. All 199 regional draws are deterministic.\n')
    shutil.copy2(Path(__file__),OUT/Path(__file__).name)
    save(OUT/'manifest.json',{'createdAt':d.now(),'files':[{'name':p.name,'bytes':p.stat().st_size,'sha256':sha(p)} for p in sorted(OUT.iterdir()) if p.is_file() and p.name!='manifest.json']})
    print(json.dumps({'status':status,'gains':point['temporal']['gain'],'spatialGains':point['spatial']['gain'],'groups':display['selected']['groups'],'completeBootstrap':len(complete),'failedGates':result['failedGates']},indent=2),flush=True)

if __name__=='__main__':
    a=argparse.ArgumentParser();a.add_argument('action',choices=['prepare','fit','bootstrap','report']);a.add_argument('--workers',type=int,default=4);args=a.parse_args()
    globals()[args.action](args.workers) if args.action=='bootstrap' else globals()[args.action]()