"""Reproduce Graviton Research Note 01's synthetic kinetic examples.

Python standard library only. These are arbitrary demonstration inputs, not
measurements or evidence of biological plausibility. No consciousness outcome.
Run beside this file to write response-checks.json beside the script.
"""
from decimal import Decimal as D
from pathlib import Path
import hashlib
import json
import math

flux, change, clearance, slope = map(D, ('0.20', '0.010', '0.10', '50'))
dx = flux * change / clearance
dy = slope * dx
assert dx == D('0.020') and dy == D('1.00')
# Independently express concentration in mol/litre rather than micromolar.
dx_molar = D('0.20e-6') * change / clearance
dy_si = D('50e-12') / D('1e-6') * dx_molar
assert dx_molar * D('1e9') == D('20')
assert dy_si * D('1e12') == D('1')
bound_x = D('.25') * D('.001') / D('.08')
bound_y = D('60') * bound_x
assert bound_x * 1000 == D('3.125') and bound_y == D('.1875')
assert bound_y < D('.5')

def analytic(t, shift=.01, rate=.1):
    return .2 * shift / rate * (-math.expm1(-rate*t))

def rk4(t, shift, rate, dt=.01):
    value=0.
    steps=round(t/dt)
    for _ in range(steps):
        def f(x): return .2*shift-rate*x
        a=f(value); b=f(value+dt*a/2); c=f(value+dt*b/2); d=f(value+dt*c)
        value += dt*(a+2*b+2*c+d)/6
    return value

errors=[]
for shift in [-.05, 0, .01, .05]:
    for rate in [.05,.1,.5]:
        for t in [0,1,10,60]:
            errors.append(abs(rk4(t,shift,rate)-analytic(t,shift,rate)))
            assert abs(analytic(t,-shift,rate)+analytic(t,shift,rate))<1e-14
assert max(errors)<1e-11
result={
    'status':'Synthetic mathematical examples; no measured molecular or neural data',
    'assumptions':['Matched formation flux and background production','First-order clearance in a well-mixed compartment','Instantaneous linear calibrated current relation over the relevant interval','Initial concentration difference is zero'],
    'reference_parameters':{'flux_micromolar_per_second':float(flux),'yield_change':float(change),'clearance_per_second':float(clearance),'slope_picoamp_per_micromolar':float(slope)},
    'steady_result':{'concentration_micromolar':float(dx),'concentration_nanomolar':float(dx*1000),'current_picoamp':float(dy)},
    'separate_exclusion_example':{'maximum_flux_micromolar_per_second':.25,'maximum_abs_yield_change':.001,'minimum_clearance_per_second':.08,'maximum_abs_slope_picoamp_per_micromolar':60,'concentration_bound_nanomolar':float(bound_x*1000),'current_bound_picoamp':float(bound_y),'illustrative_meaningful_threshold_picoamp':.5},
    'validation':{'decimal_and_SI_units_passed':True,'independent_RK4_max_error_micromolar':max(errors),'tolerance_micromolar':1e-11,'opposite_shift_symmetry_passed':True},
    'transient_reference':[{'seconds':t/2,'concentration_nanomolar':1000*analytic(t/2),'current_picoamp':50*analytic(t/2)} for t in range(121)],
    'generator_sha256':hashlib.sha256(Path(__file__).read_bytes()).hexdigest()
}
target=Path(__file__).with_name('response-checks.json')
target.write_text(json.dumps(result,indent=2)+'\n')
print(json.dumps({'steady_result':result['steady_result'],'exclusion':result['separate_exclusion_example'],'validation':result['validation']},indent=2))