"""Author-run, scoped source audit. No fitting or independent validation.
Run from any directory: python reproduce.py
Original implementation: Christopher Hill / NetPhysics, AI-assisted, MIT.
"""
from pathlib import Path
import csv, json, math, hashlib, platform
import numpy as np
import scipy
from scipy.integrate import solve_ivp
import openpyxl
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

ROOT=Path(__file__).resolve().parent
OUT=ROOT/'output'; OUT.mkdir(exist_ok=True)

def dump(name,obj):
    (OUT/name).write_text(json.dumps(obj,indent=2,allow_nan=False)+'\n')
def write_csv(name,rows):
    with (OUT/name).open('w',newline='') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
def xe(u,vb,vm):return 1/(1+math.exp(-(u-vb)/vm))
def table_tau(u,vb,vm):
    x=xe(u,vb,vm)
    return (64.1*(x*(1-x)/vm)**1.85+0.10)/1000

def end_state(n,ton,toff,tau_on,tau_off,x0,x1):
    # Eq. 15 has a_0 for pulse 1. Here n counts physical pulses starting at 1.
    a=math.exp(-ton/tau_on);b=math.exp(-toff/tau_off)
    af=a*(1-b)/(-math.expm1(-ton/tau_on-toff/tau_off))
    an=af+(a-af)*(a*b)**(n-1)
    return x1+(x0-x1)*an

def fixed(ton,toff,tau_on,tau_off,x0,x1):
    a=math.exp(-ton/tau_on);b=math.exp(-toff/tau_off)
    return x1+(x0-x1)*a*(1-b)/(-math.expm1(-ton/tau_on-toff/tau_off))

# S1: gL,gH in microSiemens; VB,Vm in volts. Shared K=64.1 ms V^gamma,
# gamma=1.85 and t0=0.10 ms. Four voltage-specific calibrations, not one fit.
PARAMS={2:(.50,9.50,1.02,.51),3:(.43,11.3,1.23,.68),4:(.43,12.4,1.37,.81),5:(.44,13.2,1.44,.89)}
wb=openpyxl.load_workbook(ROOT/'input/Fig_4c.xlsx',data_only=True)
observations=[];pulses=[];results=[];checks=[]
for sheet in wb:
    amp=int(sheet.title[:-1]);rows=[]
    for rownum,row in enumerate(sheet.values,1):
        if not all(isinstance(v,(int,float)) for v in row[:3]):continue
        t,v,i=row[:3];rows.append([t,v,i*1000]) # source mA -> microA
        observations.append(dict(trace=sheet.title,source_row=rownum,time_s=t,voltage_V=v,current_uA=i*1000))
    a=np.array(rows);assert np.all(np.diff(a[:,0])>0)
    high=a[:,1]>.9*amp
    starts=np.flatnonzero(high&~np.r_[False,high[:-1]])
    ends=np.flatnonzero(high&~np.r_[high[1:],False])
    assert len(starts)==len(ends)==11
    for n,(start,end) in enumerate(zip(starts,ends),1):
        pulses.append(dict(trace=sheet.title,pulse=n,first_high_sample_s=a[start,0],last_high_sample_s=a[end,0],last_high_current_uA=a[end,2]))
    if amp==1:continue # S1 does not supply 1-V parameters; do not invent them.
    gl,gh,vb,vm=PARAMS[amp];x0=xe(0,vb,vm);x1=xe(amp,vb,vm)
    t0=table_tau(0,vb,vm);t1=table_tau(amp,vb,vm)
    xf=fixed(.002,.002,t1,t0,x0,x1);x11=end_state(11,.002,.002,t1,t0,x0,x1)
    ipred=(gl+(gh-gl)*x11)*amp
    results.append(dict(voltage_V=amp,gL_uS=gl,gH_uS=gh,VB_V=vb,Vm_V=vm,tau_off_s=t0,tau_on_s=t1,x_equilibrium_off=x0,x_equilibrium_on=x1,x_fixed=xf,state_gain=xf/x0,current_gain=(gl+(gh-gl)*xf)/(gl+(gh-gl)*x0),predicted_pulse11_end_uA=ipred,predicted_asymptote_uA=(gl+(gh-gl)*xf)*amp,observed_pulse11_last_high_uA=a[ends[-1],2],observed_last5_high_median_uA=float(np.median(a[ends[-5:],2])),relative_difference_pulse11_percent=100*(ipred/a[ends[-1],2]-1)))
    # Independent numerical integration of on/off ODE segments.
    state=x0;err=0
    for n in range(1,12):
        state=solve_ivp(lambda t,y:(x1-y)/t1,(0,.002),[state],rtol=1e-11,atol=1e-13).y[0,-1]
        err=max(err,abs(state-end_state(n,.002,.002,t1,t0,x0,x1)))
        state=solve_ivp(lambda t,y:(x0-y)/t0,(0,.002),[state],rtol=1e-11,atol=1e-13).y[0,-1]
    checks.append(dict(check='S1 recurrence vs ODE',voltage_V=amp,max_absolute_state_error=err,pass_check=bool(err<1e-9)))
write_csv('observations.csv',observations);write_csv('pulse-endpoints.csv',pulses);write_csv('source-parameter-comparison.csv',results)
# Limiting cases: zero rest approaches x1; zero pulse holds x0.
assert abs(fixed(.002,0,.001,.001,.1,.9)-.9)<1e-15
assert abs(fixed(0,.002,.001,.001,.1,.9)-.1)<1e-15
assert abs(end_state(1,.002,.002,.003,.004,.1,.9)-(.9+(.1-.9)*math.exp(-2/3)))<1e-15

# Figure 3 caption literal values. Alternative is explicitly Figure 1's function,
# never a silently corrected reconstruction of Figure 3.
def fig1_tau(u):
    return .001/(1+math.exp((u-1)/.2)+math.exp(-(u+1)/.2))
unit_cases=[]
for ton,toff,label in [(.1,.05,'3a/3b'),(.05,.2,'3c/3d')]:
    for mode in ['literal Figure 3 caption','Figure 1 functions, Figure 3 voltages and durations']:
        if mode.startswith('literal'):t0,t1,x0,x1=.00097,.98,.05,1.
        else:t0,t1,x0,x1=fig1_tau(.2),fig1_tau(1.2),xe(.2,.5,.1),xe(1.2,.5,.1)
        unit_cases.append(dict(panel=label,interpretation=mode,pulse_s=ton,rest_s=toff,tau_off_s=t0,tau_on_s=t1,x_fixed=fixed(ton,toff,t1,t0,x0,x1)))
write_csv('figure3-unit-audit.csv',unit_cases)

# SI extension of M03-D at one geometry. Synthetic values ONLY; no specimen fit.
M=.02;K=8.;C=.096;Q=.005;beta=1.5;tstar=math.sqrt(M/K)
alphaE=.1/(K*Q**2);tau=15*tstar;cbar=C/math.sqrt(M*K);abar=alphaE*K*Q**2
# Use a smooth drive here to isolate unit mapping from switching errors.
def force(t):return K*Q*math.sin(t/tstar)
def si(t,y):
    q,v,m=y
    return [v,(force(t)-C*v-K*(1+beta*m)*q)/M,alphaE*C*v*v*(1-m)-m/tau]
def normalized(s,y):
    q,p,m=y
    return [p,math.sin(s)-cbar*p-(1+beta*m)*q,abar*cbar*p*p*(1-m)-m/15]
grid=np.linspace(0,30,601)
u=solve_ivp(si,(0,30*tstar),[0,0,0],t_eval=grid*tstar,rtol=2e-11,atol=1e-13)
v=solve_ivp(normalized,(0,30),[0,0,0],t_eval=grid,rtol=2e-11,atol=1e-13)
assert u.success and v.success
scaled=u.y/np.array([Q,Q/tstar,1])[:,None]
errors=np.max(np.abs(scaled-v.y),axis=1)
assert np.max(errors)<1e-8 and np.min(u.y[2])>=0 and np.max(u.y[2])<=1
# Energy derivative identity from the vector field, not a thermodynamic closure.
energy_errors=[]
for t,y in zip(u.t,u.y.T):
    q,vel,m=y;dq,dv,dm=si(t,y)
    derivative=M*vel*dv+K*(1+beta*m)*q*dq+.5*K*beta*q*q*dm
    ledger=force(t)*vel-C*vel*vel+.5*K*beta*q*q*dm
    energy_errors.append(abs(derivative-ledger))
units=dict(status='Synthetic dimensional consistency check; physical calibration pending',M_kg=M,K0_N_per_m=K,C_N_s_per_m=C,Q_m=Q,tstar_s=tstar,alpha_per_J=alphaE,tau_s=tau,c_normalized=cbar,alpha_normalized=abar,max_normalized_errors=dict(zip(['q','p','m'],map(float,errors))),max_energy_identity_residual_W=max(energy_errors))
dump('si-check.json',units)
write_csv('si-trajectories.csv',[dict(time_s=float(t),displacement_m=float(y[0]),velocity_m_per_s=float(y[1]),state=float(y[2])) for t,y in zip(u.t,u.y.T)])
# Cycle-average illustration: saturation + finite retention alter a frequency ratio.
def averaged(rate,T,retention):
    total=rate+1/retention
    return rate/total*(-math.expm1(-total*T))
ratios=[]
for write_rate,retention in [(1e-7,1e9),(1e-7,300),(1e-4,3000),(1e-3,3000)]:
    low=averaged(write_rate,3000,retention);highN=averaged(100*write_rate,300,retention);highT=averaged(100*write_rate,3000,retention)
    ratios.append(dict(low_frequency_write_rate_per_s=write_rate,retention_s=retention,matched_cycles_high_over_low=highN/low,matched_duration_high_over_low=highT/low))
write_csv('frequency-ratio-sensitivity.csv',ratios)
# Jacobian sanity check: Gaussian log10-y density on common positive-y support.
from scipy.integrate import quad
sigma=.3;mu=.7
mass=quad(lambda y:math.exp(-.5*((math.log10(y)-mu)/sigma)**2)/(math.sqrt(2*math.pi)*sigma*y*math.log(10)),0,np.inf)[0]
assert abs(mass-1)<1e-8
checks.extend([dict(check='SI/normalized equivalence',pass_check=True,max_absolute_normalized_error=float(max(errors))),dict(check='log10 response Jacobian integrates to one',pass_check=True,integral=mass),dict(check='zero-rest, zero-pulse and first-pulse limits',pass_check=True)])
assert all(c['pass_check'] for c in checks)
dump('checks.json',checks);dump('source-parameter-comparison.json',results)
dump('environment.json',dict(python=platform.python_version(),numpy=np.__version__,scipy=scipy.__version__,openpyxl=openpyxl.__version__,matplotlib=matplotlib.__version__))
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':11,'axes.spines.top':False,'axes.spines.right':False,'svg.hashsalt':'netphysics060'})
fig,ax=plt.subplots(figsize=(10,6.2),layout='constrained')
for amp in range(2,6):
    rr=[r for r in pulses if r['trace']==f'{amp}V'];n=np.arange(1,12);gl,gh,vb,vm=PARAMS[amp]
    xp=[end_state(int(i),.002,.002,table_tau(amp,vb,vm),table_tau(0,vb,vm),xe(0,vb,vm),xe(amp,vb,vm)) for i in n]
    line=ax.plot(n,[(gl+(gh-gl)*x)*amp for x in xp],label=f'{amp} V · source parameters')[0]
    ax.scatter(n,[r['last_high_current_uA'] for r in rr],color=line.get_color(),marker='x',s=40)
ax.set(xlabel='Pulse number',ylabel='Current (µA)',title='Published parameters and deposited measurements\nLines: ideal 2 ms on / 2 ms off, u ≈ V · Crosses: last high sample')

ax.legend(loc='upper center',bbox_to_anchor=(.5,-.16),ncol=4,fontsize=9,frameon=False);ax.grid(alpha=.15)
fig.savefig(OUT/'pulse-comparison.svg',metadata={'Date':None});fig.savefig(OUT/'pulse-comparison.png',dpi=150);plt.close(fig)
print(json.dumps({'observations':len(observations),'source_parameter_voltages':len(results),'checks_passed':len(checks),'maximum_recurrence_error':max(c.get('max_absolute_state_error',0) for c in checks),'comparison':results,'si':units},indent=2))