#!/usr/bin/env python3
"""Phase 4 adjusted descriptive and spatial-transfer analysis. NumPy + SciPy.
prepare freezes protocol, sample and folds before fit; run validates every hash.
python scripts/model-parking-phase4.py prepare|run
"""
import bisect,collections,csv,gzip,hashlib,json,math,pathlib,sys,time
import numpy as np,scipy
from scipy.optimize import minimize
from scipy.special import gammaln,digamma,xlogy
from scipy.stats import norm,chi2
ROOT=pathlib.Path(__file__).resolve().parents[1];OUT=ROOT/'public/research/parking/v4';SEED=20261010
URBAN=['urbanShare','logDensity','densitySquared','trunkShare','carShare','hgvShare','vanShare','logFlow']
EXPANDED=URBAN+['income','noQualifications','ageYoung','ageOld','detached','terraced','flats','noCars','twoPlusCars']
ENDPOINTS=['all','parked','other'];STAGES={'raw':['parking'],'urban':URBAN+['parking'],'expanded':EXPANDED+['parking'],'interaction':EXPANDED+['parking','parkingUrban']}
OUTCOME={'all':'collisions','parked':'parkedCollisions','other':'otherCollisions'}
def sha(p):return hashlib.sha256(p.read_bytes()).hexdigest()
def read(p):return json.loads(p.read_text())
def dump(name,d):(OUT/name).write_text(json.dumps(d,separators=(',',':'),ensure_ascii=False,allow_nan=False)+'\n')
def csvwrite(name,rr):
 with (OUT/name).open('w',newline='') as f:w=csv.DictWriter(f,fieldnames=list(rr[0]));w.writeheader();w.writerows(rr)
def rnd(x,n=6):return round(float(x),n) if x is not None and np.isfinite(x) else None
def deviance(y,mu):return 2*(xlogy(y,y/mu)-y+mu)
def score(theta,X,y,offset):
 eta=offset+X@theta[:-1];mu=np.exp(np.clip(eta,-700,700));r=np.exp(-theta[-1]);den=r+mu
 sb=((y-mu)*r/den)[:,None]*X
 sa=-r*(digamma(y+r)-digamma(r)+np.log(r)+1-np.log(den)-(r+y)/den)
 return np.column_stack([sb,sa])
def objective(theta,X,y,offset):
 eta=offset+X@theta[:-1]
 if np.max(eta)>700:return 1e100,np.zeros_like(theta)
 mu=np.exp(eta);r=np.exp(-theta[-1]);den=r+mu
 ll=gammaln(y+r)-gammaln(r)-gammaln(y+1)+r*(np.log(r)-np.log(den))+y*(eta-np.log(den))
 grad=-score(theta,X,y,offset).sum(axis=0);grad[1:-1]+=1e-6*theta[1:-1]
 return float(-sum(ll)+.5e-6*sum(theta[1:-1]**2)),grad

def fit(X,y,offset,clusters=None):
 start=np.zeros(X.shape[1]+1);start[0]=np.log(sum(y)/sum(np.exp(offset)));start[-1]=np.log(.1)
 opt=minimize(objective,start,args=(X,y,offset),jac=True,method='L-BFGS-B',bounds=[(None,None)]*X.shape[1]+[(np.log(1e-6),np.log(10))],options={'maxiter':2000,'ftol':1e-12,'gtol':1e-6,'maxls':50})
 g=float(np.max(np.abs(opt.jac)))
 if not opt.success or g>.08 or not np.isfinite(opt.fun):raise ValueError(f'Fit failed: {opt.message}; max gradient={g}')
 cov=None
 if clusters is not None:
  t=opt.x;p=len(t);H=np.zeros((p,p))
  for k in range(p):
   h=1e-4*max(1,abs(t[k]));plus=t.copy();minus=t.copy();plus[k]+=h;minus[k]-=h;H[:,k]=(objective(plus,X,y,offset)[1]-objective(minus,X,y,offset)[1])/(2*h)
  H=(H+H.T)/2;bread=np.linalg.pinv(H,rcond=1e-10);scores=score(t,X,y,offset);blocks=np.array([scores[np.array(clusters)==a].sum(axis=0) for a in sorted(set(clusters))]);G=len(blocks);N=len(y)
  cov=bread@blocks.T@blocks@bread*(G/(G-1))*((N-1)/(N-p));cov=(cov+cov.T)/2
 diag=dict(alpha=rnd(np.exp(opt.x[-1])),gradientMax=rnd(g),iterations=int(opt.nit),conditionNumber=rnd(np.linalg.cond(X)),columns=X.shape[1],rank=int(np.linalg.matrix_rank(X)),trainingRows=len(y),converged=True)
 return opt.x,cov,diag

def design(train,test,features,context=True):
 def numeric(rows):return np.array([[r[k] for k in features] for r in rows],float)
 a,b=numeric(train),numeric(test);mean=a.mean(axis=0);scale=a.std(axis=0);scale[scale<1e-12]=1;a=(a-mean)/scale;b=(b-mean)/scale
 years=sorted({r['year'] for r in train})[1:] if context else [];regions=sorted({r['region'] for r in train}) if context else [];shares={g:sum(r['region']==g for r in train)/len(train) for g in regions}
 def cats(rr):return np.array([[float(r['year']==y) for y in years]+[float(r['region']==g) if r['region'] in regions else shares[g] for g in regions[1:]] for r in rr]).reshape(len(rr),len(years)+max(len(regions)-1,0))
 X=np.column_stack([np.ones(len(train)),a,cats(train)]);Z=np.column_stack([np.ones(len(test)),b,cats(test)])
 return X,Z,dict(names=['intercept']+features+['year:'+str(y) for y in years]+['region:'+g for g in regions[1:]],mean=mean.tolist(),scale=scale.tolist(),trainingRegions=regions,unseenRegionPolicy='Training-authority-weighted mean of training regional effects')

def effect(theta,cov,meta,urban=None):
 if 'parking' not in meta['names']:return None
 a=np.zeros(len(theta));i=meta['names'].index('parking');a[i]=.1/meta['scale'][i-1]
 if urban is not None and 'parkingUrban' in meta['names']:
  j=meta['names'].index('parkingUrban');a[j]=.1*(urban-.5)/meta['scale'][j-1]
 logrr=float(a@theta);var=float(a@cov@a);se=math.sqrt(max(0,var));return dict(ratio=rnd(np.exp(np.clip(logrr,-700,700))),ci95=[rnd(np.exp(np.clip(logrr-1.959964*se,-700,700))),rnd(np.exp(np.clip(logrr+1.959964*se,-700,700)))],changePercentagePoints=10,urbanShare=urban)

def metrics(y,mu,miles):return dict(deviance=float(sum(deviance(y,mu))),observedExpected=float(sum(y)/sum(mu)),rateError=float(1e5*sum(abs(y-mu))/sum(miles)),count=int(sum(y)),expected=float(sum(mu)))
def comparison(rows,base,pred,endpoint):
 y=np.array([r[OUTCOME[endpoint]] for r in rows]);miles=np.array([r['motorMiles'] for r in rows]);m=metrics(y,pred,miles);b=metrics(y,base,miles);blocks=[];regional=[]
 for g in sorted({r['region'] for r in rows}):
  ix=np.array([r['region']==g for r in rows]);bm=metrics(y[ix],base[ix],miles[ix]);pm=metrics(y[ix],pred[ix],miles[ix]);blocks.append([bm['deviance'],pm['deviance']]);regional.append(dict(region=g,gainPct=rnd(100*(1-pm['deviance']/bm['deviance'])),observedExpected=rnd(pm['observedExpected']),authorities=len({r['areaId'] for r in rows if r['region']==g})))
 blocks=np.array(blocks);rng=np.random.default_rng(SEED);sums=blocks[rng.integers(0,len(blocks),(2000,len(blocks)))].sum(axis=1);gains=100*(1-sums[:,1]/sums[:,0]);annual=[]
 for yr in sorted({r['year'] for r in rows}):
  ix=np.array([r['year']==yr for r in rows]);bm=metrics(y[ix],base[ix],miles[ix]);pm=metrics(y[ix],pred[ix],miles[ix]);annual.append(dict(year=yr,gainPct=rnd(100*(1-pm['deviance']/bm['deviance'])),observedExpected=rnd(pm['observedExpected'])))
 return dict(gainPct=rnd(100*(1-m['deviance']/b['deviance'])),gain95=[rnd(x) for x in np.quantile(gains,[.025,.975])],positiveRegions=sum(r['gainPct']>0 for r in regional),regions=regional,annual=annual,observedExpected=rnd(m['observedExpected']),baselineObservedExpected=rnd(b['observedExpected']),absoluteRateError=rnd(m['rateError']),baselineAbsoluteRateError=rnd(b['rateError']),regionBlocks=len(blocks))

def bins(scores,n=10):
 s=sorted(scores.values());cuts=[s[math.ceil(len(s)*i/n)-1] for i in range(1,n)];return {k:bisect.bisect_left(cuts,v)+1 for k,v in scores.items()}
def rates(rows,scores,endpoint,expected,n=10):
 bb=bins(scores,n);overall=sum(r[OUTCOME[endpoint]] for r in rows)/sum(r['motorMiles'] for r in rows)*1e5;rr=[]
 for b in range(1,n+1):
  ix=[i for i,r in enumerate(rows) if bb[r['areaId']]==b];rs=[rows[i] for i in ix];ids={r['areaId'] for r in rs};count=sum(r[OUTCOME[endpoint]] for r in rs);m=sum(r['motorMiles'] for r in rs);e=sum(expected[i] for i in ix);lo=.5*chi2.ppf(.025,2*count) if count else 0;hi=.5*chi2.ppf(.975,2*(count+1))
  rr.append(dict(bin=b,authorities=len(ids),minShare=rnd(min((scores[a] for a in ids),default=None)),maxShare=rnd(max((scores[a] for a in ids),default=None)),collisions=count,miles=m,expected=rnd(e),rawRate=rnd(count/m*1e5) if m else None,adjustedRate=rnd(count/e*overall) if e else None,raw95=[rnd(lo/m*1e5),rnd(hi/m*1e5)] if m else [None,None],adjusted95=[rnd(lo/e*overall),rnd(hi/e*overall)] if e else [None,None],unknownParking=sum(r['unknownCollisions'] for r in rs)))
 return rr

def prepare():
 assert not (OUT/'design.json').exists(),'Design is frozen; create another version to change it.'
 rows=read(OUT/'panel.json');features=read(OUT/'features.json');required=set(EXPANDED)-{'densitySquared'};ex=[];ids=[]
 for a in sorted({r['areaId'] for r in rows}):
  ar=[r for r in rows if r['areaId']==a];missing=sorted({k for r in ar for k in required if r.get(k) is None or not np.isfinite(r[k])})
  if missing:ex.append(dict(areaId=a,name=ar[0]['name'],reasons=missing))
  else:ids.append(a)
 kept=[r for r in rows if r['areaId'] in ids]
 for r in kept:r['densitySquared']=r['logDensity']**2
 folds=[dict(id=g,trainAuthorities=[a for a in ids if next(r for r in kept if r['areaId']==a)['region']!=g],testAuthorities=[a for a in ids if next(r for r in kept if r['areaId']==a)['region']==g]) for g in sorted({r['region'] for r in kept})]
 ntids=[a for a in ids if all(r['nonTrunkMiles'] and r['nonTrunkMiles']<=r['motorMiles'] and r['unknownTrunkCount']==0 for r in kept if r['areaId']==a)]
 source_manifest=read(OUT/'input-manifest.json');files=source_manifest['files']+[dict(path=str(pathlib.Path(__file__).relative_to(ROOT)),sha256=sha(pathlib.Path(__file__)),bytes=pathlib.Path(__file__).stat().st_size)]
 protocol=dict(aim='Incremental spatial predictive information of an unvalidated parking proxy; not a causal garage-use effect.',design='Exploratory, frozen after Phase 3 descriptive results were seen. No confirmatory or future-year validation claim.',models=dict(raw='Parking alone, no year/region controls; explanatory comparison',urban='Rurality, density and density squared, trunk mileage share, car/HGV/van mileage shares, log miles per road-km, year and region, plus parking',expanded='Urban model plus income deprivation, no qualifications, young/older resident shares, detached/terraced/flat shares and zero-/two-plus-car household shares',interaction='Expanded plus parking × (urban population share − 0.5)',network='Expanded plus log(1 + junctions/km²); April 2026 snapshot, retrospective sensitivity only'),family='Negative binomial NB2, maximum likelihood; log(vehicle miles / 100,000) offset. Fixed tiny ridge 1e-6 on non-intercept coefficients for numerical stability; no hyperparameter tuning.',parking='Continuous proportion; association reported per +10 percentage points. Both off-street and inside-garage shares, all 27 Phase 2 variants. Housing/car inputs also enter the baseline; a derived parking score is not independently observed behaviour.',cohort='Same complete authorities in every stage and variant, fixed over 2021–2025. No imputation. Extra exclusions listed.',spatialValidation='Leave one of nine English regions out; all years for its authorities held together. Scaling learned in training only. Unseen-region effects use the training-authority-weighted mean. Every held-out observation gets one prediction per model. This is spatial transfer using retrospectively available inputs, not future prediction.',uncertainty='Full-fit coefficient intervals use an authority-cluster sandwich including the NB dispersion parameter; fixed-shape specification and measurement error not covered. CV gain intervals resample nine region loss blocks 2,000 times: descriptive instability intervals, not independent confirmatory inference.',metrics='Percentage reduction in held-out Poisson deviance versus the same-sample baseline, plus observed/expected and exposure-weighted absolute rate error. Positive gain improves prediction; negative gain worsens it. No scenario selected as best.',adjustedCurve='Observed / expected (expanded baseline WITHOUT parking, full-sample fitted) × the cohort pooled rate. This indirect standardisation is descriptive and in-sample. Conditional Poisson intervals ignore uncertainty in expected counts; they are not causal effect intervals.',ruralGroups='Analyst-defined population bands: urban >=80% urban residents, mixed 50–<80%, rural <50%. Quintiles within groups of at least 15 authorities; smaller groups show pooled rates only. Ties kept together. Resident rurality is not road-mileage rurality.',nonTrunk='Separate reference-model sensitivity with TRA8903 mileage × 1,000,000 and canonical STATS19 trunk flag 2. Require zero unknown trunk flags across all five years per authority; compare to all-road fits on this identical subcohort. Non-trunk is not synonymous with residential.',unknown='Parked and other endpoints are recorded lower counts. Unknown involvement remains separately tabulated and is never reclassified as other. Additional reference sensitivity assigns all unresolved cases to each endpoint separately.',limits='No new LSOA mileage. Injury collisions only. Parking use unvalidated; 2021 households and 2024 parking provision are retrospectively applied. Network snapshot is 2026. No claim of year-ahead validity, causation, or driver-level demographic effects.')
 dump('design.json',dict(version='v4',createdAt='2026-10-10',protocol=protocol,files=files,cohort=ids,rows=kept,exclusions=ex,folds=folds,nonTrunkCohort=ntids,seed=SEED,software=dict(numpy=np.__version__,scipy=scipy.__version__)))
 print(json.dumps(dict(authorities=len(ids),rows=len(kept),exclusions=ex,nonTrunkAuthorities=len(ntids),folds=[dict(region=f['id'],authorities=len(f['testAuthorities'])) for f in folds]),indent=2))

class Engine:
 def __init__(self,rows):self.rows=rows;self.cache={};self.predictions=[];self.fits=[]
 def model(self,features,endpoint,context=True,full_only=False):
  key=(tuple(features),endpoint,context)
  if key in self.cache:return self.cache[key]
  rows=self.rows;y=np.array([r[OUTCOME[endpoint]] for r in rows],float);off=np.log([r['motorMiles']/1e5 for r in rows]);X,Z,meta=design(rows,rows,features,context);theta,cov,diag=fit(X,y,off,[r['areaId'] for r in rows]);full=np.exp(off+X@theta[:-1]);pred=np.zeros(len(rows));fold_diags=[]
  if not full_only:
   for g in sorted({r['region'] for r in rows}):
    ti=[i for i,r in enumerate(rows) if r['region']!=g];vi=[i for i,r in enumerate(rows) if r['region']==g];tr=[rows[i] for i in ti];te=[rows[i] for i in vi];A,B,m=design(tr,te,features,context);th,_,dg=fit(A,y[ti],off[ti]);pred[vi]=np.exp(off[vi]+B@th[:-1]);fold_diags.append(dict(region=g,**dg,transform=m))
  result=dict(theta=theta,cov=cov,meta=meta,full=full,oof=pred,diagnostic=diag,foldDiagnostics=fold_diags)
  # Parking values vary by scenario, so cache baselines only.
  if 'parking' not in features:self.cache[key]=result
  return result

def run():
 frozen=read(OUT/'design.json')
 for f in frozen['files']:assert sha(ROOT/f['path'])==f['sha256'],'Frozen bytes changed: '+f['path']
 rows=frozen['rows'];ids=set(frozen['cohort']);source=read(OUT/'source-register.json');v3=read(ROOT/'public/research/parking/v3/analysis.json');members=list(csv.DictReader((ROOT/'public/research/parking/v3/authority-decile-membership.csv').open()));scores=collections.defaultdict(dict)
 for r in members:
  if r['areaId'] in ids:scores['|'.join(r[k] for k in ['basis','mapping','scenario','association'])][r['areaId']]=float(r['sharePct'])/100
 engine=Engine(rows);results=[];sensitivity=[];fit_archive=[];prediction_rows=[];coeff_rows=[]
 def archive_fit(rid,stage,endpoint,m):
  fit_archive.append(dict(id=rid,stage=stage,endpoint=endpoint,coefficients=m['theta'].tolist(),covariance=m['cov'].tolist(),transform=m['meta'],diagnostic=m['diagnostic'],foldDiagnostics=m['foldDiagnostics']))
 def set_parking(ss):
  for r in rows:r['parking']=ss[r['areaId']];r['parkingUrban']=r['parking']*(r['urbanShare']-.5)
 baseline={}
 for endpoint in ENDPOINTS:
  baseline[endpoint]={k:engine.model(ff,endpoint) for k,ff in [('urban',URBAN),('expanded',EXPANDED)]}
  for stage,m in baseline[endpoint].items():archive_fit('baseline',stage,endpoint,m)
  print('baseline fitted',endpoint,flush=True)
 for key,ss in sorted(scores.items()):
  basis,mapping,scenario,association=key.split('|');set_parking(ss)
  for endpoint in ENDPOINTS:
   stages={};models={}
   # Every variant receives identical full-fit and leave-region-out comparisons.
   for stage in ['expanded','interaction']:
    m=engine.model(STAGES[stage],endpoint);models[stage]=m;stages[stage]=dict(effect=effect(m['theta'],m['cov'],m['meta']),spatial=comparison(rows,baseline[endpoint]['expanded']['oof'],m['oof'],endpoint),diagnostic=m['diagnostic'])
    if stage=='interaction':stages[stage]['effectsByUrbanShare']=[dict(label=label,**effect(m['theta'],m['cov'],m['meta'],u)) for label,u in [('40% urban',.4),('70% urban',.7),('95% urban',.95)]]
    coeff_rows.append(dict(basis=basis,mapping=mapping,scenario=scenario,association=association,endpoint=endpoint,stage=stage,ratio=stages[stage]['effect']['ratio'],low=stages[stage]['effect']['ci95'][0],high=stages[stage]['effect']['ci95'][1],cvGainPct=stages[stage]['spatial']['gainPct']))
    archive_fit(key,stage,endpoint,m)
   sensitivity.append(dict(basis=basis,mapping=mapping,scenario=scenario,association=association,endpoint=endpoint,ratio=stages['expanded']['effect']['ratio'],ci95=stages['expanded']['effect']['ci95'],gainPct=stages['expanded']['spatial']['gainPct'],interactionGainPct=stages['interaction']['spatial']['gainPct']))
   if mapping=='reference':
    for stage in ['raw','urban']:
     m=engine.model(STAGES[stage],endpoint,context=stage!='raw',full_only=stage=='raw');stages[stage]=dict(effect=effect(m['theta'],m['cov'],m['meta']),diagnostic=m['diagnostic'])
     if stage=='urban':stages[stage]['spatial']=comparison(rows,baseline[endpoint]['urban']['oof'],m['oof'],endpoint)
     archive_fit(key,stage,endpoint,m)
    deciles=rates(rows,ss,endpoint,baseline[endpoint]['expanded']['full']);groups=[]
    for g in ['urban','mixed','rural']:
     ix=[i for i,r in enumerate(rows) if r['ruralGroup']==g];rr=[rows[i] for i in ix];aa={r['areaId'] for r in rr};groups.append(dict(id=g,authorities=len(aa),quintiles=rates(rr,{a:ss[a] for a in aa},endpoint,baseline[endpoint]['expanded']['full'][ix],5) if len(aa)>=15 else [],pooledRate=rnd(sum(r[OUTCOME[endpoint]] for r in rr)/sum(r['motorMiles'] for r in rr)*1e5),status='Quintiles shown' if len(aa)>=15 else 'Quintiles withheld: fewer than 15 authorities'))
    results.append(dict(id=basis+'_'+scenario+'_'+association+'_'+endpoint,basis=basis,scenario=scenario,association=association,endpoint=endpoint,stages=stages,deciles=deciles,ruralGroups=groups))
    for i,r in enumerate(rows):prediction_rows.append(dict(basis=basis,scenario=scenario,association=association,endpoint=endpoint,areaId=r['areaId'],region=r['region'],year=r['year'],observed=r[OUTCOME[endpoint]],miles=r['motorMiles'],baseline=baseline[endpoint]['expanded']['oof'][i],parking=models['expanded']['oof'][i],interaction=models['interaction']['oof'][i]))
  print('variant fitted',key,flush=True)
 # A priori reference sensitivities: network context, no London, non-trunk network and unknown coding.
 extensions=[]
 for basis in ['offStreet','garage']:
  ss=scores[basis+'|reference|reference|independent'];set_parking(ss)
  for kind in ['junctions','noLondon','nonTrunk','unknownUpper']:
   rr=[dict(r) for r in rows if kind!='noLondon' or r['region']!='London']
   if kind=='nonTrunk':
    rr=[r for r in rr if r['areaId'] in frozen['nonTrunkCohort']]
    for r in rr:r.update(motorMiles=r['nonTrunkMiles'],collisions=r['nonTrunkCount'],parkedCollisions=r['nonTrunkParked'],otherCollisions=r['nonTrunkOther'],unknownCollisions=r['nonTrunkUnknownParking'])
   if kind=='unknownUpper':
    for r in rr:r['parkedCollisions']+=r['unknownCollisions'];r['otherCollisions']+=r['unknownCollisions']
   features=EXPANDED+(['logJunctionDensity'] if kind=='junctions' else []);eng=Engine(rr)
   for endpoint in ENDPOINTS:
    if kind=='unknownUpper' and endpoint=='all':continue
    if len({r['areaId'] for r in rr})<60:extensions.append(dict(basis=basis,kind=kind,endpoint=endpoint,status='Too few complete authorities',authorities=len({r['areaId'] for r in rr})));continue
    b=eng.model(features,endpoint);p=eng.model(features+['parking'],endpoint);item=dict(basis=basis,kind=kind,endpoint=endpoint,status='Fitted',authorities=len({r['areaId'] for r in rr}),effect=effect(p['theta'],p['cov'],p['meta']),spatial=comparison(rr,b['oof'],p['oof'],endpoint))
    if kind=='nonTrunk':
     ar=[dict(r) for r in rows if r['areaId'] in frozen['nonTrunkCohort']];ae=Engine(ar);ab=ae.model(features,endpoint);ap=ae.model(features+['parking'],endpoint);item['allRoadSameSample']=dict(effect=effect(ap['theta'],ap['cov'],ap['meta']),spatial=comparison(ar,ab['oof'],ap['oof'],endpoint))
    extensions.append(item);archive_fit(basis+'_'+kind,'baseline',endpoint,b);archive_fit(basis+'_'+kind,'parking',endpoint,p)
  print('extensions fitted',basis,flush=True)
 summary=dict(version='v4',createdAt='2026-10-10',status='Exploratory rural–urban adjustment and spatial prediction checks',quality=dict(authorities=len(ids),authorityYears=len(rows),phase3Authorities=140,additionalExclusions=frozen['exclusions'],nonTrunkAuthorities=len(frozen['nonTrunkCohort']),collisions=sum(r['collisions'] for r in rows),miles=sum(r['motorMiles'] for r in rows),unknownParking=sum(r['unknownCollisions'] for r in rows),ruralGroups=[dict(id=g,authorities=len({r['areaId'] for r in rows if r['ruralGroup']==g})) for g in ['urban','mixed','rural']],parkingUseValidated=False,causalEffectEstablished=False,futureValidationRun=False,networkVintage='April 2026'),protocol=frozen['protocol'],scenarios=v3['scenarios'],associations=v3['associations'],bases=v3['bases'],endpoints=[dict(id='all',label='All injury collisions'),dict(id='parked',label='Recorded parked involvement'),dict(id='other',label='Recorded no parked involvement')],results=results,sensitivity=sensitivity,extensions=extensions,controlChecks=[dict(endpoint=e,expandedVsUrban=comparison(rows,baseline[e]['urban']['oof'],baseline[e]['expanded']['oof'],e)) for e in ENDPOINTS],sources=source,build=dict(designSha256=sha(OUT/'design.json'),scriptSha256=sha(pathlib.Path(__file__)),seed=SEED,numpy=np.__version__,scipy=scipy.__version__))
 dump('analysis.json',summary);(OUT/'fits.json.gz').write_bytes(gzip.compress(json.dumps(fit_archive,separators=(',',':'),allow_nan=False).encode(),mtime=0));csvwrite('held-out-predictions.csv',prediction_rows);csvwrite('scenario-effects.csv',coeff_rows)
 csvwrite('analysis-panel.csv',[{k:r[k] for k in ['areaId','name','region','year','ruralGroup','urbanShare','largerRuralShare','smallerRuralShare','collisions','parkedCollisions','otherCollisions','unknownCollisions','motorMiles','nonTrunkMiles','nonTrunkCount','unknownTrunkCount']+EXPANDED+['logJunctionDensity'] if k!='densitySquared'} for r in rows])
 (OUT/'model-parking-phase4.py').write_bytes(pathlib.Path(__file__).read_bytes());(OUT/'prepare-parking-phase4.py').write_bytes((ROOT/'scripts/prepare-parking-phase4.py').read_bytes())
 print(json.dumps({'quality':summary['quality'],'reference':next(r for r in results if r['id']=='offStreet_reference_independent_all')['stages'],'extensions':[r for r in extensions if r['basis']=='offStreet' and r['endpoint']=='all']},indent=2))
if __name__=='__main__':prepare() if sys.argv[1]=='prepare' else run()