#!/usr/bin/env python3
"""Frozen England LSOA21 parking scenarios. Standard library only; no network at build.
Run from repository: python scripts/build-parking-phase2.py
This is an unvalidated synthetic household model, not observed parking behaviour.
"""
import bisect, csv, gzip, hashlib, io, json, math, pathlib, re, statistics, zipfile
import xml.etree.ElementTree as ET
ROOT=pathlib.Path(__file__).resolve().parents[1]
OUT=ROOT/'public/research/parking/v2'
SCENARIOS=[
 dict(id='constrained',label='Constrained capacity',garageUse=.25,garageExtraSpaces=0,otherSpaces=1,topCars=4),
 dict(id='reference',label='Reference assumptions',garageUse=.5,garageExtraSpaces=1,otherSpaces=1,topCars=3.5),
 dict(id='ample',label='Ample capacity',garageUse=.75,garageExtraSpaces=2,otherSpaces=2,topCars=3),
]
ASSOCIATIONS=[dict(id='independent',label='Independent housing and cars'),dict(id='positive',label='More cars in higher-provision housing'),dict(id='negative',label='More cars in lower-provision housing')]
TYPES=['Detached (including bungalows)','Semi-detached (including bungalows)','Terraced (including bungalows)','Purpose-built flat','Converted/shared house','Other converted building','Commercial building','Caravan/mobile structure']
NS={'t':'urn:oasis:names:tc:opendocument:xmlns:table:1.0','x':'urn:oasis:names:tc:opendocument:xmlns:text:1.0'}
def sha(b):return hashlib.sha256(b).hexdigest()
def ods_rows(path,sheet):
 root=ET.fromstring(zipfile.ZipFile(path).read('content.xml'))
 table=next(t for t in root.findall('.//t:table',NS) if t.get('{'+NS['t']+'}name')==sheet)
 result=[]
 for row in table.findall('t:table-row',NS):
  vals=[]
  for cell in row:
   v=' '.join(''.join(p.itertext()) for p in cell.findall('.//x:p',NS))
   # Limit trailing repeated blanks, never shift a populated column.
   count=int(cell.get('{'+NS['t']+'}number-columns-repeated','1'))
   vals.extend([v]*min(count,32))
  if any(vals):result.append(vals)
 return result

def census(table):
 path=OUT/'inputs'/f'census2021-{table}.zip'
 z=zipfile.ZipFile(path);member=f'census2021-{table}-lsoa.csv'
 b=z.read(member);reader=csv.DictReader(io.StringIO(b.decode('utf-8-sig')));fields=reader.fieldnames
 result={}
 for r in reader:
  if not r['geography code'].startswith('E'):continue
  code=r['geography code'];assert code not in result
  counts=[int(r[k]) for k in fields[4:]];assert all(x>=0 for x in counts)
  assert sum(counts)==int(r[fields[3]])
  result[code]=dict(name=r['geography'],total=int(r[fields[3]]),counts=counts)
 return result,dict(member=member,memberSha256=sha(b),columns=fields)

def coupled(housing,cars,score,association):
 """Household counts with exact supplied margins; sorted coupling is a stress test."""
 total=sum(cars);hs=[total*x/sum(housing) for x in housing]
 if association=='independent':return [[h*c/total for c in cars] for h in hs]
 order=sorted(range(len(hs)),key=lambda i:(score[i],i),reverse=association=='negative')
 out=[[0.0]*len(cars) for _ in hs];remaining=list(map(float,cars));j=0
 for i in order:
  h=hs[i]
  while h>1e-8 and j<len(cars):
   take=min(h,remaining[j]);out[i][j]+=take;h-=take;remaining[j]-=take
   if remaining[j]<1e-8:j+=1
 assert max([abs(sum(row)-hs[i]) for i,row in enumerate(out)]+[abs(sum(row[j] for row in out)-cars[j]) for j in range(len(cars))])<1e-6
 return out

def measure(joint,garage,other,scenario):
 n=[0,1,2,scenario['topCars']];u=scenario['garageUse'];d=scenario['garageExtraSpaces'];o=scenario['otherSpaces']
 vehicles=off=inside=0.0
 for i,row in enumerate(joint):
  for j,h in enumerate(row):
   c=n[j];vehicles+=h*c
   # Garage use is an assumed fraction of garage-provision households with cars.
   # Usable garage: one car inside, then d additional spaces; otherwise d only.
   gcap=u*(min(c,1)+min(max(c-1,0),d))+(1-u)*min(c,d)
   off+=h*(garage[i]*gcap+other[i]*min(c,o))
   inside+=h*garage[i]*u*min(c,1)
 assert -1e-6<=inside<=off+1e-6 and off<=vehicles+1e-6
 return dict(vehicles=vehicles,off=off,garage=inside,residual=vehicles-off,
  offPct=100*off/vehicles if vehicles else None,garagePct=100*inside/vehicles if vehicles else None)

def rank_deciles(values):
 valid=sorted(v for v in values if v is not None)
 if not valid:return [None]*len(values),[]
 cuts=[valid[math.ceil(len(valid)*i/10)-1] for i in range(1,10)]
 return [bisect.bisect_left(cuts,x)+1 if x is not None else None for x in values],cuts

def write_json(path,obj):path.write_text(json.dumps(obj,separators=(',',':'),ensure_ascii=False,allow_nan=False)+'\n')
def rnd(n,d=3):return round(n,d) if n is not None else None

def build():
 OUT.mkdir(parents=True,exist_ok=True);(OUT/'groups').mkdir(exist_ok=True)
 housing,hmeta=census('ts044');cars,cmeta=census('ts045');assert housing.keys()==cars.keys()
 ehs=ods_rows(ROOT/'public/research/parking/v1/inputs/ehs-da2201.ods','2024')
 ehs={r[0]:r for r in ehs if len(r)>9 and r[1] and re.fullmatch(r'[0-9.]+',r[1])}
 def rates(k):return float(ehs[k][1])/100,float(ehs[k][2])/100
 pb=[ehs[k] for k in ['purpose built flat, low rise','purpose built flat, high rise']]
 weights=[float(r[8].replace(',','')) for r in pb]
 flat=tuple(sum(float(r[j])*w for r,w in zip(pb,weights))/sum(weights)/100 for j in [1,2])
 reference=[rates('detached'),rates('semi-detached'),rates('all terrace'),flat,rates('converted flat'),rates('converted flat'),rates('converted flat'),rates('all dwellings')]
 # These scenarios expose incompatible typologies, not a confidence interval.
 bungalow=rates('bungalow');low=[];high=[]
 for i,r in enumerate(reference):
  choices=[r,bungalow] if i<3 else [rates('purpose built flat, low rise'),rates('purpose built flat, high rise')] if i==3 else [r]
  low.append(min(choices,key=sum));high.append(max(choices,key=sum))
 # No matching EHS type for commercial/mobile dwellings: stress-test 0%-100% off street.
 for i in [6,7]:low[i]=(0,0);high[i]=(1,0)
 mappings=[('reference',reference),('lower-provision',low),('higher-provision',high)]
 model_ids=[s['id']+'_'+a['id'] for s in SCENARIOS for a in ASSOCIATIONS]
 variants=[(m,s,a) for m in mappings for s in SCENARIOS for a in ASSOCIATIONS]
 ids=sorted(housing);rows=[];results=[[] for _ in variants];qa=dict(duplicateKeys=0,unmatchedHousing=0,unmatchedCars=0,componentTotalFailures=0,householdTotalMismatchRows=0,maximumHouseholdDifference=0,zeroVehicleAreas=0)
 for code in ids:
  h,c=housing[code],cars[code];assert h['name']==c['name'];delta=h['total']-c['total']
  qa['householdTotalMismatchRows']+=int(delta!=0);qa['maximumHouseholdDifference']=max(qa['maximumHouseholdDifference'],abs(delta))
  if c['counts'][1:]==[0,0,0]:qa['zeroVehicleAreas']+=1
  provision=sum(n*sum(r) for n,r in zip(h['counts'],reference))/h['total']*100
  group=re.sub(r' \d{3}[A-Z]$','',h['name']);assert group!=h['name']
  row=dict(id=code,name=h['name'],group=group,households=c['total'],housingHouseholds=h['total'],housingCounts=h['counts'],carHouseholds=c['counts'],householdDifference=delta,
   provisionPct=rnd(provision),fallbackPct=rnd(100*(h['counts'][6]+h['counts'][7])/h['total']),observedOffStreetPct=None,observedGaragePct=None,residentialRoadKm=None,streetVehiclesPerKm=None,vehicleMiles=None,collisions=None,ratePer100k=None)
  for v,((mapping,rs),s,a) in enumerate(variants):
   joint=coupled(h['counts'],c['counts'],[sum(r) for r in rs],a['id'])
   result=measure(joint,[r[0] for r in rs],[r[1] for r in rs],s);results[v].append(result)
  rows.append(row)
 cuts=[];ranks=[]
 for variant in results:
  ranks_v,cut=rank_deciles([r['offPct'] for r in variant]);ranks.append(ranks_v);cuts.append(cut)
 for j,row in enumerate(rows):
  row['models']={model_ids[i]:[rnd(results[i][j]['offPct']),rnd(results[i][j]['garagePct']),rnd(results[i][j]['residual']),ranks[i][j],rnd(results[i][j]['vehicles'])] for i in range(9)}
  for field,key in [('offPct','offStreetRange'),('garagePct','garageRange')]:
   vals=[r[j][field] for r in results if r[j][field] is not None];row[key]=[rnd(min(vals)),rnd(max(vals))] if vals else [None,None]
  valid=[r[j] for r in ranks if r[j] is not None];row['decileRange']=[min(valid),max(valid)] if valid else [None,None]
  row['sameDecileVariants']=sum(r[j]==ranks[3][j] for r in ranks) if valid else 0
 groups=[]
 for group in sorted(set(r['group'] for r in rows)):
  rs=[r for r in rows if r['group']==group];gid=sha(group.encode())[:12];path=OUT/'groups'/f'{gid}.json';write_json(path,rs)
  groups.append(dict(id=gid,name=group,count=len(rs),path='/research/parking/v2/groups/'+path.name,sha256=sha(path.read_bytes())))
 summaries=[]
 for i,model in enumerate(model_ids):
  rr=results[i];v=sum(r['vehicles'] for r in rr);off=sum(r['off'] for r in rr);garage=sum(r['garage'] for r in rr)
  deciles=[]
  for d in range(1,11):
   selected=[j for j in range(len(rows)) if ranks[i][j]==d];vehicle=sum(rr[j]['vehicles'] for j in selected);of=sum(rr[j]['off'] for j in selected)
   deciles.append(dict(decile=d,areas=len(selected),households=sum(rows[j]['households'] for j in selected),minPct=rnd(min((rr[j]['offPct'] for j in selected),default=None)),maxPct=rnd(max((rr[j]['offPct'] for j in selected),default=None)),weightedOffStreetPct=rnd(100*of/vehicle) if vehicle else None,vehicleMiles=None,collisions=None,ratePer100k=None))
  summaries.append(dict(id=model,vehicles=rnd(v),offStreetPct=rnd(100*off/v),garagePct=rnd(100*garage/v),residualVehicles=rnd(v-off),cuts=[rnd(x,8) for x in cuts[i]],deciles=deciles))
 nts=next(r for r in ods_rows(ROOT/'public/research/parking/v1/inputs/nts0908.ods','NTS0908') if r[0].startswith('2024') and r[1]=='All areas')
 sources=[]
 defs=[('Census TS044','inputs/census2021-ts044.zip','https://www.nomisweb.co.uk/output/census/2021/census2021-ts044.zip','2021-03-21',hmeta),('Census TS045','inputs/census2021-ts045.zip','https://www.nomisweb.co.uk/output/census/2021/census2021-ts045.zip','2021-03-21',cmeta),('EHS DA2201','../v1/inputs/ehs-da2201.ods','https://assets.publishing.service.gov.uk/media/6a4d8ccb895317eadc35450c/DA2201_Parking_and_mains_gas_-_dwellings.ods','2024',dict(sheet='2024')),('NTS0908','../v1/inputs/nts0908.ods','https://assets.publishing.service.gov.uk/media/6a9ed429474b8101ece46446/nts0908.ods','2024',dict(sheet='NTS0908')),('EHS DA2203','inputs/ehs-da2203.ods','https://assets.publishing.service.gov.uk/media/6a4d8cfe92f65594cb029995/DA2203_Parking_and_mains_gas_-_households.ods','2024',dict(sheet='2024'))]
 for title,path,url,period,meta in defs:
  p=OUT/path;sources.append(dict(title=title,path='/'+str(p.resolve().relative_to(ROOT/'public')),url=url,period=period,sha256=sha(p.read_bytes()),bytes=p.stat().st_size,extracted=meta,license='Open Government Licence v3.0'))
 benchmark=dict(source='NTS0908',year=2024,garagePct=float(nts[2]),offStreetPct=float(nts[2])+float(nts[3]),streetPct=float(nts[4]),otherPct=float(nts[5]),sampleVehicles=int(nts[7].replace(',','')),usedToFit=False,localValidation=False,note='Rounded national vehicle-use shares; 2024 survey versus synthetic 2021 household population. Other locations are not all on street. No confidence intervals available in this table.')
 occupied=next(r for r in ods_rows(OUT/'inputs/ehs-da2203.ods','2024') if r[0]=='all households' and len(r)>2 and r[1])
 qa.update(lsoaCount=len(rows),housingHouseholds=sum(h['total'] for h in housing.values()),carHouseholds=sum(c['total'] for c in cars.values()),fallbackAbove5Pct=sum(r['fallbackPct']>5 for r in rows),sameDecileAll27=sum(r['sameDecileVariants']==27 for r in rows),medianDecileSpan=statistics.median(r['decileRange'][1]-r['decileRange'][0] for r in rows if r['decileRange'][0] is not None))
 methodology=dict(geography='England; Census 2021 LSOAs',reference='2021 household counts with 2024 EHS provision rates; no annual backcast',householdReconciliation='Scale TS044 type counts proportionally to each TS045 household total. Preserve both published totals and their difference.',joint='Independent within each LSOA by default. Positive/negative rank couplings preserve both margins, concentrating more cars in higher/lower-provision types; they are stress tests, not joint Census observations.',capacity='For n vehicles, garage-use probability u and d extra spaces, garage-provision off-street capacity is u[min(n,1)+min(max(n−1,0),d)]+(1−u)min(n,d). Other off-street provision accommodates min(n,o). Apply these to synthetic household counts.',residual='Total modeled household vehicles minus off-street vehicles. This is street-or-other demand, not observed street parking. No per-km measure without matched residential centreline geometry.',deciles='Rank eligible English LSOAs by unrounded modeled vehicle off-street share, equal area weights. Cutpoints at ceil(N×d/10). Ties stay together in the lower decile; counts may differ. Re-rank separately for every scenario. Deciles are exposure candidates, not collision findings.',ranges='27 variants = 3 capacity/use scenarios × 3 housing–car associations × 3 type mappings. Range is a sensitivity envelope, not a probability or confidence interval, and does not span every possible model.',mapping='House categories use EHS detached/semi/terraced rates, although Census includes bungalows. Purpose-built flats use the EHS population-weighted low/high-rise mix. Converted categories and commercial flats use converted-flat rates; mobile homes use all-dwelling rates. Alternative mappings substitute bungalow/flat extrema and 0–100% provision for commercial/mobile categories.',validation='National NTS benchmark held out from construction; numerical proximity is not validation. No aligned independent LSOA overnight vehicle-use validation, and no local error, correlation or holdout score is claimed.',interpretation='A housing/car-derived predictor cannot independently identify a parking effect beyond those ingredients. EHS rates describe provision, not capacity or use. Garage-use assumptions are not estimated from observed local behaviour.')
 payload=dict(version='v2',createdAt='2026-10-10',status='Experimental parking scenarios · independent validation outstanding',modelReady=False,analysisRun=False,phase2Complete=False,groups=groups,scenarios=SCENARIOS,associations=ASSOCIATIONS,modelOrder=model_ids,modelColumns=['offStreetPct','garageUsePct','streetOrOtherVehicles','decile','modeledVehicles'],national=summaries,quality=qa,benchmark=benchmark,provisionBenchmark=dict(source='EHS DA2203',householdOffStreetPct=float(occupied[1])+float(occupied[2]),modeledPct=rnd(sum(r['provisionPct']*r['households'] for r in rows)/sum(r['households'] for r in rows))),sources=sources,housingTypes=TYPES,mapping=[dict(type=t,garagePct=rnd(g*100),otherOffStreetPct=rnd(o*100)) for t,(g,o) in zip(TYPES,reference)],methodology=methodology,build=dict(script='build-parking-phase2.py',scriptSha256=sha(pathlib.Path(__file__).read_bytes())),validationRegister=[dict(title='Coventry · Earlsdon parking beat survey',url='https://www.coventry.gov.uk/roads-highways-pavements/earlsdon-traffic-data/7',checkedAt='2026-10-10',status='Reviewed; not comparable',detail='Hourly beats, 06:00–24:00, weekday and Saturday, selected streets. Published percentages measure occupied spaces, not resident-vehicle parking shares. No matched LSOA denominator or full overnight coverage.'),dict(title='Bristol · March 2025 parking survey package',url='https://www.bristol.gov.uk/ask/projects/south-bristol-liveable-neighbourhood/key-dates-and-reports-sbln',checkedAt='2026-10-10',status='Identified; download blocked',detail='Council describes daytime/night-time occupancy surveys. ZIP retrieval returned HTTP 403 in this run. Raw observations and exact LSOA/road matching were not verified.')])
 write_json(OUT/'summary.json',payload)
 # One-row-per-LSOA reference extract, preserving nulls as empty CSV fields.
 fields=['lsoa21cd','name','households','housing_households','household_difference','provision_proxy_pct','fallback_households_pct','modelled_vehicles','modelled_off_street_pct','modelled_garage_use_pct','street_or_other_vehicles','reference_decile','off_street_sensitivity_min','off_street_sensitivity_max','garage_sensitivity_min','garage_sensitivity_max','decile_sensitivity_min','decile_sensitivity_max','observed_off_street_pct','observed_garage_pct','residential_road_km','street_vehicles_per_km','vehicle_miles','collisions','rate_per_100000_miles']
 with (OUT/'lsoa-reference-estimates.csv').open('w',newline='') as f:
  writer=csv.writer(f);writer.writerow(fields)
  for r in rows:
   m=r['models']['reference_independent'];writer.writerow([r['id'],r['name'],r['households'],r['housingHouseholds'],r['householdDifference'],r['provisionPct'],r['fallbackPct'],m[4],m[0],m[1],m[2],m[3],*r['offStreetRange'],*r['garageRange'],*r['decileRange'],None,None,None,None,None,None,None])
 buf=io.StringIO();w=csv.writer(buf);w.writerow(['lsoa21cd','mapping','scenario','association','modelled_vehicles','off_street_pct','garage_use_pct','street_or_other_vehicles','decile'])
 for j,r in enumerate(rows):
  for v,((m,_),s,a) in enumerate(variants):
   rr=results[v][j];w.writerow([r['id'],m,s['id'],a['id'],rnd(rr['vehicles']),rnd(rr['offPct']),rnd(rr['garagePct']),rnd(rr['residual']),ranks[v][j]])
 (OUT/'all-scenarios.csv.gz').write_bytes(gzip.compress(buf.getvalue().encode(),mtime=0))
 dictionary=dict(version='v2',rowGrain='One English Census 2021 LSOA in the reference CSV',referenceModel='reference_independent with reference housing mapping',nulls='Blank CSV cells and JSON null mean unavailable; they are not zero.',inputs='Census TS044 housing and TS045 household car availability; 2021 counts',modelled_vehicles='Household cars/vans: 1*n1 + 2*n2 + 3.5*n3plus in reference scenario; not an observed exact count',modelled_off_street_pct='100 * modelled off-street vehicles / modelled household cars-vans; unvalidated',modelled_garage_use_pct='100 * modelled inside-garage vehicles / modelled household cars-vans; unvalidated',street_or_other_vehicles='Modelled household vehicles minus modelled off-street vehicles; cannot be assumed on street',provision_proxy_pct='Housing-weighted EHS garage + other off-street provision percentage; household denominator',fallback_households_pct='Commercial and mobile housing categories using proxy mappings',household_difference='Published TS044 total minus TS045 total; housing mix is scaled to TS045 total',reference_decile='Equal-LSOA-weight national rank of unrounded reference off-street share; 1 low, 10 high; ties kept together',sensitivity='Min/max over 27 explicitly defined scenarios; not confidence intervals or exhaustive bounds',unavailable=['observed_off_street_pct','observed_garage_pct','residential_road_km','street_vehicles_per_km','vehicle_miles','collisions','rate_per_100000_miles'],allScenariosRows=len(rows)*len(variants),groupJsonModelColumns=payload['modelColumns'],sourceFiles='See summary.json for source URLs, archive paths, observation dates and SHA-256 hashes',limitations=methodology)
 write_json(OUT/'data-dictionary.json',dictionary)
 (OUT/'build-parking-phase2.py').write_bytes(pathlib.Path(__file__).read_bytes())
 print(json.dumps(dict(quality=qa,national=[{k:r[k] for k in ['id','offStreetPct','garagePct']} for r in summaries],benchmark=benchmark),indent=2))
if __name__=='__main__':build()