"""Fictional parcel experiment. Transparent teaching models, not production ML.
Python 3.12 + numpy 2.2.6. No training on validation/test; no internet or database.
Run python ml_reference.py --source . --output experiment.json [--final-test]
Without --final-test, held-out test targets are not scored. Fresh output required.
"""
import argparse,csv,json,hashlib
from pathlib import Path
import numpy as np

FEATURES=['distance_km','load_ratio','zone']
FORBIDDEN=['actual_hours','late','delivered_at','post_delivery_rating']

def read(path):
    with Path(path).open(encoding='utf-8',newline='') as f:return list(csv.DictReader(f))

def fit_preprocessor(rows):
    medians=[];means=[];scales=[]
    for k in FEATURES[:2]:
        vals=np.array([float(r[k]) for r in rows if r[k]!=''])
        if not len(vals): raise ValueError('Entirely missing training feature: '+k)
        if not np.isfinite(vals).all(): raise ValueError('Non-finite training feature: '+k)
        med=float(np.median(vals));a=np.array([float(r[k]) if r[k]!='' else med for r in rows])
        medians.append(med);means.append(float(a.mean()));scales.append(float(a.std()) or 1.)
    return {'medians':medians,'means':means,'scales':scales,'zones':sorted({r['zone'] for r in rows})}

def transform(rows,p):
    # Unknown zones use all-zero indicators; report their count separately.
    return np.array([[( (float(r[k]) if r[k]!='' else p['medians'][i])-p['means'][i])/p['scales'][i] for i,k in enumerate(FEATURES[:2])]+[float(r['zone']==z) for z in p['zones'][1:]] for r in rows])

def fit_stump(x,y):
    # Depth-one classification tree. Enumerate midpoints, minimise weighted Gini.
    candidates=[]
    for j in range(x.shape[1]):
        vals=np.unique(x[:,j])
        for t in (vals[:-1]+vals[1:])/2:
            mask=x[:,j]<=t
            left=float(y[mask].mean());right=float(y[~mask].mean())
            loss=float((mask.sum()*2*left*(1-left)+(~mask).sum()*2*right*(1-right))/len(y))
            candidates.append((loss,j,float(t),left,right))
    if not candidates:return (0.,0,float('inf'),float(y.mean()),float(y.mean()))
    return min(candidates)

def probabilities(x,tree):return np.where(x[:,tree[1]]<=tree[2],tree[3],tree[4])

def metrics(y,p):
    y=np.array(y,dtype=int);p=np.array(p,dtype=int)
    tp=int(((y==1)&(p==1)).sum());fp=int(((y==0)&(p==1)).sum());fn=int(((y==1)&(p==0)).sum());tn=int(((y==0)&(p==0)).sum())
    return {'n':len(y),'positives':int(y.sum()),'tp':tp,'fp':fp,'fn':fn,'tn':tn,'accuracy':float((y==p).mean()),'precision':tp/(tp+fp) if tp+fp else None,'recall':tp/(tp+fn) if tp+fn else None}

def train(rows):
    p=fit_preprocessor(rows);x=transform(rows,p);h=np.array([float(r['actual_hours']) for r in rows]);y=np.array([int(r['late']) for r in rows]);coef=np.linalg.lstsq(np.c_[np.ones(len(x)),x],h,rcond=None)[0]
    return p,coef,fit_stump(x,y),float(h.mean()),int(y.mean()>=.5)

def evaluate(rows,fit,threshold=.5):
    p,coef,tree,mean,majority=fit;x=transform(rows,p);h=np.array([float(r['actual_hours']) for r in rows]);y=np.array([int(r['late']) for r in rows]);pred=np.c_[np.ones(len(x)),x]@coef
    return {'regression_mae_hours':float(np.abs(h-pred).mean()),'baseline_mae_hours':float(np.abs(h-mean).mean()),'classification':metrics(y,probabilities(x,tree)>=threshold),'baseline_classification':metrics(y,np.full(len(y),majority)),'unknown_zone_rows':sum(r['zone'] not in p['zones'] for r in rows)}


if __name__ == '__main__':
    pilot=read('pilot_parcels.csv')
    print(json.dumps(evaluate(pilot[40:],train(pilot[:40])),indent=2))