"""Reproducible model checks, not physical experiments or a formal proof assistant.
Run: python3 code/verify_models.py
All results are generated mathematical values; no human or experimental data.
"""
from pathlib import Path
from fractions import Fraction as F
import csv
import json
import math

ROOT=Path(__file__).resolve().parents[1]
RESULTS=ROOT/'results'


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


def rk4_constant(y,u,tau,dt):
    def f(v):return (u-v)/tau
    a=f(y);b=f(y+dt*a/2);c=f(y+dt*b/2);d=f(y+dt*c)
    return y+dt*(a+2*b+2*c+d)/6


def check_buffer():
    rows=[]
    for tau in [0.2,1.0,3.0]:
        amplitude=0.7;duration=0.8
        for steps in [100,200]:
            dt=duration/steps;y=0.;error=0.
            for i in range(steps):
                y=rk4_constant(y,amplitude,tau,dt)
                exact=amplitude*(1-math.exp(-(i+1)*dt/tau))
                error=max(error,abs(y-exact))
            peak=amplitude*(1-math.exp(-duration/tau))
            rows.append(dict(data_origin='numerical ODE verification only',tau=tau,input_amplitude=amplitude,pulse_duration=duration,rise_steps=steps,analytic_peak=peak,numerical_end_of_rise=y,max_absolute_rise_error=error,units='consistent arbitrary model units'))
            assert peak<amplitude and error<1e-7
    for tau in [0.2,1.,3.]:
        coarse,fine=[r for r in rows if r['tau']==tau]
        assert fine['max_absolute_rise_error']<coarse['max_absolute_rise_error']
    # Positive-part threshold closure checked by independent time quadrature.
    closure=[]
    amplitude=duration=tau=1.
    dt=0.0005;end=16.
    for threshold in [0.,0.1,0.5,0.7,1.2]:
        total=0.
        # Midpoint quadrature avoids endpoint conventions on the rectangular pulse.
        for i in range(int(end/dt)):
            t=(i+0.5)*dt
            y=1-math.exp(-t) if t<=1 else (1-math.exp(-1))*math.exp(-(t-1))
            total+=max(y-threshold,0.)*dt
        rigid=max(1-threshold,0.)
        assert total<=rigid+1e-6
        closure.append(dict(data_origin='model quadrature; not measured release',closure='positive part',threshold=threshold,rigid_integrated_release=rigid,buffer_integrated_release=total,units='arbitrary model units'))
    threshold=.1
    start=-math.log(1-threshold)
    end=1+math.log((1-math.exp(-1))/threshold)
    binary_duration=end-start
    assert binary_duration>1
    matched_work_peak=math.exp(.5)-math.exp(-.5)
    assert matched_work_peak>1
    write_csv('H0836_ODE_checks.csv',rows)
    write_csv('H0836_convex_release_checks.csv',closure)
    return dict(max_ODE_error=max(r['max_absolute_rise_error'] for r in rows),threshold_counterexample_buffer_duration=binary_duration,threshold_counterexample_rigid_duration=1.0,matched_work_buffer_peak=matched_work_peak,matched_work_rigid_peak=1.0,matched_work_volume_is_also_matched=False)


def check_skin():
    rows=[]
    for b in [F(1,10**6),F(1,1000),F(1,100),F(1),F(100)]:
        for skin in [F(1,1000),F(1,100),F(1,10),F(1)]:
            for ratio in [F(1),F(2),F(10),F(100)]:
                hole=skin*ratio;delta=hole-skin
                for phi in [F(1,1000),F(1,100),F(1,10),F(1,2)]:
                    g=skin+phi*delta
                    t0=1/b+1/skin;t1=1/b+1/g
                    reduction=1-t1/t0
                    algebra=b*phi*delta/(g*(b+skin))
                    assert reduction==algebra
                    condition=delta*(b-phi*(b+skin))>skin*(b+skin)
                    assert (reduction>phi)==condition
                    assert t1<=t0
                    rows.append(dict(data_origin='exact rational model calculation',bulk_conductance=str(b),skin_conductance=str(skin),hole_conductance=str(hole),open_fraction=str(phi),half_time_ratio=str(t1/t0),relative_delay_reduction=str(reduction),reduction_exceeds_open_fraction=condition,units='consistent arbitrary model units'))
    b=F(1);skin=F(1,100);hole=F(1);phi=F(1,100)
    g=skin+phi*(hole-skin)
    D=b*phi*(hole-skin)/(g*(b+skin))
    assert D==F(9900,20099)
    # Explicit bulk-limited counterexample.
    bslow=F(1,10**6)
    Dslow=bslow*phi*(hole-skin)/(g*(bslow+skin))
    assert Dslow<phi
    write_csv('H0346_exact_parameter_checks.csv',rows)
    return dict(exact_cases_checked=len(rows),example_reduction_exact=str(D),example_reduction_float=float(D),example_open_fraction=float(phi),bulk_limited_reduction=float(Dslow),example_units='arbitrary conductance units; not measured material parameters')


def check_layout():
    # Exact coverage computation on [0,1] for sorted rational sites.
    def radius(sites):
        sites=sorted(sites)
        return max([sites[0],1-sites[-1]]+[(b-a)/2 for a,b in zip(sites,sites[1:])])
    distributed=[F(2*i+1,18) for i in range(9)]
    edge=[F(0),F(1,16),F(1,8),F(3,16),F(1,4),F(3,4),F(13,16),F(7,8),F(1)]
    hd=radius(distributed);he=radius(edge)
    assert hd==F(1,18) and he==F(1,4)
    # A pair of fields and opposite bounded noise give identical observations.
    L=F(2);epsilon=F(1,100);M=F(10)
    at=lambda x,S:min(abs(x-s) for s in S)
    adversary=[]
    for name,sites,h in [('distributed',distributed,hd),('edge',edge,he)]:
        for s in sites:
            fplus=M+epsilon+L*at(s,sites)
            fminus=M-epsilon-L*at(s,sites)
            assert fplus-epsilon==fminus+epsilon==M
        error=L*h+epsilon
        adversary.append(dict(data_origin='exact rational minimax certificate',layout=name,reference_count=len(sites),coverage_radius=str(h),L=str(L),noise_bound=str(epsilon),minimax_error=str(error),units='arbitrary distance and field units'))
    # 3x3 square grid theorem: exact squared geometric radii.
    square_grid_h2=F(1,18);boundary_lower_h2=F(1,4)
    assert square_grid_h2<boundary_lower_h2
    write_csv('H0199_exact_coverage_checks.csv',adversary)
    return dict(one_dimensional_distributed_radius=str(hd),one_dimensional_edge_radius=str(he),square_grid_squared_radius=str(square_grid_h2),square_boundary_squared_radius_lower_bound=str(boundary_lower_h2),square_grid_radius=math.sqrt(float(square_grid_h2)),square_boundary_radius_lower_bound=.5,guarantee='worst-case bounded-error point-sampling model only')


def main():
    RESULTS.mkdir(exist_ok=True)
    summary=dict(data_origin='mathematical derivations and model checks only',physical_experiments=0,formal_proof_assistant_used=False,checks={'H0199':check_layout(),'H0836':check_buffer(),'H0346':check_skin()},interpretation='Analytical proofs are in proofs/. Finite computational checks are diagnostics, not substitutes for those proofs or for physical experiments.')
    (RESULTS/'Verification Summary.json').write_text(json.dumps(summary,indent=2))
    print(json.dumps(summary,indent=2))


if __name__=='__main__':main()