"""Conditional model checks, not physical experiments. Standard library only."""
from pathlib import Path
from fractions import Fraction as F
import math,csv,json
ROOT=Path(__file__).resolve().parents[1]
OUT=ROOT/'results'


def solve(a,b):
    rows=[[F(x) for x in row]+[F(v)] for row,v in zip(a,b)]
    n=len(b)
    for col in range(n):
        pivot=next((i for i in range(col,n) if rows[i][col]),None)
        if pivot is None:raise ValueError('Singular: constraints do not identify a unique result')
        rows[col],rows[pivot]=rows[pivot],rows[col]
        divisor=rows[col][col];rows[col]=[x/divisor for x in rows[col]]
        for i in range(n):
            if i!=col:
                f=rows[i][col];rows[i]=[x-f*y for x,y in zip(rows[i],rows[col])]
    return [row[-1] for row in rows]


def thermal_recover(n,edges,anchors,powers):
    L=[[F(0) for _ in range(n)] for _ in range(n)]
    for i,j,k in edges:
        k=F(k);L[i][i]+=k;L[j][j]+=k;L[i][j]-=k;L[j][i]-=k
    unknown=[i for i in range(n) if i not in anchors]
    a=[[L[i][j] for j in unknown] for i in unknown]
    b=[F(powers.get(i,0))-sum(L[i][j]*F(t) for j,t in anchors.items()) for i in unknown]
    vals=solve(a,b)
    result={**{i:F(t) for i,t in anchors.items()},**dict(zip(unknown,vals))}
    residual=[sum(L[i][j]*result[j] for j in range(n))-F(powers.get(i,0)) for i in unknown]
    assert all(x==0 for x in residual)
    return result


def network_checks():
    tests=[('chain_hidden_boundary',4,[(0,1,1),(1,2,1),(2,3,1)],{0:0},{3:4},{0:F(0),1:F(4),2:F(8),3:F(12)}),
           ('nonuniform_chain',4,[(0,1,1),(1,2,2),(2,3,4)],{0:0,3:14},{},{0:F(0),1:F(8),2:F(12),3:F(14)}),
           ('weighted_star',3,[(0,2,1),(1,2,3)],{0:0,1:10},{},{0:F(0),1:F(10),2:F(15,2)}),
           ('equal_star_negative_control',3,[(0,2,1),(1,2,1)],{0:0,1:10},{},{0:F(0),1:F(10),2:F(5)})]
    rows=[]
    for name,n,edges,anchors,powers,expected in tests:
        actual=thermal_recover(n,edges,anchors,powers);assert actual==expected
        for node in sorted(actual):rows.append(dict(data_origin='exact constructed thermal-network model',case=name,node=node,anchor=node in anchors,temperature_exact=str(actual[node]),units='arbitrary consistent thermal units'))
    try:thermal_recover(2,[(0,1,1)],{},{});raise AssertionError('Unanchored graph incorrectly identified')
    except ValueError:pass
    # Anchor-only error cannot exceed the anchor error bound in the harmonic case.
    base=thermal_recover(4,[(0,1,1),(1,2,2),(2,3,4)],{0:0,3:14},{})
    err=F(1,10)
    for e0,e3 in [(-err,-err),(-err,err),(err,-err),(err,err)]:
        noisy=thermal_recover(4,[(0,1,1),(1,2,2),(2,3,4)],{0:e0,3:14+e3},{})
        assert all(abs(noisy[i]-base[i])<=err for i in [1,2])
    # Identifiability failure: one boundary sum permits many pairs.
    assert F(2)+F(8)==F(4)+F(6)==F(10)
    # Unknown power error can be strongly amplified by a weak link.
    weak_link_error=F(1,100)/F(1,1000);assert weak_link_error==10
    write('N0755_exact_network_cases.csv',rows)
    return dict(exact_cases=len(tests),unanchored_case_rejected=True,anchor_error_extremes_checked=4,weak_link_temperature_error_exact=str(weak_link_error),interpretation='Closure and known conductances required; numbers are not physical measurements.')


def step_ode(tf,ts,steps,end=20.):
    dt=end/steps;x=y=0.;maxerr=0.;values=[]
    def advance(v,tau):
        f=lambda z:(1-z)/tau
        a=f(v);b=f(v+dt*a/2);c=f(v+dt*b/2);d=f(v+dt*c)
        return v+dt*(a+2*b+2*c+d)/6
    for i in range(steps):
        x=advance(x,tf);y=advance(y,ts);t=(i+1)*dt
        exact=math.exp(-t/ts)-math.exp(-t/tf)
        maxerr=max(maxerr,abs((x-y)-exact))
        if i%max(1,steps//200)==0:values.append(dict(data_origin='ODE model output; not measured swelling',time=t,fast_path=x,slow_path=y,differential=x-y,analytic_differential=exact,units='arbitrary concentration and displacement units'))
    return maxerr,values


def swelling_checks():
    tf,ts=1.,4.
    tpeak=tf*ts*math.log(ts/tf)/(ts-tf)
    peak=math.exp(-tpeak/ts)-math.exp(-tpeak/tf)
    assert abs(peak-.75*4**(-1/3))<1e-12
    coarse,_=step_ode(tf,ts,1000)
    fine,rows=step_ode(tf,ts,2000)
    assert fine<coarse and fine<1e-8
    freqs=[]
    for i in range(241):
        w=10**(-4+8*i/240)
        z=1j*w
        H=1/(1+z*tf)-1/(1+z*ts)
        analytic=w*(ts-tf)/math.sqrt((1+(w*tf)**2)*(1+(w*ts)**2))
        assert abs(abs(H)-analytic)<1e-12
        freqs.append(dict(data_origin='linear transfer model; not sensor calibration',angular_frequency=w,gain=abs(H),units='arbitrary consistent time and displacement units'))
    optimal=1/math.sqrt(tf*ts)
    H=1/(1+1j*optimal*tf)-1/(1+1j*optimal*ts)
    assert abs(abs(H)-(ts-tf)/(ts+tf))<1e-12
    assert freqs[0]['gain']<.001 and freqs[-1]['gain']<.001
    # Constant-input mismatch defeats DC cancellation.
    mismatch_static=(F(101,100)-F(1))*F(10)
    assert mismatch_static==F(1,10)
    # Equal paths cancel every frequency, independent of input amplitude.
    assert 1/(1+2j)-1/(1+2j)==0
    write('N0751_step_ODE.csv',rows);write('N0751_frequency_response.csv',freqs)
    return dict(step_peak_time=tpeak,step_peak_amplitude=peak,optimal_angular_frequency=optimal,maximum_frequency_gain=abs(H),coarse_RK4_error=coarse,fine_RK4_error=fine,static_gain_mismatch_output_exact=str(mismatch_static),equal_paths_cancel=True,high_frequencies_attenuate=True,interpretation='Calibrated equal gains and valid linear kinetics are assumptions, not measurements.')


def write(name,rows):
    with (OUT/name).open('w',newline='') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)


def main():
    OUT.mkdir(exist_ok=True)
    result=dict(data_origin='mathematical model checks only',physical_experiments=0,formal_proof_assistant_used=False,N0755=network_checks(),N0751=swelling_checks())
    (OUT/'Verification Summary.json').write_text(json.dumps(result,indent=2));print(json.dumps(result,indent=2))


if __name__=='__main__':main()