#!/usr/bin/env python3
"""Phase 3 exploratory authority decile rates; no causal/adjusted findings.
Requires numpy and scipy. Run from repo root: python scripts/build-parking-phase3.py
All input paths and bytes are locked by public/research/parking/v3/input-manifest.json.
"""
import bisect,collections,csv,gzip,hashlib,io,json,math,pathlib
import numpy as np
import scipy
from scipy.stats import chi2
ROOT=pathlib.Path(__file__).resolve().parents[1]
OUT=ROOT/'public/research/parking/v3'
YEARS=list(range(2021,2026));PERIODS=['pooled']+[str(y) for y in YEARS]
BASES=['offStreet','garage'];ENDPOINTS=['all','parked'];BOOTSTRAPS=2000;SEED=20261010

def sha(b):return hashlib.sha256(b).hexdigest()
def dump(name,value): (OUT/name).write_text(json.dumps(value,separators=(',',':'),ensure_ascii=False,allow_nan=False)+'\n')
def rnd(x,n=8):return round(float(x),n) if x is not None else None

def rate(count,miles):return count/miles*100000 if miles is not None and miles>0 else None

def poisson_interval(count,miles):
 if miles is None or miles<=0:return [None,None]
 if count<0 or count!=int(count):raise ValueError('Collision counts must be non-negative integers')
 lo=0 if count==0 else .5*chi2.ppf(.025,2*count)
 hi=.5*chi2.ppf(.975,2*(count+1))
 return [rnd(rate(lo,miles)),rnd(rate(hi,miles))]

def ranks(scores):
 # Boundaries use equal authority weights and keep identical scores together.
 valid=sorted(v for v in scores.values() if v is not None)
 if not valid:return {},[]
 cuts=[valid[math.ceil(len(valid)*i/10)-1] for i in range(1,10)]
 return {k:bisect.bisect_left(cuts,v)+1 for k,v in scores.items() if v is not None},cuts

def pooled(rows,endpoint):
 if not rows:return dict(count=0,motorMiles=None,rate=None,poisson95=[None,None],unknownParking=0,parkingUnknownUpperRate=None)
 miles=sum(r['motorMiles'] for r in rows);n=sum(r['collisions' if endpoint=='all' else 'parkedCollisions'] for r in rows);unknown=sum(r['unknownCollisions'] for r in rows)
 return dict(count=n,motorMiles=miles,rate=rnd(rate(n,miles)),poisson95=poisson_interval(n,miles),unknownParking=unknown,parkingUnknownUpperRate=rnd(rate(n+unknown,miles)) if endpoint=='parked' else None)

def contrast(low,high,endpoint):
 l=pooled(low,endpoint);h=pooled(high,endpoint)
 rr=(h['count']/h['motorMiles'])/(l['count']/l['motorMiles']) if l['count'] and l['motorMiles'] and h['motorMiles'] else None
 if l['count'] and h['count'] and rr:
  se=math.sqrt(1/l['count']+1/h['count']);ci=[rr*math.exp(-1.959963984540054*se),rr*math.exp(1.959963984540054*se)]
 else:ci=[None,None]
 # Resample whole authorities within the two fixed extreme bins. Every year's
 # counts and miles for a sampled authority travel together in pooled periods.
 def cluster(rows):
  out=collections.defaultdict(lambda:[0.,0.])
  for r in rows:
   out[r['areaId']][0]+=r['collisions' if endpoint=='all' else 'parkedCollisions'];out[r['areaId']][1]+=r['motorMiles']
  return np.array([out[k] for k in sorted(out)],dtype=float)
 a,b=cluster(low),cluster(high);boot=[None,None]
 if len(a)>1 and len(b)>1:
  rng=np.random.default_rng(SEED);la=a[rng.integers(0,len(a),(BOOTSTRAPS,len(a)))].sum(axis=1);hb=b[rng.integers(0,len(b),(BOOTSTRAPS,len(b)))].sum(axis=1)
  good=(la[:,0]>0)&(la[:,1]>0)&(hb[:,1]>0)
  bs=(hb[good,0]/hb[good,1])/(la[good,0]/la[good,1]);boot=[rnd(x) for x in np.quantile(bs,[.025,.975])] if len(bs) else boot
 missing_bounds=[None,None]
 if endpoint=='parked' and l['motorMiles'] and h['motorMiles'] and l['count']:
  # Unknown-involvement sensitivity, not a statistical confidence interval.
  missing_bounds=[(h['count']/h['motorMiles'])/((l['count']+l['unknownParking'])/l['motorMiles']),((h['count']+h['unknownParking'])/h['motorMiles'])/(l['count']/l['motorMiles'])]
 return dict(ratio=rnd(rr),poisson95=[rnd(x) for x in ci],authorityBootstrap95=boot,unknownAllocationRatioRange=[rnd(x) for x in missing_bounds],lowCount=l['count'],highCount=h['count'],lowRate=l['rate'],highRate=h['rate'])

def load_csv(path):
 b=path.read_bytes();b=gzip.decompress(b) if path.suffix=='.gz' else b
 return list(csv.DictReader(io.StringIO(b.decode('utf-8-sig'))))

def build():
 manifest=json.loads((OUT/'input-manifest.json').read_text())
 for f in manifest['files']:
  b=(ROOT/f['path']).read_bytes();assert sha(b)==f['sha256'],f['path'];assert len(b)==f['bytes']
 context=json.loads((OUT/'inputs/authority-context.json').read_text());v1=json.loads((ROOT/'public/research/parking/v1/audit.json').read_text());v2=json.loads((ROOT/'public/research/parking/v2/summary.json').read_text())
 areas={r['id']:r for r in context['areas']};observations={(r['areaId'],r['year']):r for r in context['observations']};coding={(r['areaId'],r['year']):r for r in v1['coverage']}
 assert len(observations)==len(coding)==153*5
 source={s['id']:s for s in context['sources']};traffic_rows=load_csv(ROOT/'public'/source['traffic']['archivePath'].lstrip('/'));traffic={(r['local_authority_code'],int(r['year'])):r for r in traffic_rows}
 assert len(traffic)==len(traffic_rows),'Duplicate authority/year traffic keys'
 lu=load_csv(ROOT/'public'/context['lookup']['archivePath'].lstrip('/'));lu=[r for r in lu if r['LSOA21CD'].startswith('E')];lookup={r['LSOA21CD']:r['UTLA23CD'] for r in lu};assert len(lookup)==len(lu)==33755
 lu_counts=collections.Counter(lookup.values());lookup_names={r['UTLA23CD']:r['UTLA23NM'] for r in lu}
 panel=[];excluded=[];by_area=collections.defaultdict(list)
 for a in areas.values():
  for year in YEARS:
   key=(a['id'],year);o=observations[key];c=coding[key];t=traffic.get(key);reasons=list(o['quality'])
   assert c['collisions']==o['collisions'];assert c['parkedCollisions']+c['otherCollisions']+c['unknownCollisions']==c['collisions']
   if a['id'] not in lu_counts:reasons.append('No exact UTLA23 lookup code; no name matching or boundary apportionment')
   if not t or float(t['all_motor_vehicles'])<=0:reasons.append('No positive exact-code DfT authority/year mileage')
   if o['motorMiles'] is None:reasons.append('Upstream collision/mileage reconciliation did not pass')
   if not reasons:assert o['motorMiles']==float(t['all_motor_vehicles'])
   r=dict(areaId=a['id'],areaName=a['name'],region=a['region'],year=year,collisions=c['collisions'],parkedCollisions=c['parkedCollisions'],otherCollisions=c['otherCollisions'],unknownCollisions=c['unknownCollisions'],sourceMotorMiles=float(t['all_motor_vehicles']) if t else None,motorMiles=o['motorMiles'] if not reasons else None,sourceEligible=not reasons,reasons=sorted(set(reasons)))
   by_area[a['id']].append(r);panel.append(r)
 cohort=sorted(k for k,rs in by_area.items() if all(r['sourceEligible'] for r in rs));assert len(cohort)>=100
 for r in panel:
  r['included']=r['areaId'] in cohort
  if not r['included']:r['rate']=None
  else:r['rate']=rate(r['collisions'],r['motorMiles'])
 for a in areas.values():
  if a['id'] not in cohort:excluded.append(dict(id=a['id'],name=a['name'],reasons=sorted({q for r in by_area[a['id']] for q in r['reasons']}),failedYears=[r['year'] for r in by_area[a['id']] if not r['sourceEligible']],collisions=sum(r['collisions'] for r in by_area[a['id']]),lsoas=lu_counts.get(a['id'],0)))
 # Aggregate every frozen Phase 2 model by VEHICLES, never by averaging LSOA %.
 acc=collections.defaultdict(lambda:[0.,0.,0.,0]);seen=collections.Counter();model_keys=set()
 with gzip.open(ROOT/'public/research/parking/v2/all-scenarios.csv.gz','rt') as f:
  for r in csv.DictReader(f):
   lsoa=r['lsoa21cd'];assert lsoa in lookup;key='|'.join([r['mapping'],r['scenario'],r['association']]);model_keys.add(key);a=lookup[lsoa];v=float(r['modelled_vehicles']);assert v>0
   target=acc[a,key];target[0]+=v;target[1]+=v*float(r['off_street_pct'])/100;target[2]+=v*float(r['garage_use_pct'])/100;target[3]+=1;seen[lsoa]+=1
 assert len(seen)==33755 and set(seen.values())=={27} and len(model_keys)==27
 models=sorted(model_keys,key=lambda k:(not k.startswith('reference|'),k))
 features={k:{a:{'vehicles':v[0],'offStreet':100*v[1]/v[0],'garage':100*v[2]/v[0],'lsoas':v[3]} for (a,key),v in acc.items() if key==k} for k in models}
 for key,ff in features.items():
  for a,f in ff.items():assert f['lsoas']==lu_counts[a]
 results=[];sensitivity=[];memberships=[];decile_rows=[];contrast_rows=[]
 sample=[r for r in panel if r['included']]
 for basis in BASES:
  for key in models:
   mapping,scenario,association=key.split('|');scores={a:features[key][a][basis] for a in cohort};bins,cuts=ranks(scores)
   for a in cohort:memberships.append(dict(areaId=a,basis=basis,mapping=mapping,scenario=scenario,association=association,sharePct=rnd(scores[a]),modeledVehicles=rnd(features[key][a]['vehicles'],3),lsoas=lu_counts[a],decile=bins[a]))
   periods={}
   for period in PERIODS:
    rr=sample if period=='pooled' else [r for r in sample if r['year']==int(period)];outcomes={}
    for endpoint in ENDPOINTS:
     deciles=[]
     for d in range(1,11):
      members=[a for a in cohort if bins[a]==d];rs=[r for r in rr if bins[r['areaId']]==d];p=pooled(rs,endpoint);vehicles=sum(features[key][a]['vehicles'] for a in members)
      row=dict(decile=d,authorities=len(members),authorityYears=len(rs),minShare=rnd(min((scores[a] for a in members),default=None)),maxShare=rnd(max((scores[a] for a in members),default=None)),vehicleWeightedShare=rnd(sum(scores[a]*features[key][a]['vehicles'] for a in members)/vehicles) if vehicles else None,**p)
      deciles.append(row);decile_rows.append(dict(basis=basis,mapping=mapping,scenario=scenario,association=association,period=period,endpoint=endpoint,**{k:v for k,v in row.items() if k!='poisson95'},poissonLow=p['poisson95'][0],poissonHigh=p['poisson95'][1]))
     low=[r for r in rr if bins[r['areaId']]==1];high=[r for r in rr if bins[r['areaId']]==10];comparison=contrast(low,high,endpoint)
     outcomes[endpoint]=dict(deciles=deciles,contrast=comparison)
     contrast_rows.append(dict(basis=basis,mapping=mapping,scenario=scenario,association=association,period=period,endpoint=endpoint,ratio=comparison['ratio'],bootstrapLow=comparison['authorityBootstrap95'][0],bootstrapHigh=comparison['authorityBootstrap95'][1]))
    periods[period]=outcomes
   sensitivity.append(dict(basis=basis,mapping=mapping,scenario=scenario,association=association,periods={p:{e:periods[p][e]['contrast']['ratio'] for e in ENDPOINTS} for p in PERIODS}))
   if mapping=='reference':
    # Separate geography sensitivity: remove London and re-rank remaining authorities.
    outside=[a for a in cohort if areas[a]['region']!='London'];nb,nc=ranks({a:scores[a] for a in outside});nrows=[r for r in sample if r['areaId'] in outside]
    no_london={e:contrast([r for r in nrows if nb[r['areaId']]==1],[r for r in nrows if nb[r['areaId']]==10],e) for e in ENDPOINTS}
    results.append(dict(id=basis+'_'+scenario+'_'+association,basis=basis,scenario=scenario,association=association,cuts=[rnd(c) for c in cuts],members=[dict(id=a,name=areas[a]['name'],region=areas[a]['region'],sharePct=rnd(scores[a]),decile=bins[a],lsoas=lu_counts[a]) for a in cohort],periods=periods,noLondon=dict(authorities=len(outside),definition='London removed, remaining authorities re-ranked; pooled 2021–2025',outcomes=no_london)))
 coverage=[]
 for year in YEARS:
  total=[r for r in panel if r['year']==year];inc=[r for r in total if r['included']]
  coverage.append(dict(year=year,totalCollisions=sum(r['collisions'] for r in total),includedCollisions=sum(r['collisions'] for r in inc),excludedCollisions=sum(r['collisions'] for r in total if not r['included']),includedMotorMiles=sum(r['motorMiles'] for r in inc),includedAuthorities=len(cohort),excludedAuthorities=len(excluded)))
 source_register=[dict(title='Phase 1 collision and parked-involvement audit',path='/research/parking/v1/audit.json',url='https://www.gov.uk/government/statistical-data-sets/road-safety-open-data',period='2021–2025'),dict(title='Phase 2 frozen parking scenarios',path='/research/parking/v2/all-scenarios.csv.gz',url='/data-lab/parking?version=v2',period='2021 households; 2024 provision')]+[dict(title=s['title'],path=s['archivePath'],url=s['url'],period=s['period']) for s in context['sources']]+[dict(title='ONS LSOA21 → UTLA23 best-fit lookup',path=context['lookup']['archivePath'],url=context['lookup']['url'],period='LSOA 2021 / authority 2023')]
 methods=dict(scope='Exploratory English authority analysis on ALL public roads. Numerator: all police-reported injury collisions, including non-motor-only collisions; denominator: DfT all-motor vehicle miles. This is an area-level event frequency measure, not per-driver collision probability or a matched residential/motor-only cohort.',geography='LSOA21 estimates are assigned using the ONS UTLA23 best-fit lookup. Exact authority codes join to current-highway-authority collision counts and year-specific traffic. No LSOA mileage is inferred. Barnsley and Sheffield code differences remain excluded rather than name-matched.',cohort='Balanced authority cohort with valid parking lookup and collision/mileage reconciliation in every year 2021–2025. An authority failing any year is excluded from every year, avoiding changing composition. No outcome magnitude is used to select deciles.',parking='For each authority, sum modelled off-street (or inside-garage) vehicles over LSOAs and divide by total modelled vehicles. Reuse all 27 frozen Phase 2 variants at their published numeric precision. These parking shares are unvalidated proxies.',deciles='Rank the fixed included authorities by parking share with equal authority weights. Cutpoints at ceil(N*d/10); ties stay together in the lower decile. Freeze the same membership across years within each scenario. These are AUTHORITY deciles, not the Phase 2 LSOA deciles.',rates='100,000 * summed collisions / summed motor vehicle miles on the same included authority-years. Never average authority rates. Unknown parking status remains in total collisions and mileage; recorded parked involvement is a lower count, with an upper sensitivity count assuming every unresolved case involved a parked vehicle.',uncertainty='Per-decile 95% exact Poisson (Garwood) intervals condition on fixed traffic estimates and ignore spatial clustering/overdispersion. D10/D1 also has a 95% percentile authority bootstrap: 2,000 samples within the two fixed bins, with all years for an authority resampled together. Neither interval covers parking measurement error, traffic estimation error or systematic under-reporting.',contrast='D10/D1 = highest-parking decile rate divided by lowest-parking decile rate. Poisson ratio interval uses log-rate approximation. No hypothesis-selection correction or causal interpretation. The scenario envelope is not a confidence interval.',london='Additional pooled sensitivity removes London and re-ranks the remainder. This checks geographic concentration, not confounding adjustment. Urbanity, road mix, deprivation, housing and vehicle ownership still require Phase 4 controls.',time='Retrospective descriptive comparison: 2024 EHS rates are applied to 2021 household populations, then compared with 2021–2025 outcomes. Annual consistency is not out-of-time predictive validation. 2021 is pandemic-affected.',limits='STATS19 excludes damage-only collisions. A parked vehicle being involved does not establish contact, fault or causation. Missing involvement coding varies by authority/year. Residential LSOA frequency and independently observed garage-use effects remain untested.')
 quality=dict(authorities=len(cohort),authorityYears=len(sample),excludedAuthorities=len(excluded),inputLsoas=len(seen),includedLsoas=sum(lu_counts[a] for a in cohort),lookupOnlyCodes=[dict(id=a,name=lookup_names[a],lsoas=lu_counts[a]) for a in sorted(set(lu_counts)-set(areas))],unmatchedLsoaKeys=0,scenarioRows=sum(seen.values()),collisionCountReconciliationFailures=0,originalLsoaResidentialTestComplete=False,parkingUseValidated=False,adjustedModelRun=False,descriptiveAnalysisRun=True)
 payload=dict(version='v3',createdAt='2026-10-10',baseEvidenceRelease=v1['baseEvidenceRelease'],status='Exploratory authority results; unadjusted and unvalidated parking predictor',years=YEARS,periods=PERIODS,scenarios=v2['scenarios'],associations=v2['associations'],bases=[dict(id='offStreet',label='Modelled off-street vehicle share'),dict(id='garage',label='Modelled inside-garage vehicle share')],endpoints=[dict(id='all',label='All reported injury collisions'),dict(id='parked',label='Recorded parked-vehicle involvement')],quality=quality,coverage=coverage,excluded=excluded,results=results,sensitivity=sensitivity,methods=methods,sources=source_register,build=dict(script='build-parking-phase3.py',scriptSha256=sha(pathlib.Path(__file__).read_bytes()),inputManifestSha256=sha((OUT/'input-manifest.json').read_bytes()),numpyVersion=np.__version__,scipyVersion=scipy.__version__,bootstrapSamples=BOOTSTRAPS,seed=SEED))
 dump('analysis.json',payload);dump('authority-year-panel.json',panel)
 def csv_write(name,rr):
  with (OUT/name).open('w',newline='') as f:
   writer=csv.DictWriter(f,fieldnames=list(rr[0]));writer.writeheader();writer.writerows(rr)
 csv_write('authority-year-panel.csv',[{**r,'reasons':'; '.join(r['reasons'])} for r in panel]);csv_write('decile-rates.csv',decile_rows);csv_write('authority-decile-membership.csv',memberships);csv_write('scenario-contrasts.csv',contrast_rows)
 dump('data-dictionary.json',dict(version='v3',units=dict(motorMiles='Vehicle miles, not million miles; road travel within the authority, including non-residents',rate='Collision events per 100,000 motor vehicle miles',sharePct='Modelled vehicle share in percent, not observed use',ratio='D10 rate / D1 rate'),nulls='Null or blank means unavailable. Zero collisions with positive mileage is a valid zero and has a positive upper confidence limit.',grain=dict(authorityYear='One current authority code × year; excluded rows retained',decile='Parking measure × mapping × scenario × allocation × period × outcome × authority decile',membership='One authority per parking-measure/scenario combination, frozen over all years'),definitions=methods,notes='Authority bootstrap intervals apply only to the D10/D1 contrast. Decile chart intervals are conditional Poisson. Exact unrounded source miles drive calculations; display miles rounded to millions. No LSOA rate is supplied.'))
 (OUT/'build-parking-phase3.py').write_bytes(pathlib.Path(__file__).read_bytes())
 ref=next(r for r in results if r['id']=='offStreet_reference_independent')
 print(json.dumps(dict(quality=quality,reference=ref['periods']['pooled'],noLondon=ref['noLondon']),indent=2))
if __name__=='__main__':build()