"""Reproducible NUMERICAL SIMULATION, not camera or participant measurements.
Run: python run_simulation.py   (requires numpy and matplotlib)
All adjustable assumptions are in parameters.json. No data are fitted to real observations.
"""
from pathlib import Path
import json, csv, platform, hashlib
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

ROOT = Path(__file__).resolve().parent
OUT = ROOT / 'results'
OUT.mkdir(exist_ok=True)
P = json.loads((ROOT / 'parameters.json').read_text(encoding='utf-8'))
rng = np.random.Generator(np.random.PCG64(P['seed']))
checks = []
def check(name, condition, value=None):
    checks.append({'name': name, 'passed': bool(condition), 'value': value})
    if not condition: raise AssertionError(name)
def save_csv(name, rows):
    with (OUT / name).open('w', encoding='utf-8-sig', newline='') as f:
        w = csv.DictWriter(f, fieldnames=list(rows[0])); w.writeheader(); w.writerows(rows)
def save_json(name, value):
    (OUT / name).write_text(json.dumps(value, ensure_ascii=False, indent=2, allow_nan=False), encoding='utf-8')

# A. Project 3D vertices, then sample their silhouettes at pixel centres.
# The analytic width ratio is NOT used by the projection/raster measurement.
A = P['projection']
K = np.array([[A['f_px'],0,A['width_px']/2],[0,A['f_px'],A['height_px']/2],[0,0,1.]])
def corners(cx, z):
    w,h = A['card_width_cm']/2,A['card_height_cm']/2
    return np.array([[cx-w,-h,z],[cx+w,-h,z],[cx+w,h,z],[cx-w,h,z]])
def project(vertices, phase):
    uvw = vertices @ K.T
    uv = uvw[:,:2] / uvw[:,2:]
    uv[:,0] += phase
    return uv
def raster_width(uv):
    # Axis-aligned rectangles; the central horizontal row crosses both silhouettes.
    centres = np.arange(A['width_px']) + .5
    occupied = (centres >= uv[:,0].min()) & (centres < uv[:,0].max())
    indices = np.flatnonzero(occupied)
    if not len(indices): raise ValueError('Object has no occupied pixel')
    return int(indices[-1]-indices[0]+1)
rows=[]; continuous_errors=[]
for depth in [0., A['depth_cm']]:
    for z in A['z_cm']:
        for swapped in [False,True]:
            for phase in np.arange(A['phases']) / A['phases']:
                sign = -1 if swapped else 1
                front = project(corners(-sign*A['centre_offset_cm'],z),phase)
                back = project(corners(sign*A['centre_offset_cm'],z+depth),phase)
                for uv in [front,back]:
                    assert uv[:,0].min()>0 and uv[:,0].max()<A['width_px']
                    assert uv[:,1].min()>0 and uv[:,1].max()<A['height_px']
                continuous = np.ptp(front[:,0])/np.ptp(back[:,0])
                theory = 1+depth/z
                continuous_errors.append(abs(continuous-theory))
                wf,wb = raster_width(front),raster_width(back)
                rows.append(dict(data_source='NUMERICAL_SIMULATION',z_cm=z,depth_cm=depth,
                    swapped=int(swapped),phase_px=float(phase),front_width_px=wf,back_width_px=wb,
                    continuous_ratio=float(continuous),analytic_ratio=theory,raster_ratio=wf/wb,
                    raster_error_percentage_points=100*(wf/wb-theory)))
save_csv('projection_all.csv',rows)
summary_A=[]
for z in A['z_cm']:
    r=[x for x in rows if x['z_cm']==z and x['depth_cm']==A['depth_cm']]
    vals=np.array([x['raster_ratio'] for x in r])
    summary_A.append(dict(data_source='NUMERICAL_SIMULATION',z_cm=z,
        analytic_ratio=1+A['depth_cm']/z,virtual_configurations=len(r),
        raster_min=float(vals.min()),raster_max=float(vals.max()),
        max_abs_error_percentage_points=float(max(abs(x['raster_error_percentage_points']) for x in r))))
save_csv('projection_summary.csv',summary_A)
check('3D vertex projection agrees with analytic ratio', max(continuous_errors)<1e-12, max(continuous_errors))
control=max(abs(x['raster_ratio']-1) for x in rows if x['depth_cm']==0)
check('Equal-depth continuous ratio is one',all(abs(x['continuous_ratio']-1)<1e-12 for x in rows if x['depth_cm']==0))
# Integer widths need not match at different horizontal subpixel alignments.
check('All raster dimensions positive',all(x['front_width_px']>0 and x['back_width_px']>0 for x in rows))

# B. Propagate hypothetical calibration and pixel uncertainty with Monte Carlo.
# Five-standard-deviation truncation ensures bounded, positive input diameters.
B=P['uncertainty']; n=B['draws']; truncation_redraws=0
def noise(scale):
    global truncation_redraws
    x=rng.standard_normal(n)
    bad=np.abs(x)>B['normal_cutoff_sd']
    while bad.any():
        truncation_redraws+=int(bad.sum()); x[bad]=rng.standard_normal(int(bad.sum()))
        bad=np.abs(x)>B['normal_cutoff_sd']
    return scale*x
z0=B['Z0_cm']; C=B['C_px_cm']; p0=C/z0
cal_z=z0+noise(B['u_Z0_cm']); cal_p=p0+noise(B['u_p0_px']); cal_C=cal_z*cal_p
c=np.hypot(B['u_Z0_cm']/z0,B['u_p0_px']/p0)
raw={'Z0_draw_cm':cal_z,'p0_draw_px':cal_p,'C_draw_px_cm':cal_C,
     'Z_truth_cm':np.array(B['Z_cm']),'sigma_p_px':np.array(B['sigma_p_px'])}
summary_B=[]; p_all=[]; z_all=[]
for sigma in B['sigma_p_px']:
    ps=[]; zs=[]
    for Z in B['Z_cm']:
        p=C/Z+noise(sigma); estimate=cal_C/p
        check(f'Positive inputs sigma={sigma}, Z={Z}',bool((p>0).all() and (cal_C>0).all()))
        ps.append(p); zs.append(estimate)
        rel=(estimate/Z-1)*100
        analytic=np.hypot(c,sigma*Z/C)*100
        observed=float(np.std(estimate,ddof=1)/Z*100)
        summary_B.append(dict(data_source='NUMERICAL_SIMULATION',sigma_p_px=sigma,Z_cm=Z,
            computational_draws=n,mean_estimated_cm=float(estimate.mean()),
            bias_pct=float(rel.mean()),relative_sd_pct=observed,
            first_order_relative_u_pct=float(analytic),mape_pct=float(np.mean(np.abs(rel))),
            error_q025_pct=float(np.quantile(rel,.025)),error_q975_pct=float(np.quantile(rel,.975))))
    p_all.append(ps);z_all.append(zs)
raw['p_draw_px']=np.array(p_all);raw['Z_estimated_cm']=np.array(z_all)
np.savez_compressed(OUT/'monte_carlo_draws.npz',**raw)
save_csv('uncertainty_summary.csv',summary_B)
check('Monte Carlo output dimensions',raw['Z_estimated_cm'].shape==(len(B['sigma_p_px']),len(B['Z_cm']),n))
# Independently obtain the derivatives numerically, not from the propagation formula.
Ztest=80.; x=np.array([z0,p0,C/Ztest]); u=np.array([B['u_Z0_cm'],B['u_p0_px'],.5])
fn=lambda v:v[0]*v[1]/v[2]
gradient=[]
for j in range(3):
    h=1e-5*abs(x[j]); xp=x.copy(); xm=x.copy(); xp[j]+=h;xm[j]-=h
    gradient.append((fn(xp)-fn(xm))/(2*h))
fd=np.linalg.norm(np.array(gradient)*u)/Ztest
formula=np.hypot(c,.5*Ztest/C)
check('Finite-difference uncertainty agrees with analytic derivative',abs(fd-formula)<1e-9,float(abs(fd-formula)))

# B2. Hypothetical effective scale drift. This is not an autofocus measurement.
drift=[]
for beta in [-.03,0,.03]:
    for Z in B['Z_cm']:
        factor=1+beta*(Z-z0)/(max(B['Z_cm'])-z0)
        estimate=Z/factor
        drift.append(dict(data_source='NUMERICAL_SIMULATION',beta_at_120_cm=beta,
            Z_cm=Z,scale_factor=factor,estimated_cm=estimate,bias_pct=100*(estimate/Z-1)))
save_csv('scale_drift.csv',drift)

# C. Conditional feasible intervals, checked against direct inequalities on a grid.
D=P['range']; summary_C=[]
grid=np.arange(.01,600.001,.01)
for case in D['cases']:
    Ccase=case['C_px_cm']; delta=case['depth_cm']; tolerance=case['tolerance']
    cc=np.hypot(B['u_Z0_cm']/z0,B['u_p0_px']/(Ccase/z0))
    lower=D['d_cm']+delta/D['epsilon']; pix=D['k']*Ccase/D['Pmin_px']
    blocked=cc>=tolerance
    unc=0. if blocked else Ccase/D['sigma_p_px']*np.sqrt(tolerance**2-cc**2)
    upper=min(pix,unc); empty=bool(blocked or lower>upper)
    valid=(grid>D['d_cm']) & (delta/(grid-D['d_cm'])<=D['epsilon']+1e-12) & (D['k']*Ccase/grid>=D['Pmin_px']-1e-12) & (np.hypot(cc,D['sigma_p_px']*grid/Ccase)<=tolerance+1e-12)
    points=grid[valid]
    check('Grid emptiness agrees: '+case['case'],empty==(len(points)==0))
    if len(points):
        check('Grid boundaries agree within 0.01 cm: '+case['case'],abs(points[0]-lower)<=.010001 and abs(points[-1]-upper)<=.010001)
    summary_C.append(dict(data_source='NUMERICAL_SIMULATION',case=case['case'],C_px_cm=Ccase,
        depth_cm=delta,tolerance_pct=100*tolerance,calibration_u_pct=float(cc*100),
        lower_cm=lower,pixel_upper_cm=pix,uncertainty_upper_cm=float(unc),
        calibration_blocked=bool(blocked),empty=empty,upper_cm=float(upper)))
save_csv('interval_summary.csv',summary_C)

# Figures use calculated results; every figure identifies its provenance.
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False,
    'axes.grid':True,'grid.alpha':.16,'figure.facecolor':'#fcfcf9','axes.facecolor':'#fcfcf9','savefig.facecolor':'#fcfcf9'})
def save(fig,name):
    fig.savefig(OUT/(name+'.png'),dpi=180,bbox_inches='tight'); plt.close(fig)
fig,ax=plt.subplots(figsize=(7.4,3.65))
line=np.linspace(25,120,400)
ax.plot(line,100*A['depth_cm']/line,color='#17665b',label='Analytic pinhole model')
x=np.array(A['z_cm']); lo=np.array([r['raster_min'] for r in summary_A]); hi=np.array([r['raster_max'] for r in summary_A]); mid=(lo+hi)/2
ax.errorbar(x,100*(mid-1),yerr=50*(hi-lo),fmt='o',capsize=4,color='#b56431',label='Raster min–max (32 configurations / distance)')
ax.set(xlabel='Front plane distance z (cm)',ylabel='Ratio change from orthographic (%)',title='A  |  Numerical simulation: two parallel cards')
ax.legend(fontsize=8);save(fig,'projection')
fig,ax=plt.subplots(figsize=(7.4,3.65))
for sigma,color in zip(B['sigma_p_px'],['#819d37','#17665b','#b56431']):
    r=[r for r in summary_B if r['sigma_p_px']==sigma]
    ax.plot([v['Z_cm'] for v in r],[v['first_order_relative_u_pct'] for v in r],color=color,label=f'Assumed sigma = {sigma:.1f} px; first order')
    ax.scatter([v['Z_cm'] for v in r],[v['relative_sd_pct'] for v in r],s=25,color=color)
ax.plot([],[],'.',color='#222',label='Dots: 20,000 computational draws / condition')
ax.set(xlabel='True simulated distance Z (cm)',ylabel='Relative standard uncertainty / SD (%)',title='B  |  Numerical simulation: distance uncertainty')
ax.legend(fontsize=8);save(fig,'uncertainty')
fig,ax=plt.subplots(figsize=(7.4,3.5))
for i,r in enumerate(summary_C):
    label=f"{r['case']}  C={r['C_px_cm']}, depth={r['depth_cm']}, limit={r['tolerance_pct']:g}%"
    ax.plot([r['lower_cm'],r['lower_cm']],[i-.18,i+.18],color='#b56431',lw=2)
    ax.plot([r['upper_cm'],r['upper_cm']],[i-.18,i+.18],color='#17665b',lw=2)
    if not r['empty']:ax.plot([r['lower_cm'],r['upper_cm']],[i,i],color='#17665b',lw=6,solid_capstyle='round')
    else:ax.text(max(r['lower_cm'],r['upper_cm'])+4,i,'No overlap',va='center',fontsize=8,color='#a34a30')
ax.set_yticks(range(len(summary_C)),[f"{r['case']}  |  C={r['C_px_cm']}; depth={r['depth_cm']}; limit={r['tolerance_pct']:g}%" for r in summary_C])
ax.set(xlabel='Conditional distance interval Z (cm)',xlim=(0,220),title='C  |  Numerical simulation: feasible intervals')
ax.invert_yaxis();save(fig,'intervals')

# A virtual raster scene, not a camera photograph.
z=60.; xx=np.arange(A['width_px'])+.5; yy=np.arange(A['height_px'])+.5
canvas=np.ones((A['height_px'],A['width_px'],3))*.96
for cx,depth,color in [(-A['centre_offset_cm'],z,(.09,.4,.35)),(A['centre_offset_cm'],z+A['depth_cm'],(.7,.39,.18))]:
    uv=project(corners(cx,depth),.25)
    mask=(yy[:,None]>=uv[:,1].min())&(yy[:,None]<uv[:,1].max())&(xx[None,:]>=uv[:,0].min())&(xx[None,:]<uv[:,0].max())
    canvas[mask]=color
fig,ax=plt.subplots(figsize=(7.4,3.3));ax.imshow(canvas);ax.grid(False);ax.set_axis_off();ax.set_title('VIRTUAL IMAGE | z = 60 cm; depth separation = 10 cm')
save(fig,'virtual_scene')
summary={'data_source':'NUMERICAL_SIMULATION','seed':P['seed'],'projection':summary_A,'uncertainty':summary_B,
         'intervals':summary_C,'scale_drift':drift,'equal_depth_raster_max_error_pp':100*control,
         'projection_max_continuous_error':max(continuous_errors),'truncation_redraws':truncation_redraws,
         'calibration_relative_u':float(c),'checks_passed':len(checks),
         'runtime':{'python':platform.python_version(),'numpy':np.__version__,'matplotlib':matplotlib.__version__}}
save_json('summary.json',summary); save_json('checks.json',checks)
hashes={str(p.relative_to(ROOT)):hashlib.sha256(p.read_bytes()).hexdigest() for p in [ROOT/'run_simulation.py',ROOT/'parameters.json']+sorted(OUT.glob('*.csv'))+[OUT/'monte_carlo_draws.npz',OUT/'summary.json',OUT/'checks.json']}
save_json('sha256.json',hashes)
print(json.dumps({'checks_passed':len(checks),'projection':summary_A,'sigma_05':[r for r in summary_B if r['sigma_p_px']==.5],'intervals':summary_C,'runtime':summary['runtime']},indent=2))