"""Preserving native function while measuring a shared cofactor pool."""
from fractions import Fraction as F
from pathlib import Path
import argparse,csv,hashlib,json,platform
from dataclasses import asdict
import numpy as np
from scipy.optimize import brentq
from reactor import LinearRate,MichaelisMenten,Reporter,AssociationSchedule,CofactorReactor
from certificates import DesignBox,RecoveryCertificate,slope_sandwich

# EDITABLE INPUTS: dimensionless synthetic design requirements, not fitted rates.
REGENERATION=1.
NATIVE_USE=1.
COFACTOR_POOL=1.
REPORTER_TOTAL=1.
REPORTER_ACTIVITY=.08
CATALYTIC_RELEASE=20.
DISSOCIATION=1.
PREPARATION_OFFSET=0.
ACQUISITION='8'
RECOVERY='5'
FADING_RATE='21'
RECOVERY_TARGET='0.000003'
PRESERVATION_ALLOWANCE='0.05'
OBSERVED_RECORDS=('0.30549','0.41151','0.361526225','-0.1','0.6')
FUNCTIONAL_THRESHOLD='0.6'
CONCENTRATION_SCALE_MICROMOLAR='1'
TIME_SCALE_SECONDS='60'
VOLUME_MICROLITRES='100'
INTERIOR_DRAWS_PER_CLASS=8
SEED=17092026
# Primitive association is derived from u, beta, c and reporter total; do not
# vary it independently while asserting that effective activity remains fixed.
UNCERTAINTY=DesignBox()
MANUSCRIPT_SHA256='bd27ad3fee701238aedb2442336de04b53187feb6e7b9cd23072c85c3491adbf'


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'))
    out=parser.parse_args().output;out.mkdir(parents=True,exist_ok=True)
    def dump(name,data):
        (out/name).write_text(json.dumps(data,indent=2,default=lambda x:{'exact':str(x),'decimal':float(x)} if isinstance(x,F) else str(x))+'\n',encoding='utf-8')
    def table(name,header,rows):
        with (out/name).open('w',newline='',encoding='utf-8') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
    T=float(F(ACQUISITION));tau=float(F(RECOVERY))
    def nominal(c=CATALYTIC_RELEASE,p=REGENERATION,k=NATIVE_USE,background=0):
        return CofactorReactor(LinearRate(p),LinearRate(k),Reporter(REPORTER_TOTAL,c,DISSOCIATION,REPORTER_ACTIVITY),COFACTOR_POOL,background)
    comparisons=[];traces={}
    for name,c,delay in [('slow',1.,0.),('fast',20.,0.),('late',20.,T/2)]:
        schedule=AssociationSchedule(T,((delay,T),));t,y,row=nominal(c).run(schedule,tau,PREPARATION_OFFSET)
        _,_,hi=nominal(c,p=2.).run(schedule,tau,PREPARATION_OFFSET)
        row.update(name=name,release=c,delay=delay,nominal_signal_gap_minus_002=hi['acquisition_signal']-row['acquisition_signal']-.02)
        comparisons.append(row);traces[name]=(t,y)
        table(f'trajectory_{name}.csv',['time',*CofactorReactor.columns],zip(t,*y.T))
    configured=nominal().run(AssociationSchedule(T,((0.,T),),float(F(FADING_RATE))),tau,PREPARATION_OFFSET)[2]
    schedules=[]
    for name,windows,gamma in [('constant',((0.,T),),None),('late',((T/2,T),),None),
                               ('pulses',((0.,T/4),(T/2,3*T/4)),None),('fade',((0.,T),),5.),('late-fade',((T/4,T),),1.)]:
        _,_,r=nominal().run(AssociationSchedule(T,windows,gamma),40.,PREPARATION_OFFSET)
        schedules.append(dict(name=name,**r))
    # Equal final product, adjusted exposure; chemistry still determines the cost.
    fixed_product=[]
    for c in [.2,1.,5.,20.]:
        reactor=nominal(c)
        def output(duration):return reactor.run(AssociationSchedule(duration,((0.,duration),)),30.,PREPARATION_OFFSET,samples=31)[2]['complete_product_enclosure_numeric'][0]
        target=.1;right=1.
        while output(right)<target and right<128:right*=2
        if output(right)<target:fixed_product.append(dict(release=c,status='duration-budget-exhausted'));continue
        duration=brentq(lambda d:output(d)-target,1e-5,right,xtol=1e-10)
        result=reactor.run(AssociationSchedule(duration,((0.,duration),)),30.,PREPARATION_OFFSET,samples=101)[2]
        fixed_product.append(dict(release=c,duration=duration,target_product=target,**result))
    box=UNCERTAINTY;general=box.contract(F(ACQUISITION),False,F(PRESERVATION_ALLOWANCE))
    rounded=box.contract(F(ACQUISITION),True,F(PRESERVATION_ALLOWANCE))
    rec=RecoveryCertificate(box);deadlines=[]
    for gamma,n0,n1,wait in [('0.5','0.5','0.5','22.7'),('1','1','1','11.2'),('2','1.85','2','6.8'),('5','1.93','5','5.2'),('10','1.93','10','4.9'),('21','1.93','12','4.9')]:
        gamma,n0,n1,wait=map(F,(gamma,n0,n1,wait));r=rec.fading(gamma,wait,n0,n1);previous=rec.fading(gamma,wait-F(1,10),n0,n1)
        deadlines.append(dict(gamma=gamma,nu0=n0,nu1=n1,wait=wait,**r,certified=r['loss_upper']<F(RECOVERY_TARGET),
                              previous_envelope_above_target=previous['envelope_lower']>F(RECOVERY_TARGET)))
    fading=[];base=nominal().run(AssociationSchedule(T,((0.,T),)),40.,PREPARATION_OFFSET)[2]
    for gamma in [1.,5.,21.]:
        row=nominal().run(AssociationSchedule(T,((0.,T),),gamma),40.,PREPARATION_OFFSET)[2]
        row.update(gamma=gamma,additional_complete_loss_numeric=row['integrated_plus_tail_loss_enclosure_numeric'][0]-base['integrated_plus_tail_loss_enclosure_numeric'][0]);fading.append(row)
    rng=np.random.default_rng(SEED);interior=[]
    for label,prange in [('low',box.low),('high',box.high)]:
        for i in range(INTERIOR_DRAWS_PER_CLASS):
            draw=lambda iv:rng.uniform(float(iv[0]),float(iv[1]))
            p=draw(prange);k=draw(box.native);C=draw(box.pool);u=draw(box.activity);c=draw(box.release);beta=draw(box.dissociation);R=draw(box.reporter);s=rng.uniform(-float(box.preparation),float(box.preparation))
            t,y,row=CofactorReactor(LinearRate(p),LinearRate(k),Reporter(R,c,beta,u),C).run(AssociationSchedule(T,((0.,T),)),tau,s)
            # Both extreme gain/noise records are checked, not a lucky noise draw.
            records=[F(str(row['acquisition_signal']))*box.gain[0]-box.error,F(str(row['acquisition_signal']))*box.gain[1]+box.error]
            interior.append(dict(label=label,parameters=dict(p=p,k=k,C=C,u=u,c=c,beta=beta,R=R,s=s),
                classified_endpoints=[rounded.classify(v,True) for v in records],**row))
    inference=[]
    for y in OBSERVED_RECORDS:
        r=rounded.inverse(F(y));flux=r['flux'];threshold=F(FUNCTIONAL_THRESHOLD)
        decision='incompatible' if flux is None else ('below' if flux[1]<threshold else 'above' if flux[0]>threshold else 'unresolved')
        inference.append(dict(record=F(y),**r,binary_with_promise=rounded.classify(F(y),True),binary_without_promise=rounded.classify(F(y)),functional_decision=decision))
    invrows=[]
    for y in np.linspace(.2,.5,180):
        r=rounded.inverse(F(str(y)))
        if r['p'] is not None:invrows.append([y,*map(float,r['p']),*map(float,r['flux'])])
    table('inference_sweep.csv',['record','p_lower','p_upper','flux_lower','flux_upper'],invrows)
    nonlinear=[]
    for scale in [1,10,100,1000]:
        K=F(10*scale);lo=K*K/(K+1)**2;sandwich=slope_sandwich(lo,1,lo,1,20)
        law=MichaelisMenten(float(K),float(K));model=CofactorReactor(law,law,Reporter(),1.)
        t,y,r=model.run(AssociationSchedule(),40.)
        # This is a long finite horizon. No linear tail formula is applied to MM.
        nonlinear.append(dict(affinity=K,slope_range=(lo,F(1)),cost_sandwich=sandwich,finite_horizon_ratio=float(y[-1,3]/y[-1,2]),final_deficit=float(y[-1,4]-y[-1,0]),evidence='Finite-horizon numerical ratio; infinite-horizon sandwich is a theorem under slope and integrability premises.'))
    alias=[];aliaspaths=[]
    for k,l in [(.5,1.5),(1.5,.5)]:
        t,y,r=nominal(k=k,background=l).run(AssociationSchedule(T,((0.,T),),1.),20.)
        aliaspaths.append(y[:,:3]);alias.append(dict(native=k,background=l,native_stationary_flux=k*REGENERATION*COFACTOR_POOL/(REGENERATION+k+l),signal=r['acquisition_signal']))
    # Native deficits do differ; identical x,b,q alone are the observation alias.
    alias_difference=float(np.max(abs(aliaspaths[0]-aliaspaths[1])))
    def contract_record(c):return dict(rounded=c.rounded,AL=c.AL,AU=c.AU,U=c.U,L=c.L,margin=c.margin,threshold=c.threshold,certified=c.certified)
    summary=dict(constants=box.constants(),box=asdict(box),rounded=contract_record(rounded),unrounded=contract_record(general),
        configured_trajectory=configured,ideal_recovery_upper=rec.ideal(F(RECOVERY)),fast_fading_recovery=rec.fading(21,5,F('1.93'),12),
        comparison=comparisons,schedules=schedules,fixed_product=fixed_product,fading=fading,interior=interior,
        alias=dict(cases=alias,max_observable_difference=alias_difference,note='The source equations depend on native plus background consumption. Native-specific output is not identifiable from these common records.'),
        illustrative_units=dict(acquisition_seconds=F(ACQUISITION)*F(TIME_SCALE_SECONDS),recovery_seconds=F(RECOVERY)*F(TIME_SCALE_SECONDS),
            donor_cap_micromolar=box.high[1]*box.pool[1]*(F(ACQUISITION)+F(RECOVERY))*F(CONCENTRATION_SCALE_MICROMOLAR),
            native_cosubstrate_cap_micromolar=box.native[1]*box.pool[1]*(F(ACQUISITION)+F(RECOVERY))*F(CONCENTRATION_SCALE_MICROMOLAR),
            reporter_cosubstrate_cap_micromolar=box.activity[1]*(1+box.preparation)*box.high[1]*box.pool[1]/(box.high[1]+box.native[0])*(F(ACQUISITION)+F(RECOVERY))*F(CONCENTRATION_SCALE_MICROMOLAR),
            low_complete_loss_picomoles=base['integrated_plus_tail_loss_enclosure_numeric'][0]*float(F(CONCENTRATION_SCALE_MICROMOLAR)*F(VOLUME_MICROLITRES)),
            note='Illustrative scale conversion and finite material caps; not stock recipes or empirical calibration.'))
    dump('summary.json',summary);dump('recovery.json',deadlines);dump('inference.json',inference);dump('nonlinear.json',nonlinear)
    plot(out,traces,invrows,deadlines)
    print(f'Uniform suppression upper: {float(box.constants()["suppression"]):.10%}')
    print(f'Rounded residual detector margin: {float(rounded.margin):.12f}; joint design certified: {rounded.certified}')
    print('Native loss is integrated independently; fading tails include remaining association.',flush=True)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest();here=Path(__file__).parent
    dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,python=platform.python_version(),source_sha256=digest(Path(__file__)),module_sha256={n:digest(here/n) for n in ['reactor.py','certificates.py']},output_sha256={p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}))


def plot(out,traces,inverse,deadlines):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    for name,(t,y) in traces.items():
        axs[0].plot(t,100*(y[:,4]-y[:,0])/y[:,4],label=name)
        axs[1].plot(y[:,2],y[:,3],label=name)
    axs[0].axhline(5,color='black',ls=':');axs[0].set(xlabel='Time including recovery',ylabel='Native-flux suppression (%)',title='Native-flux suppression by reporter design')
    axs[1].set(xlabel='Accumulated reporter product',ylabel='Independently integrated native-output loss',title='Native-output loss versus reporter product')
    for ax in axs:ax.grid(alpha=.2);ax.legend(fontsize=8)
    fig.savefig(out/'storage.png',dpi=180);fig.savefig(out/'storage.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained');iv=np.array(inverse)
    axs[0].fill_between(iv[:,0],iv[:,3],iv[:,4],alpha=.25);axs[0].plot(iv[:,0],iv[:,3]);axs[0].plot(iv[:,0],iv[:,4]);axs[0].axhline(float(F(FUNCTIONAL_THRESHOLD)),color='black',ls=':')
    axs[0].set(xlabel='Bounded-error reporter record',ylabel='Outer interval for stationary native flux',title='Native-flux bounds from reporter\nmeasurements')
    axs[1].plot([float(r['gamma']) for r in deadlines],[float(r['wait']) for r in deadlines],'o-');axs[1].axhline(5,color='black',ls=':')
    axs[1].set(xscale='log',xlabel='Association fade rate',ylabel='Certified recovery wait',title='Required recovery time versus reporter\nshutoff rate')
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'inference-recovery.png',dpi=180);fig.savefig(out/'inference-recovery.svg');plt.close(fig)


if __name__=='__main__':main()