"""Reproduce the national named-highway fatal-crash counts from official FARS annual ZIPs.
Python 3 stdlib only. Usage: python reproduce.py INPUT_ZIP_DIR OUTPUT_DIR interstate-miles-by-state.csv
Download FARS2022NationalCSV.zip through FARS2024NationalCSV.zip from
https://static.nhtsa.gov/nhtsa/downloads/FARS/{year}/National/FARS{year}NationalCSV.zip
Same inclusion rules as the Texas analysis, applied to every state and D.C.
Counts are crash burden. Deaths per 100 miles uses miles of road, not traffic.
Large trucks: vehicle.csv BODY_TYP 60-64, 66, 67, 71, 72 or 78 (trucks FARS describes as over
10,000 lbs GVWR; buses, motor homes and code 79, weight unknown, are left out).
"""
import csv, io, json, re, zipfile, pathlib, sys, collections
YEARS=(2022,2023,2024)
# FARS writes state highways as SR-. Some states pad the number with a space (US- 7).
# Hawaii's Interstates are H-1, H-2, H-3.
# Lettered U.S. and state routes (US-31W, SR-25A, SR-99E, New York's SR-908M parkways)
# are separately signed highways, so a letter attached to the number keeps them apart.
# Interstate branches (I-35E, I-35W) stay with the parent, because an attached letter
# on an Interstate is often just a direction (I-80E for eastbound).
PREFIX=re.compile(r'^(I|US|SR)-\s*(H-?)?0*(\d+)([A-Z]?)(?=$|[^0-9])')
EXCLUDE=re.compile(r'\b(?:RAMP|FRONTAGE|FRNTG|FRONT|SERVICE|SERV|SVRD|SV|FEEDER|ACCESS|CONNECTOR|BUSINESS|BUS|BU|SPUR)\b')
LARGE_TRUCK={'60','61','62','63','64','66','67','71','72','78'}
RATE_MIN_MILES=25;RATE_MIN_DEATHS=15;SHARE_MIN_CRASHES=40
def load(z,table):
    name=next(n for n in z.namelist() if n.lower().split('/')[-1]==table)
    raw=z.read(name)
    try:text=raw.decode('utf-8-sig')
    except UnicodeDecodeError:text=raw.decode('cp1252')
    return csv.DictReader(io.StringIO(text))
def classify(r):
    road=r['TWAY_ID'].strip().upper()
    if r['ROUTE'] not in('1','2','3'):return None,None,'not_interstate_us_or_state_system'
    if EXCLUDE.search(road) or r['RELJCT2']=='20':return None,None,'explicit_ramp_frontage_business_or_connector'
    m=PREFIX.match(road)
    if not m:return None,None,'unmatched_primary_road_name'
    if {'I':'1','US':'2','SR':'3'}[m[1]]!=r['ROUTE']:return None,None,'route_system_name_disagreement'
    num=str(int(m[3]));label=None
    if m[1]!='I' and m[4]:
        num+=m[4]
        if r['STATE']=='36' and m[1]=='SR' and 900<=int(m[3])<=999:label=road[m.end():].strip(' -/') or None
    return m[1]+'-'+('H' if m[2] else '')+num,label,'included'
def run(src,out,miles_csv):
    out.mkdir(parents=True,exist_ok=True)
    with open(miles_csv,newline='') as f:mrows=list(csv.DictReader(f))
    miles={(r['state'],r['route']):float(r['miles']) for r in mrows};shared={(r['state'],r['route']):r['sharedWith'] for r in mrows}
    G={};labels=collections.defaultdict(collections.Counter);states=collections.defaultdict(lambda:collections.Counter());cls=collections.Counter();seen=set()
    for year in YEARS:
        z=zipfile.ZipFile(src/f'FARS{year}NationalCSV.zip')
        trucks={(v['STATE'],v['ST_CASE']) for v in load(z,'vehicle.csv') if v['BODY_TYP'] in LARGE_TRUCK}
        for r in load(z,'accident.csv'):
            key=(year,r['STATE'],r['ST_CASE']);assert key not in seen;seen.add(key)
            st=r['STATENAME'];f=int(r['FATALS']);assert f>0;tk=(r['STATE'],r['ST_CASE']) in trucks
            s=states[st];s['fatalCrashes']+=1;s['deaths']+=f;s['largeTruckFatalCrashes']+=tk
            route,label,reason=classify(r);cls[reason]+=1
            if not route:continue
            s['matchedFatalCrashes']+=1
            g=G.setdefault((st,route),collections.Counter());g['fatalCrashes']+=1;g['deaths']+=f;g['largeTruckFatalCrashes']+=tk;g['deathsInLargeTruckCrashes']+=f*tk;g[f'deaths{year}']+=f
            if label:labels[(st,route)][label]+=1
    rows=[]
    for (st,route),g in G.items():
        row={'state':st,'route':route,'name':labels[(st,route)].most_common(1)[0][0].title() if labels[(st,route)] else '','fatalCrashes':g['fatalCrashes'],'deaths':g['deaths'],'largeTruckFatalCrashes':g['largeTruckFatalCrashes'],'deathsInLargeTruckCrashes':g['deathsInLargeTruckCrashes'],'largeTruckShare':round(g['largeTruckFatalCrashes']/g['fatalCrashes'],3),**{f'deaths{y}':g[f'deaths{y}'] for y in YEARS},'miles':'','deathsPer100Miles':'','rateRanked':0,'sharesMilesWith':shared.get((st,route),'') if route.startswith('I-') else ''}
        mi=miles.get((st,route))
        if mi and route.startswith('I-'):
            row['miles']=mi;row['deathsPer100Miles']=round(g['deaths']/mi*100,1);row['rateRanked']=int(mi>=RATE_MIN_MILES and g['deaths']>=RATE_MIN_DEATHS)
        rows.append(row)
    rows.sort(key=lambda r:(r['state'],-r['deaths'],-r['fatalCrashes'],r['route']))
    by=collections.defaultdict(list)
    for r in rows:by[r['state']].append(r)
    state_rows=[]
    for st in sorted(by):
        rs=by[st];a=rs[0];b=rs[1] if len(rs)>1 else None
        rated=sorted([r for r in rs if r['rateRanked']],key=lambda r:-r['deathsPer100Miles'])
        tk=max(rs,key=lambda r:(r['largeTruckFatalCrashes'],r['deaths']))
        state_rows.append({'state':st,'deadliestHighway':a['route'],'name':a['name'],'deaths':a['deaths'],'fatalCrashes':a['fatalCrashes'],'largeTruckFatalCrashes':a['largeTruckFatalCrashes'],'largeTruckShare':a['largeTruckShare'],'tiedOnDeaths':int(bool(b) and b['deaths']==a['deaths']),'runnerUp':b['route'] if b else '','runnerUpDeaths':b['deaths'] if b else '','deadliestInterstatePerMile':rated[0]['route'] if rated else '','deathsPer100Miles':rated[0]['deathsPer100Miles'] if rated else '','mostLargeTruckCrashesHighway':tk['route'] if tk['largeTruckFatalCrashes'] else '','mostLargeTruckCrashes':tk['largeTruckFatalCrashes'],'statewideFatalCrashes':states[st]['fatalCrashes'],'statewideDeaths':states[st]['deaths'],'statewideLargeTruckShare':round(states[st]['largeTruckFatalCrashes']/states[st]['fatalCrashes'],3)})
    rate_rows=sorted([{k:r[k] for k in('state','route','miles','deaths','deathsPer100Miles','fatalCrashes','largeTruckFatalCrashes','sharesMilesWith')} for r in rows if r['rateRanked']],key=lambda r:-r['deathsPer100Miles'])
    corr=collections.defaultdict(lambda:collections.Counter())
    for r in rows:
        if r['route'].startswith('I-') and r['miles']!='':
            c=corr[r['route']];c['states']+=1;c['miles']+=r['miles'];c['deaths']+=r['deaths'];c['fatalCrashes']+=r['fatalCrashes'];c['largeTruckFatalCrashes']+=r['largeTruckFatalCrashes']
    corridor_rows=sorted([{'route':k,'states':c['states'],'miles':round(c['miles'],2),'deaths':c['deaths'],'fatalCrashes':c['fatalCrashes'],'largeTruckFatalCrashes':c['largeTruckFatalCrashes'],'largeTruckShare':round(c['largeTruckFatalCrashes']/c['fatalCrashes'],3),'deathsPer100Miles':round(c['deaths']/c['miles']*100,1)} for k,c in corr.items()],key=lambda r:-r['deaths'])
    share_rows=sorted([{k:r[k] for k in('state','route','fatalCrashes','largeTruckFatalCrashes','largeTruckShare','deathsInLargeTruckCrashes')} for r in rows if r['fatalCrashes']>=SHARE_MIN_CRASHES],key=lambda r:(-r['largeTruckShare'],-r['fatalCrashes']))
    us={'fatalCrashes':sum(s['fatalCrashes'] for s in states.values()),'deaths':sum(s['deaths'] for s in states.values()),'largeTruckFatalCrashes':sum(s['largeTruckFatalCrashes'] for s in states.values()),'matchedFatalCrashes':cls['included'],'routeGroups':len(rows)}
    assert sum(cls.values())==us['fatalCrashes'] and sum(r['fatalCrashes'] for r in rows)==cls['included']
    for name,data in [('route-state.csv',rows),('deadliest-highway-by-state.csv',state_rows),('interstate-rate-by-state.csv',rate_rows),('interstate-corridors.csv',corridor_rows),('large-truck-share.csv',share_rows)]:
        with (out/name).open('w',newline='') as f:
            w=csv.DictWriter(f,fieldnames=list(data[0]));w.writeheader();w.writerows(data)
    (out/'us-highways.json').write_text(json.dumps({'years':list(YEARS),'source':'NHTSA FARS','retrieved':'2026-09-18','milesSource':'FHWA Interstate Route Log and Finder List, as of January 2026','minimums':{'rateMiles':RATE_MIN_MILES,'rateDeaths':RATE_MIN_DEATHS,'shareCrashes':SHARE_MIN_CRASHES},'us':us,'classificationCounts':dict(cls),'states':state_rows,'rates':rate_rows[:25],'corridors':corridor_rows[:15],'truckShare':share_rows[:15],'topByDeaths':[{k:r[k] for k in('state','route','deaths','fatalCrashes','largeTruckFatalCrashes','largeTruckShare')} for r in sorted(rows,key=lambda r:-r['deaths'])[:15]]},separators=(',',':')))
    print(json.dumps({'us':us,'classification':dict(cls)},indent=1))
if __name__=='__main__':run(pathlib.Path(sys.argv[1]),pathlib.Path(sys.argv[2]),sys.argv[3])