#!/usr/bin/env python3
"""Phase 6 fail-closed exposure preparation. No inferred mileage from population.

audit writes the real input inventory. run requires an independently prepared manifest,
complete road network, EPSG:27700 LSOA polygons, events and traffic validation observations.
The executable geometry engine is tested with labelled synthetic fixtures only.
"""
import argparse, calendar, csv, gzip, hashlib, json, math
from collections import Counter, defaultdict
from datetime import datetime, timezone
from pathlib import Path
import shapely
from shapely.geometry import shape, Point
from shapely.strtree import STRtree
ROOT=Path(__file__).resolve().parents[1]
OUT=ROOT/'public/research/social-factors/lsoa-exposure-v1'
TRANSPORT_PREDICTORS={'road_class','speed_limit','lanes','junction_density','road_length','network_connectivity','road_function','count_history','vehicle_mix'}
STRATA=[c+'|'+u for c in ['motorway','A','B','C','unclassified'] for u in ['urban','rural']]
def save(path,x):path.parent.mkdir(parents=True,exist_ok=True);path.write_text(json.dumps(x,indent=2,allow_nan=False)+'\n')
def read(path):return json.loads(path.read_text())
def sha(path):return hashlib.sha256(path.read_bytes()).hexdigest()
def stamp():return datetime.now(timezone.utc).isoformat()
def protocol():
    return {'id':'social-lsoa-exposure-2026-10-10-v1','geography':'England LSOA21; a new scale, separate from authority coefficients',
      'status':'Engineering protocol; no LSOA outcome model fitted','crs':'EPSG:27700 metres; full-resolution compatible boundaries required for allocation. Display centroids/generalised maps cannot allocate road mileage.',
      'tolerances':{s:{'maximumAbsoluteBias':.10 if s.split('|')[0] in ['motorway','A'] else .15,'maximumWeightedAbsoluteError':.20 if s.split('|')[0] in ['motorway','A'] else .30,
        'minimumIndependentTrafficObservations':30,'minimumHeldOutRegions':3,'minimum95IntervalCoverage':.90,'maximumAuthorityReconciliationDifference':.05 if s.split('|')[0] in ['motorway','A'] else .10} for s in STRATA},
      'thresholdRationale':'Prospective research acceptance tolerances, not DfT standards. Stratum-specific transport checks are frozen before LSOA outcome analysis. No outcome-informed relaxation. Domain review can propose a new version before outcome access.',
      'networkGate':'Every intended physical segment/year must appear exactly once and be fully assigned, with known measured/estimated AADT and uncertainty bounds. Missing and zero traffic are different. No arbitrary aggregate coverage threshold.',
      'mileage':'Two-way all-motor AADT × segment length in metres / 1609.344 × actual calendar days. Input adapter must consolidate directional counts; no extra factor of two.',
      'collisionGate':'One event ID/year, unique polygon assignment. Record border/near-border, positional and unmatched events. Ambiguous events are unassigned and block release. A major-road study additionally needs a validated unique event-to-covered-link match.',
      'trafficModelPredictorsAllowed':sorted(TRANSPORT_PREDICTORS),'prohibited':'Social features, deciles, deprivation, education, crime and collision outcomes must not enter the primary mileage model.',
      'independence':'Validation sites and entire holdout regions must be disjoint from traffic-model training observations and regions; evaluate weighted absolute error, aggregate bias and interval coverage by road stratum and year.',
      'reconciliation':'Compare original segment sums to official authority × year × road-class × rurality totals. Optional documented stratum factor preserves original estimates and intervals. Agreement after scaling is not independent validation.',
      'exposureUncertainty':'Carry lower/upper AADT into mileage intervals. Require >=95% of reconciled mileage within relative half-width 25%; then repeat separate LSOA collision models across exposure scenarios. No predictive release until that sensitivity and Phases 3–5 pass.',
      'toleranceMetres':{'allocationAbsolute':.01,'nearBorder':10},'humanReviewRecorded':False}
def gate(checks):return bool(checks) and all(v is True for v in checks.values())
def valid_number(x):return isinstance(x,(int,float)) and math.isfinite(x)
def geometry_index(boundaries):
    ids=[];polys=[];seen=set()
    for f in boundaries:
        p=f['properties'];code=p['lsoa21'];g=shape(f['geometry'])
        if code in seen or not code.startswith('E010') or len(code)!=9:raise ValueError('Invalid/duplicate LSOA code')
        if not g.is_valid or g.is_empty or g.geom_type not in ['Polygon','MultiPolygon']:raise ValueError('Invalid LSOA polygon')
        seen.add(code);ids.append(code);polys.append(g)
    if not ids:raise ValueError('No LSOA polygons')
    return ids,polys,STRtree(polys)
def split_segments(boundaries,segments,years):
    ids,polys,tree=geometry_index(boundaries);parts=[];audit=[];seen=set();geometries=[];geometry_keys=[]
    for f in segments:
        p=f['properties'];year=p['year'];key=(p['physical_link_id'],year)
        if key in seen:raise ValueError('Duplicate physical road link/year (possible duplicated direction)')
        seen.add(key)
        if year not in years:raise ValueError('Road year outside protocol')
        if p['flow_basis']!='both_directions':raise ValueError('Consolidate directional flow before allocation')
        if p['stratum'] not in STRATA:raise ValueError('Unknown road/rurality stratum')
        if p['traffic_status'] not in ['counted','estimated','missing']:raise ValueError('Unknown traffic provenance')
        line=shape(f['geometry'])
        if not line.is_valid or line.is_empty or line.geom_type not in ['LineString','MultiLineString'] or line.length<=0:raise ValueError('Invalid road geometry')
        geometries.append(line);geometry_keys.append(key)
        length=0.;n=0
        for ix in tree.query(line,predicate='intersects'):
            piece=line.intersection(polys[ix]);metres=piece.length
            if metres<=1e-8:continue
            length+=metres;n+=1;flow=p.get('aadt');lo=p.get('aadt_lower95');hi=p.get('aadt_upper95')
            valid=all(valid_number(x) for x in [flow,lo,hi]) and 0<=lo<=flow<=hi and p['traffic_status']!='missing'
            factor=metres/1609.344*(366 if calendar.isleap(year) else 365)
            parts.append({'physicalLinkId':p['physical_link_id'],'lsoa21':ids[ix],'authority':p['authority'],'year':year,'stratum':p['stratum'],
                'metres':metres,'trafficStatus':p['traffic_status'],'motorMilesOriginal':flow*factor if valid else None,
                'lowerMiles':lo*factor if valid else None,'upperMiles':hi*factor if valid else None})
        audit.append({'physicalLinkId':p['physical_link_id'],'year':year,'sourceMetres':line.length,'allocatedMetres':length,
            'pieces':n,'exactAccounting':abs(length-line.length)<=max(.01,line.length*1e-8)})
    # Detect positive-length overlaps even where duplicate IDs were renamed. Crossings at a point are allowed.
    tree=STRtree(geometries);overlaps=[]
    for i,line in enumerate(geometries):
        for j in tree.query(line,predicate='intersects'):
            if j>i and geometry_keys[i][1]==geometry_keys[j][1] and line.intersection(geometries[j]).length>.01:
                overlaps.append([list(geometry_keys[i]),list(geometry_keys[j])])
    return parts,audit,overlaps
def assign_collisions(boundaries,events,years,scope='all-road',covered_links=None):
    ids,polys,tree=geometry_index(boundaries);rows=[];seen=set()
    for event in events:
        key=(event['event_id'],event['year'])
        if key in seen:raise ValueError('Duplicate collision event/year')
        seen.add(key)
        if event['year'] not in years:raise ValueError('Collision year outside protocol')
        row=dict(event,lsoa21=None,assignment='missing-position',nearBorder=None)
        if valid_number(event.get('easting')) and valid_number(event.get('northing')):
            point=Point(event['easting'],event['northing']);matches=list(tree.query(point,predicate='intersects'))
            row['assignment']='unmatched' if not matches else 'ambiguous-border' if len(matches)>1 else 'assigned'
            if len(matches)==1:
                i=matches[0];row['lsoa21']=ids[i];row['nearBorder']=point.distance(polys[i].boundary)<=10
                if scope=='major-road-only' and (event.get('physical_link_id'),event['year']) not in (covered_links or set()):
                    row.update(lsoa21=None,assignment='unmatched-covered-network')
        rows.append(row)
    return rows
def traffic_checks(observations,manifest):
    predictors=set(manifest['trafficModel']['predictors']);allowed=bool(predictors) and predictors<=TRANSPORT_PREDICTORS
    train_sites=set(manifest['trafficModel']['trainingSites']);train_regions=set(manifest['trafficModel']['trainingRegions']);rows=[];seen=set()
    for o in observations:
        key=(o['site'],o['year'])
        if key in seen:raise ValueError('Duplicate traffic validation observation')
        seen.add(key)
        if o['stratum'] not in STRATA or o['year'] not in manifest['years']:raise ValueError('Traffic validation outside intended scope')
        if not all(valid_number(o[k]) for k in ['observed','predicted','lower95','upper95']) or not 0<=o['lower95']<=o['predicted']<=o['upper95'] or o['observed']<0:raise ValueError('Invalid traffic observation')
    independent=bool(observations) and bool(train_sites) and bool(train_regions) and all(o['site'] not in train_sites and o['region'] not in train_regions for o in observations)
    for year in manifest['years']:
        for s in manifest['intendedStrata']:
            group=[o for o in observations if o['year']==year and o['stratum']==s];total=sum(o['observed'] for o in group);p=protocol()['tolerances'][s]
            bias=(sum(o['predicted'] for o in group)/total-1) if total else None
            error=sum(abs(o['predicted']-o['observed']) for o in group)/total if total else None
            coverage=sum(o['lower95']<=o['observed']<=o['upper95'] for o in group)/len(group) if group else None
            passed=independent and allowed and len({o['site'] for o in group})>=p['minimumIndependentTrafficObservations'] and len({o['region'] for o in group})>=p['minimumHeldOutRegions'] and total>0 and abs(bias)<=p['maximumAbsoluteBias'] and error<=p['maximumWeightedAbsoluteError'] and coverage>=p['minimum95IntervalCoverage']
            rows.append({'year':year,'stratum':s,'observations':len(group),'bias':bias,'weightedAbsoluteError':error,'intervalCoverage':coverage,'passed':passed})
    return rows,independent,allowed
def run(manifest_path,out):
    manifest=read(manifest_path);base=manifest_path.parent
    if manifest['crs']!='EPSG:27700' or manifest['boundaryVersion']!='LSOA21':raise ValueError('Incompatible coordinate system or boundaries')
    if not manifest['years'] or not manifest['intendedStrata'] or not set(manifest['intendedStrata'])<=set(STRATA):raise ValueError('Missing intended years/strata')
    if manifest['scope'] not in ['all-road','major-road-only']:raise ValueError('Unknown study scope')
    if manifest['scope']=='all-road' and set(manifest['intendedStrata'])!=set(STRATA):raise ValueError('All-road scope requires every road/rurality stratum')
    if manifest['scope']=='major-road-only' and any(s.split('|')[0] not in ['motorway','A'] for s in manifest['intendedStrata']):raise ValueError('Major-road scope includes a minor road')
    inputs={}
    for name,record in manifest['files'].items():
        path=base/record['path']
        if sha(path)!=record['sha256']:raise ValueError('Input hash mismatch: '+name)
        inputs[name]=read(path)
    boundaries=inputs['boundaries']['features'];roads=inputs['segments']['features']
    parts,allocation,overlaps=split_segments(boundaries,roads,manifest['years'])
    links={(r['properties']['physical_link_id'],r['properties']['year']) for r in roads}
    expected={(r['physical_link_id'],r['year']) for r in inputs['networkInventory']}
    inventory_unique=len(expected)==len(inputs['networkInventory'])
    events=assign_collisions(boundaries,inputs['collisions'],manifest['years'],manifest['scope'],links)
    validation,independent,predictors=traffic_checks(inputs['trafficValidation'],manifest)
    totals=defaultdict(float);complete=defaultdict(lambda:True)
    for p in parts:
        key=(p['authority'],p['year'],p['stratum']);totals[key]+=p['motorMilesOriginal'] or 0;complete[key]&=p['motorMilesOriginal'] is not None
    reconciliations=[];factors={};official=set()
    for r in inputs['officialTotals']:
        key=(r['authority'],r['year'],r['stratum'])
        if key in official:raise ValueError('Duplicate official authority/stratum/year')
        official.add(key);raw=totals.get(key);target=r['motorMiles']
        if not valid_number(target) or target<0:raise ValueError('Invalid official miles')
        difference=raw/target-1 if raw is not None and target>0 else 0 if raw==target==0 else None
        factor=target/raw if manifest.get('reconcile') and raw and r.get('explanation') else 1.
        factors[key]=factor
        passed=difference is not None and abs(difference)<=protocol()['tolerances'][r['stratum']]['maximumAuthorityReconciliationDifference'] and bool(r.get('explanation')) and complete[key]
        reconciliations.append(dict(r,originalMotorMiles=raw,relativeDifference=difference,scalingFactor=factor,passed=passed))
    panel=defaultdict(lambda:{'motorMilesOriginal':0.,'motorMiles':0.,'lowerMiles':0.,'upperMiles':0.,'missingSegments':0,'collisions':0})
    for code in [f['properties']['lsoa21'] for f in boundaries]:
        for year in manifest['years']:panel[code,year]
    for p in parts:
        row=panel[p['lsoa21'],p['year']];factor=factors.get((p['authority'],p['year'],p['stratum']),1)
        if p['motorMilesOriginal'] is None:row['missingSegments']+=1
        else:
            row['motorMilesOriginal']+=p['motorMilesOriginal'];row['motorMiles']+=p['motorMilesOriginal']*factor
            row['lowerMiles']+=p['lowerMiles']*factor;row['upperMiles']+=p['upperMiles']*factor
    for e in events:
        if e['lsoa21']:panel[e['lsoa21'],e['year']]['collisions']+=1
    feature_ids={r['lsoa21'] for r in read(ROOT/'public/research/social-factors/deciles-v1/lsoa-profiles.json')}
    checks={'realData':manifest.get('synthetic') is False,'completeInventory':bool(expected) and inventory_unique and links==expected and bool(manifest.get('inventoryCompletenessEvidence')),
        'fullNetworkAllocation':bool(allocation) and all(a['exactAccounting'] for a in allocation) and not overlaps,
        'knownTrafficAndUncertainty':bool(parts) and all(p['motorMilesOriginal'] is not None for p in parts),
        'allEventsAssignedOnce':len(events)==manifest['expectedCollisionCount'] and all(e['assignment']=='assigned' for e in events),
        'positionalUncertaintyAudited':bool(manifest.get('positionalUncertaintyEvidence')),
        'independentTrafficValidation':independent and bool(validation) and all(v['passed'] for v in validation),'transportPredictorsOnly':predictors,
        'officialTotalsComplete':bool(reconciliations) and set(totals)==official and all(r['passed'] for r in reconciliations),
        'featureGeographyMatch':set(k[0] for k in panel)<=feature_ids,'sourceProvenance':all(r.get('publisher') and r.get('edition') and r.get('url') for r in manifest['files'].values())}
    miles=sum(p['motorMiles'] for p in panel.values());stable=sum(p['motorMiles'] for p in panel.values() if p['motorMiles']>0 and (p['upperMiles']-p['lowerMiles'])/(2*p['motorMiles'])<=.25)
    checks['exposureUncertaintyBounded']=bool(miles) and stable/miles>=.95
    passed=gate(checks)
    rows=[dict(lsoa21=code,year=year,**r,eligible=passed and r['motorMiles']>0 and r['missingSegments']==0,
               collisionRatePer100k=1e5*r['collisions']/r['motorMiles'] if passed and r['motorMiles']>0 and r['missingSegments']==0 else None) for (code,year),r in sorted(panel.items())]
    for row in rows:
        row['knownMotorMiles']=row['motorMiles']
        if row['missingSegments']:
            for field in ['motorMilesOriginal','motorMiles','lowerMiles','upperMiles']:row[field]=None
    save(out/'exposure-audit.json',{'checks':checks,'exposureGatePassed':passed,'predictionLayerAvailable':False,'remaining':'A new LSOA protocol and Phases 3–5 are still required after exposure acceptance.',
        'scope':manifest['scope'],'inputManifestSha256':sha(manifest_path),'collisions':dict(Counter(e['assignment'] for e in events)),'nearBorder':sum(e['nearBorder'] is True for e in events),
        'duplicateGeometries':overlaps,'stableMileageShare':stable/miles if miles else None,'trafficValidation':validation,'reconciliation':reconciliations})
    for name,data in [('segment-parts',parts),('allocation-audit',allocation),('collision-assignments',events),('candidate-panel',rows)]:save(out/(name+'.json'),data)
    return passed
def audit():
    OUT.mkdir(parents=True,exist_ok=True)
    p=protocol();p['registeredAt']=stamp();save(OUT/'protocol.json',p)
    inventory=[
      ('Social profiles','available','33,755 exact-code LSOA21 profiles; 2021 education and 2025 deprivation. Descriptive only.','public/research/social-factors/deciles-v1/feature-audit.json'),
      ('Collision records','partial','Archived STATS19 records contain positions and published LSOA codes; unique event-to-LSOA21 allocation and positional audit remain to be run with compatible boundaries.','public/research/parking/v1/lsoa-coding-coverage.csv'),
      ('LSOA21 allocation boundaries','missing','Full compatible polygons not loaded into this exposure pipeline; population centroids are display-only.',None),
      ('Historical all-road segments','missing','The 2026 road-network snapshot used elsewhere is not a year-matched 2018–2025 network.',None),
      ('Major-road flow/link matching','partial','DfT publishes major-road link estimates; no audited link-to-LSOA mileage allocation is present.',None),
      ('Minor-road AADT','missing','DfT count points cover a sample. Unsampled-road estimates and independent holdout validation are absent.',None),
      ('Road-class × rurality totals','missing','No matched authority/year/stratum reconciliation set is present.',None),
      ('Exposure uncertainty','missing','No segment-level uncertainty bounds or independent traffic validation set is present.',None)]
    rows=[{'input':n,'status':s,'detail':detail,'artifact':path,'sha256':sha(ROOT/path) if path and (ROOT/path).exists() else None} for n,s,detail,path in inventory]
    result={'id':p['id'],'auditedAt':stamp(),'exposureGatePassed':False,'predictionLayerAvailable':False,'phase6Status':'Pipeline implemented and fixture-tested; real LSOA frequency modelling blocked by exposure',
        'profileRows':33755,'eligibleFrequencyRows':0,'inventory':rows,'gates':[{'name':n,'passed':False,'reason':r} for n,r in [
            ('Complete intended-road network','No matched complete historical segment inventory'),('Independent traffic validation','No validated unsampled minor-road allocation'),('Authority/stratum reconciliation','No compatible year × road class × rurality totals'),('Collision allocation','Compatible polygon and positional audit pending'),('Exposure uncertainty','No validated segment uncertainty scenarios')]],
        'assumptionsRejected':['Population as vehicle-miles','Zero exposure for missing flow','Count-point totals as a complete minor-road network','Redistributing authority miles by social bands','Backdating the 2026 road network'],
        'sources':[{'publisher':'DfT','url':'https://roadtraffic.dft.gov.uk/count-points','checkedAt':stamp(),'finding':'Major-road links and a sample of minor-road locations.'},
                   {'publisher':'DfT','url':'https://roadtraffic.dft.gov.uk/downloads','checkedAt':stamp(),'finding':'Individual link/small-area estimates have greater uncertainty and may use older counts.'}],
        'nextAction':'Acquire a year-matched network and transport-only minor-road traffic model with independent holdouts, then run this exposure gate before opening a separate LSOA outcome study.',
        'majorRoadOption':'Possible only as a separate version after matching collision events and vehicle-miles to exactly the covered major-road links; not implemented as an all-road proxy.',
        'software':{'shapely':shapely.__version__},'codeSha256':sha(Path(__file__))}
    save(OUT/'audit.json',result);save(ROOT/'db/research/social-lsoa-exposure-v1.json',result)
    print(result['phase6Status'])
if __name__=='__main__':
    parser=argparse.ArgumentParser();parser.add_argument('action',choices=['audit','run']);parser.add_argument('--manifest',type=Path);parser.add_argument('--out',type=Path);args=parser.parse_args()
    if args.action=='audit':audit()
    else:
        if not args.manifest or not args.out:parser.error('run requires --manifest and --out')
        if args.out.exists():parser.error('Output exists: use a new candidate directory')
        print('Exposure gate:',run(args.manifest,args.out))