#!/usr/bin/env python3
"""Reproduce the published mathematical audit and all derived numerical tables.

Python 3.10+; standard library only; no network, third-party packages or registry access.
Run: python check-digit-audit.py --out ./reproduced --verify-serials
Input: reference-inputs.json in this script's directory (or pass --inputs).
The algorithm is the ISO-format procedure described in GDV Container Handbook 3.3,
including the crucial final mapping of raw remainder 10 to printed check digit 0.
The published article distinguishes that public explanation from licensed ISO text.

The substitution audit counts directed character-change PATTERNS, not observed errors.
The checksum-state grid separately tests each pattern against all 11 possible starting
raw remainders. Its detection percentage is an equally weighted mathematical grid,
not a measured rate, not a fleet-weighted rate and not an OCR performance estimate.
The full-space count uses an exact integer residue convolution: it does not iterate
17.576 billion strings, does not sample them and does not assume random fleet prefixes.
"""
from __future__ import annotations
import argparse
import csv
from collections import Counter
from decimal import Decimal, ROUND_HALF_UP
from itertools import combinations, product
import json
from pathlib import Path
import re
import string
from typing import Iterable

LETTERS=string.ascii_uppercase
DIGITS=string.digits
EQUIPMENT='UJZ'
WEIGHTS=tuple(1 << p for p in range(10))
VALUES:dict[str,int]={}
v=10
for letter in LETTERS:
    while v%11==0:
        v+=1
    VALUES[letter]=v
    v+=1
VALUES.update({c:int(c) for c in DIGITS})
BODY_RE=re.compile(r'[A-Z]{3}[UJZ][0-9]{6}\Z')


def raw_remainder(body:str)->int:
    """Validate structure and return the weighted sum modulo 11 (0..10)."""
    if not BODY_RE.fullmatch(body):
        raise ValueError('Expected 3 uppercase owner letters, U/J/Z, and 6 digits')
    return sum(VALUES[c]*w for c,w in zip(body,WEIGHTS))%11


def check_digit(body:str)->int:
    """A remainder of 10 prints as 0. Never discard that final step."""
    return raw_remainder(body)%10


def test_patterns(deltas:Iterable[int])->dict[str,int|float]:
    ds=list(deltas)
    same=sum(d%11==0 for d in ds)
    zero_alias=sum(d%11 in (1,10) for d in ds)
    undetected=sum(r%10==((r+d)%11)%10 for d in ds for r in range(11))
    states=11*len(ds)
    return {'directed_change_patterns':len(ds),'always_invisible_same_remainder_patterns':same,
            'additional_zero_alias_patterns':zero_alias,'checksum_states_tested':states,
            'undetected_checksum_states':undetected,
            'detected_checksum_states':states-undetected,
            'detection_pct_equal_weight_checksum_states':round(100*(states-undetected)/states,9) if states else 0.0}


def position_deltas(p:int)->list[int]:
    chars=LETTERS if p<3 else EQUIPMENT if p==3 else DIGITS
    return [(VALUES[b]-VALUES[a])*WEIGHTS[p] for a in chars for b in chars if a!=b]


def write_csv(out:Path,name:str,rows:list[dict])->None:
    if not rows:
        raise ValueError(f'Empty output {name}')
    with (out/name).open('w',encoding='utf-8',newline='') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)


def main()->None:
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--inputs',type=Path,default=Path(__file__).with_name('reference-inputs.json'))
    parser.add_argument('--out',type=Path,default=Path('reproduced'))
    parser.add_argument('--verify-serials',action='store_true',help='Independently enumerate all one million six-digit serials and compare with the exact residue convolution')
    args=parser.parse_args()
    inputs=json.loads(args.inputs.read_text(encoding='utf-8'))
    out=args.out;out.mkdir(parents=True,exist_ok=True)
    date=inputs['verified_date'];sm={s['source_id']:s for s in inputs['sources']}
    groups={r:''.join(c for c in LETTERS if VALUES[c]%11==r) for r in range(11)}
    letter_rows=[{'letter':c,'assigned_value':VALUES[c],'value_mod_11':VALUES[c]%11,'behaviour_group_members':groups[VALUES[c]%11],'source_id':'S08','verified_date':date} for c in LETTERS]
    pair_rows=[]
    for a,b in combinations(LETTERS,2):
        if VALUES[a]%11==VALUES[b]%11:
            pair_rows.append({'letter_a':a,'value_a':VALUES[a],'letter_b':b,'value_b':VALUES[b],'value_difference':VALUES[b]-VALUES[a],'shared_mod_11':VALUES[a]%11,'effect':'Always preserves the raw remainder and printed digit at a fixed owner-letter position','verified_date':date})
    assert len(pair_rows)==22
    subs=[{'position':p+1,'field':'owner' if p<3 else 'equipment' if p==3 else 'serial',**test_patterns(position_deltas(p)),'verified_date':date} for p in range(10)]
    aggregates=[]
    for name,positions in [('owner_all',range(3)),('equipment',range(3,4)),('serial_all',range(4,10)),('body_all',range(10))]:
        aggregates.append({'scope':name,**test_patterns(d for p in positions for d in position_deltas(p)),'verified_date':date})
    trans=[]
    for p in range(9):
        chars=LETTERS if p<2 else EQUIPMENT if p==2 else '' if p==3 else DIGITS
        ds=[(VALUES[b]-VALUES[a])*(WEIGHTS[p]-WEIGHTS[p+1]) for a in chars for b in chars if a!=b]
        trans.append({'position_pair':f'{p+1}-{p+2}',**test_patterns(ds),'format_scope':'All swaps cross letter/digit fields and are structurally invalid; excluded from numerical count' if not chars else ('Only U/J/Z in the third owner position keep the fourth position structurally valid' if p==2 else 'Distinct characters; unchanged same-character swaps excluded'),'verified_date':date})
    assert aggregates[-1]['directed_change_patterns']==2496
    assert aggregates[-1]['always_invisible_same_remainder_patterns']==132
    assert aggregates[-1]['additional_zero_alias_patterns']==472
    # Exact serial residue distribution. Counts are integers, not probabilities.
    serial=[1]+[0]*10
    for w in WEIGHTS[4:]:
        nxt=[0]*11
        for r,count in enumerate(serial):
            for digit in range(10):
                nxt[(r+w*digit)%11]+=count
        serial=nxt
    direct_verification='Not requested; use --verify-serials for independent one-million-serial enumeration'
    if args.verify_serials:
        direct=Counter(sum(d*w for d,w in zip(ds,WEIGHTS[4:]))%11 for ds in product(range(10),repeat=6))
        assert [direct[r] for r in range(11)]==serial
        direct_verification='Passed: independently enumerated all 1,000,000 serials and matched all 11 residue counts'
    prefixes=Counter((VALUES[a]+2*VALUES[b]+4*VALUES[c]+8*VALUES['U'])%11 for a,b,c in product(LETTERS,repeat=3))
    totals=[sum(count*serial[(r-pr)%11] for pr,count in prefixes.items()) for r in range(11)]
    assert sum(totals)==26**3*10**6==17_576_000_000
    assert totals[10]==1_597_818_177
    prefix10=Counter()
    for pr,count in prefixes.items():
        prefix10[serial[(10-pr)%11]]+=count
    assert prefix10=={90909:15983,90910:1593}
    residue_rows=[{'raw_remainder':r,'printed_digit':r%10,'hypothetical_U_identifier_bodies':totals[r],'population_bodies':sum(totals),'verified_date':date} for r in range(11)]
    validation=[]
    for number,sid,label in [('BICU1234565','S09','BIC published theoretical example'),('SUDU3070079','S08','GDV published example'),('TEXU4521496','S08','GDV published example'),('BICU0000140','DERIVED','Constructed example; not a claim about a registered or physical unit'),('BICU2000140','DERIVED','Constructed single-serial-digit change; not a claim about a registered or physical unit')]:
        expected=int(number[-1]); actual=check_digit(number[:-1])
        assert actual==expected
        validation.append({'container_number':number,'reference_check_digit':expected,'raw_remainder':raw_remainder(number[:-1]),'recomputed_check_digit':actual,'match':actual==expected,'example_basis':label,'source_id':sid,'source_url':sm[sid]['url'] if sid in sm else '', 'verified_date':date})
    assert validation[-2]['raw_remainder']==0 and validation[-1]['raw_remainder']==10
    worked=[]
    for p,c in enumerate('BICU123456'):
        worked.append({'position':p+1,'character':c,'assigned_value':VALUES[c],'weight':WEIGHTS[p],'product':VALUES[c]*WEIGHTS[p],'verified_date':date})
    assert sum(x['product'] for x in worked)==5494
    # Reproduce every carrier comparison from the preserved metric source values.
    c=inputs['carrier_code_examples']; calc=inputs['calculation_inputs']
    for row in c:
        assert row['max_gross_kg']-row['tare_kg']==row['max_payload_kg']
    kg_per_lb=Decimal(calc['kg_per_avoirdupois_lb'])
    def pounds(kg:int)->Decimal:return Decimal(kg)/kg_per_lb
    def nearest_lb(kg:int)->int:return int(pounds(kg).quantize(Decimal('1'),rounding=ROUND_HALF_UP))
    metrics=[
        ('owner_always_invisible_directed_pairs',44,'44 of 26*25 ordered letter substitutions at a fixed owner position'),
        ('owner_always_invisible_pct',float(Decimal(44)/Decimal(650)*100),'Pattern share, not an empirical error rate or total undetected share'),
        ('unordered_indistinguishable_letter_pairs',22,'Six groups of three and four groups of two letters'),
        ('raw_remainder_ten_bodies',totals[10],'All 26^3 prefixes x 1,000,000 serials, fixed equipment U; no registration filter'),
        ('raw_remainder_ten_pct',float(Decimal(totals[10])/Decimal(sum(totals))*100),'Not unissuable: final rule prints 0'),
        ('volume_increase_pct',float((Decimal(str(c[3]['capacity_cbm']))/Decimal(str(c[0]['capacity_cbm']))-1)*100),'45-foot HC example versus 20-foot standard example'),
        ('volume_ratio',float(Decimal(str(c[3]['capacity_cbm']))/Decimal(str(c[0]['capacity_cbm']))),'45-foot HC volume divided by 20-foot standard volume'),
        ('payload_change_45_vs_20_kg',c[3]['max_payload_kg']-c[0]['max_payload_kg'],'45-foot HC example minus 20-foot standard example'),
        ('payload_range_kg',max(x['max_payload_kg'] for x in c)-min(x['max_payload_kg'] for x in c),'Maximum minus minimum across four examples'),
        ('payload_range_pct_of_minimum',float(Decimal(1050)/Decimal(27700)*100),'Range divided by smallest example payload, not the mean'),
        ('tare_4800_kg_converted_lb',float(pounds(4800)),'Computed using NIST exact kg/lb; source separately prints 10,552 lb'),
        ('tare_4800_kg_rounded_lb',nearest_lb(4800),'Nearest whole pound; not an overwritten source cell'),
        ('iso_668_base_threshold_rounded_lb',nearest_lb(36000),'36,000 kg base-text threshold, nearest whole pound'),
        ('iso_threshold_headroom_under_80000_lb',80000-nearest_lb(36000),'Illustrative subtraction, not a road permit or actual load determination'),
        ('40hc_payload_example_kg',32500-3900,'Metric maximum gross minus metric tare'),
        ('40hc_vs_40standard_payload_change_kg',28600-28750,'Same 32,500 kg gross rating; 150 kg greater tare')]
    summary=[{'metric':n,'value':v,'denominator_or_qualification':q,'verified_date':date} for n,v,q in metrics]
    ratings=[]
    scenarios=[('20-foot carrier example',30480,67197,80000,'S22;S20','Published carrier pounds; comparison assumes an 80,000 lb combination scenario'),('40-foot standard/high-cube examples',32500,71650,80000,'S23;S24;S20','Published carrier pounds; no tractor/chassis mass measured'),('ISO 668:2020 base-text threshold',36000,nearest_lb(36000),80000,'S13;S28;S20','Base standard clause 5.2.2; later amendment not reviewed'),('80,000 lb reference combination',None,80000,80000,'S20','Comparison baseline, not a universal Iowa ceiling'),('Ordinary pneumatic-tire single/tandem axle limits',None,None,None,'S20','20,000 / 34,000 lb; applicable exceptions and all axle/group limits still matter'),('Qualifying noninterstate six-axle combination at 60 ft',None,90000,90000,'S20','22,803 lb difference against 67,197 lb example; statutory conditions apply'),('Qualifying noninterstate seven-axle combination at 62 ft',None,96000,96000,'S20','16,634 lb difference against rounded 79,366 lb threshold; statutory conditions apply')]
    for name,kg,lb,limit,sids,note in scenarios:
        difference=(limit-lb) if name in [x[0] for x in scenarios[:3]] else (22803 if 'six-axle' in name else 16634 if 'seven-axle' in name else '')
        ratings.append({'scenario':name,'marked_or_limit': 'marked/example' if kg else 'legal reference','kg':kg if kg is not None else '', 'lb':lb if lb is not None else '', 'comparison_combination_lb':limit if limit is not None else '', 'remaining_lb_in_stated_comparison':difference,'source_ids':sids,'note':note,'verified_date':date})
    fines=[]
    for row in calc['iowa_fine_rows']:
        over=row['pounds_over']; base=row['schedule_fine']
        if over>20000:assert Decimal(base)==Decimal(2200)+Decimal(over-20000)/10
        fines.append({'pounds_over_legal_gross':over,'axle_schedule_fine_usd':base,'gross_weight_fine_usd':float(Decimal(base)/2),'source_id':'S20','source_url':sm['S20']['url'],'qualification':'Section 321.463(11) schedule only; not a total adjudicated penalty; subsection 11(d) permits other penalties','verified_date':date})
    datasets={'check-digit-audit':letter_rows,'blind-spot-pairs':pair_rows,'single-character-audit':subs,'single-character-summary':aggregates,'transposition-audit':trans,'validation':validation,'worked-check-digit':worked,'raw-remainder-distribution':residue_rows,'check-digit-audit-summary':summary,'rating-vs-iowa-limits':ratings,'iowa-gross-weight-fines':fines,'marking-register':inputs['marking_register'],'marking-size-specs':inputs['marking_size_specs'],'carrier-code-examples':c,'sources':inputs['sources']}
    for name,rows in datasets.items():write_csv(out,name+'.csv',rows)
    bundle={'name':'Shipping Container Markings Reference Set','version':inputs['version'],'verified_date':date,'methodology':{'algorithm_reference':sm['S08']['url'],'licensed_current_iso_text_reviewed':False,'checksum_rule':'raw weighted sum modulo 11; printed digit is raw modulo 10','identifier_space':'26^3 possible owner prefixes, fixed equipment U, six digits 000000–999999; hypothetical, not registered fleet','full_space_count_method':'Exact integer convolution of prefix and serial residue distributions; no sampling','serial_bruteforce_check':direct_verification,'checksum_state_grid':'Each directed change pattern x all eleven starting raw remainders, equally weighted; not an empirical or fleet-weighted error rate','count_of_prefixes_by_remainder10_serial_count':dict(prefix10),'synthetic_examples':'BIC-prefixed constructed counterexamples make no registration, ownership or physical-existence claim','calculation_rounding':'kg/lb factor exact from NIST; converted pounds rounded to nearest whole pound only where identified; carrier-supplied imperial cells retained unchanged'},'tables':datasets}
    (out/'shipping-container-markings-dataset-v1.0.1.json').write_text(json.dumps(bundle,ensure_ascii=False,indent=2)+'\n',encoding='utf-8')
    (out/'reproduction-checks.json').write_text(json.dumps({'status':'passed','serial_independent_check':direct_verification,'hypothetical_bodies':sum(totals),'raw_remainder_ten':totals[10],'unordered_same_remainder_letter_pairs':len(pair_rows),'body_substitution_patterns':2496,'published_examples_validated':3,'synthetic_zero_alias_examples_validated':2,'carrier_metric_payload_arithmetic_checks':4,'derived_fine_rows':len(fines)},indent=2)+'\n')
    print(f'Passed. Wrote {len(datasets)} CSV tables and JSON outputs to {out}')
    print(direct_verification)

if __name__=='__main__':
    main()