"""From ACR candidate algebra to a usable, conditioned concentration guarantee."""
# EDIT THESE INPUTS. Reactor concentrations: micromolar; time: minutes.
REACTOR_RATES=('1/20','1/20','1/20','1/20','3/20')
DILUTION='1/100'
FEED_CONCENTRATION='10'
NOMINAL_LOAD='1/50'
REQUIRED_B_FLOOR='1'
B_DERIVATIVE_ERROR='1/1000'  # a derivative bound, not concentration measurement noise
RATE_RELATIVE_UNCERTAINTY='1/100'
REACTOR_INITIAL_DISPLACEMENT=('1/50','-1/100','1/100')
TRAJECTORY_MINUTES=600
ENVZ_RATES=('1',)*9
ENVZ_TOTAL_X='3/2'
ENVZ_TOTAL_Y='9/2'
RELEASE_RATES=('1','1')
RELEASE_LEAKAGE=('1/100','1/100')

import argparse,csv,hashlib,json,platform
from pathlib import Path
import numpy as np
import scipy
import sympy as s
from certificates import Interval,ResidualCertificate,BlockCandidates
from models import Reaction,MassActionModel,LoadedReactor,EnvZOmpR,envz_box_floors,PrivateRelease
from invariant_region import InvariantEllipsoid,P
import paper_checks
R=s.Rational
MANUSCRIPT_SHA256='91f82500b08f1947437c95b688e681504e38648ef338bb8ae46dfda58439d776'

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,x):(out/name).write_text(json.dumps(x,indent=2,default=lambda x:x.tolist() if isinstance(x,np.ndarray) else str(x))+'\n')
    def table(name,rows):
        with (out/name).open('w',newline='') as f:w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
    replay=paper_checks.replay();dump('source_arithmetic_replay.json',replay)
    # Original ideals only: no radical, saturation or conservation equations.
    a,b,c,d=s.symbols('a b c d');fields=[(-a*b+b-a*c+3*c,b*(a-2),c*(a-3)+b),(-a*a+c+d-a*b+b-a,d-b*b,a*a-c,c-d)];algebra=[]
    for label,variables,field,point,positive_branch in [('boundary_adds_candidate',(a,b,c),fields[0],(2,1,1),'(2,tau,tau), tau>0: ACR a=2'),('nonpositive_point_supplies_only_candidate',(a,b,c,d),fields[1],(1,1,1,1),'(tau,tau,tau^2,tau^2), tau>0: no ACR')]:
        candidates=BlockCandidates(variables,a).compute(field);algebra.append(dict(name=label,**candidates,positive_geometry=positive_branch,identity_at_positive_point=[str(s.expand(f).subs(dict(zip(variables,point)))) for f in field]))
    u,v,t=s.symbols('u v t');ideal=[u*t-v,v*t-u];mixed=s.groebner(ideal,u,v,t,order=lambda e:(e[0]+e[1],e[2],e[0]));coeff=[]
    for p in mixed.polys:coeff.append(max(s.Poly(p.as_expr(),u,v).terms(),key=lambda term:(sum(term[0]),term[0][0]))[1])
    mixed_roots=sorted({r for c0 in coeff for r in s.Poly(c0,t).real_roots() if r>0});dump('algebraic_candidates_and_positive_geometry.json',algebra);dump('monomial_order_comparison.json',dict(mixed_basis=[str(p.as_expr()) for p in mixed.polys],mixed_coefficients=list(map(str,coeff)),mixed_candidates=list(map(str,mixed_roots)),block=BlockCandidates((u,v,t),t).compute(ideal),scope='The mixed order eliminates u,v but is not a block order. Its empty candidate list misses the actual ACR value 1.'))
    reactor=LoadedReactor(REACTOR_RATES,DILUTION,FEED_CONCENTRATION,NOMINAL_LOAD);model=reactor.model();eq=reactor.equilibrium(REQUIRED_B_FLOOR);dump('configured_reactor_equilibrium.json',eq)
    if eq['status']!='positive_equilibrium':raise ValueError('Configured reactor has no positive operating point. The exact reason is saved; choose admissible parameters for the transient demo.')
    a,b,c=model.variables;k,D,ci,ell=reactor.parameters();alpha=eq['alpha'];cert=ResidualCertificate(model.variables,0,alpha,k[0]*b,(s.Integer(0),s.Integer(1),s.Integer(0)));box=(Interval(0,ci),Interval(R(REQUIRED_B_FLOOR),ci),Interval(0,ci));floor=k[0]*R(REQUIRED_B_FLOOR);eps=(0,R(B_DERIVATIVE_ERROR),0)
    good=cert.evaluate(model.field,box,eps,floor,(-a/50,0,0));dump('reactor_residual_certificate.json',good)
    bad=ResidualCertificate(model.variables,0,alpha,2*k[0]*b,(0,1,0)).evaluate(model.field,box,eps,floor);nofloor=cert.evaluate(model.field,(box[0],Interval(0,ci),box[2]),eps,floor);dump('rejected_and_unresolved_certificates.json',dict(wrong_identity=bad,unsupported_floor=nofloor))
    direct=cert.evaluate(model.field,box,eps,floor,(0,b/100,0));dump('direct_multiplier_load_certificate.json',direct)
    # Disturbances that cancel after q dot d are simplified before enclosure.
    shared=ResidualCertificate((a,b),0,1,b,(1,1));signed=shared.evaluate((b*(a-1),s.Integer(0)),(Interval(0,3),Interval(1,2)),(R(1,1000),R(1,1000)),1,(a*b,-a*b));dump('signed_load_cancellation.json',signed)
    sweep=[]
    for load in [R(j,10000) for j in range(401)]:
        r=LoadedReactor(REACTOR_RATES,DILUTION,FEED_CONCENTRATION,load);e=r.equilibrium(REQUIRED_B_FLOOR);sweep.append(dict(load=float(load),status=e['status'],setpoint=float(e['alpha']),b_pool=float(e['state'][1]) if e['status']=='positive_equilibrium' else '',floor_satisfied=e.get('floor_satisfied',False),local_stability=e.get('locally_stable',False)))
    table('load_operating_envelope.csv',sweep)
    rb=reactor_box(reactor,R(RATE_RELATIVE_UNCERTAINTY));dump('reactor_parameter_box.json',rb)
    ellipsoid=InvariantEllipsoid(model,eq['state']);proof=ellipsoid.verify();dump('invariant_ellipsoid_certificate.json',proof)
    initial=tuple(x+R(v) for x,v in zip(eq['state'],REACTOR_INITIAL_DISPLACEMENT));inside=ellipsoid.contains(initial) if proof['status']=='certified_forward_invariant' else False
    times,Y=model.integrate(initial,TRAJECTORY_MINUTES);_,Y2=model.integrate(initial,TRAJECTORY_MINUTES,method='BDF');center=np.array(eq['state'],float);U=Y-center[:,None];V=np.einsum('it,ij,jt->t',U,np.array(P,float),U);diagnostics=dict(initial_inside_certified_region=inside,solver_difference=float(np.max(abs(Y-Y2))),minimum_b=float(min(Y[1])),initial_V=float(V[0]),maximum_V=float(max(V)),final_error=float(np.linalg.norm(Y[:,-1]-center)),scope='Numerical transient. The fixed-point ellipsoid theorem applies only if its exact checks and initial membership both pass.')
    table('full_reactor_trajectory.csv',[dict(time_min=t,a=Y[0,i],b=Y[1,i],c=Y[2,i],V=V[i]) for i,t in enumerate(times)]);dump('reactor_trajectory_checks.json',diagnostics)
    extinction=[dict(b=str(R(1,10**j)),a_error='1',b_residual=str(k[0]*R(1,10**j))) for j in range(1,10)];dump('small_residual_extinction_counterexample.json',extinction)
    env=EnvZOmpR(ENVZ_RATES);em=env.model();ee=env.equilibrium(ENVZ_TOTAL_X,ENVZ_TOTAL_Y);dump('envz_equilibrium.json',ee)
    if ee['status']=='positive_equilibrium':
        dump('envz_regularity.json',em.regularity(ee['state']));es=np.array([float(v) for v in ee['state']]);ei=es.copy();shift=min(es[0],es[1])/100;ei[0]+=shift;ei[1]-=shift;et,EY=em.integrate(ei,60);poolX=EY[[0,1,2,5,6]].sum(axis=0);poolY=EY[[3,4,5,6]].sum(axis=0);table('envz_full_trajectory.csv',[dict(time_min=t,**{name:EY[j,i] for j,name in enumerate(em.species)}) for i,t in enumerate(et)]);dump('envz_conservation_checks.json',dict(X_total_error=float(max(abs(poolX-float(R(ENVZ_TOTAL_X))))),Y_total_error=float(max(abs(poolY-float(R(ENVZ_TOTAL_Y))))),scope='Numerical trajectory, no global stability assertion.'))
    eb=envz_box_floors((Interval('99/100','101/100'),)*9,Interval(1,2),Interval(4,5));dump('envz_uniform_equilibrium_floors.json',eb)
    x=em.variables;ek,ealpha,*_=env.constants();q1=(0,0,-(ek[7]+ek[8]),ek[7]+ek[8],0,0,ek[8]);q2=(0,0,-ek[6]*x[4],ek[6]*x[4],0,0,ek[2]);domain=tuple(Interval('1/10','5') for _ in x)
    c1=ResidualCertificate(x,4,ealpha,ek[6]*ek[8]*x[1],q1);c2=ResidualCertificate(x,4,ealpha,ek[6]*ek[8]*x[6],q2)
    dump('envz_two_measured_domain_certificates.json',dict(first=c1.evaluate(em.field,domain,(R(1,1000),)*7,ek[6]*ek[8]/10),second=c2.evaluate(em.field,domain,(R(1,1000),)*7,ek[6]*ek[8]/10),scope='A separately stipulated measured-state box, not the equilibrium-only pool floors.'))
    release=release_demo(out,dump,table)
    plot(out,sweep,eq,times,Y,V,release,extinction)
    print('Source symbolic/rational checks:',replay['count'],replay['status'],flush=True)
    print('Candidate lists:',[a['candidates'] for a in algebra],'; mixed/block order:',mixed_roots,'/ [1]',flush=True)
    print('Loaded reactor:',tuple(map(str,eq['state'])),'; floor limit',eq['floor_load_limit'],'; positive-branch limit',eq['critical_load'],flush=True)
    print('Conditional concentration error:',good['absolute_error_bound'],'micromolar; ellipsoid:',proof['status'],'; initial inside:',inside,flush=True)
    print('Equilibrium floors do not automatically hold along transients. Private release preserves steady-state geometry, not trajectories or storage. Illustrative rates; Lean not rerun.',flush=True)
    here=Path(__file__).parent;digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest();dump('run_metadata.json',dict(manuscript_sha256=MANUSCRIPT_SHA256,source_sha256=digest(Path(__file__)),python=platform.python_version(),sympy=s.__version__,scipy=scipy.__version__,numpy=np.__version__,module_sha256={p.name:digest(p) for p in sorted(here.glob('*.py')) if p.name not in ['example.py','test_example.py']},output_sha256={p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}))

def reactor_box(reactor,uncertainty):
    k,D,ci,ell=reactor.parameters()
    if not 0<=uncertainty<1:raise ValueError('Relative uncertainty must lie in [0,1).')
    box=lambda value:Interval(value*(1-uncertainty),value*(1+uncertainty));kb=list(map(box,k));cb=box(ci);lb=box(ell);alpha=(kb[1]+kb[2]+D)/kb[0];delta=kb[4]+D-kb[3]*alpha;N=cb-alpha*(1+lb/D)
    if delta.lo<=0 or N.lo<=0:return dict(status='not_certified_uniform_positive_region',alpha=alpha.strings(),delta=delta.strings(),pool=N.strings())
    r=delta/kb[2];floor=r.lo*N.lo/(1+r.lo);return dict(status='certified_equilibrium_region',alpha=alpha.strings(),delta_lower=str(delta.lo),pool_lower=str(N.lo),ratio_lower=str(r.lo),b_floor=str(floor),nominal_setpoint_uncertainty=str(max(abs(alpha.lo-(k[1]+k[2]+D)/k[0]),abs(alpha.hi-(k[1]+k[2]+D)/k[0]))))

def release_demo(out,dump,table):
    # A+B -> 3B, B -> A, B -> empty. Positive steady states: a=1, b>0.
    old=MassActionModel(('a','b'),[Reaction((1,1),(0,3),1),Reaction((0,1),(1,0),1),Reaction((0,1),(0,0),1)]);release=PrivateRelease(0,RELEASE_RATES);new,W=release.apply(old);a,b,z1,z2=new.variables;F=a*b;lam=tuple(map(R,RELEASE_RATES));reconstruction={z1:F/lam[0],z2:F/lam[1]};defects=[s.simplify(f.subs(reconstruction)) for f in new.field];identity=s.expand((a-1)*b-new.field[1]/2-new.field[2]*R(3,2)-new.field[3])
    if any(s.expand(a-b)!=0 for a,b in zip(defects[:2],old.field)) or any(defects[2:]) or identity!=0:raise ArithmeticError('Private-release quotient or transported identity failed.')
    dump('release_equilibrium_and_residual_identity.json',dict(old_field=list(map(str,old.field)),extended_field=list(map(str,new.field)),old_equations_recovered=True,transported_identity='(a-1)b = new_b_derivative/2 + 3*private_1_derivative/2 + private_2_derivative',quotient_scope='Triangular elimination in the original steady-state ideal; no fast-release approximation.',zero_divisor_candidates_old=BlockCandidates(old.variables,old.variables[0]).compute(old.field),zero_divisor_candidates_extended=BlockCandidates(new.variables,new.variables[0]).compute(new.field)))
    costs=[]
    for l in map(R,['1/100','1/10','1','10']):
        chain=PrivateRelease(0,(l,l));c=chain.steady_costs(1,weights=(3,2));costs.append(dict(release_rate=float(l),unit_weight_storage=float(chain.steady_costs(1)['weighted_storage']),cargo_weighted_storage=float(c['weighted_storage'])))
    table('release_storage_sweep.csv',costs);leak=PrivateRelease(0,RELEASE_RATES,RELEASE_LEAKAGE);leak_model,LW=leak.apply(old);lc=leak.steady_costs(1);dump('release_leakage.json',lc)
    # At a positive full steady state a'=b(1-a) forces a=1. Any nonzero
    # leak makes the b-balance strictly negative: there is no positive state.
    yields=lc['product_yields'];effective=R(1,2)*((1-yields[0])+2*(1-yields[1]));dump('release_leakage_existence_check.json',dict(certificate_defect_per_input_flux=effective,positive_equilibrium_exists=bool(effective==0),reason='a balance forces a=1; b balance is -b times a strictly positive release loss when any stage leaks. This example demonstrates why a conditional leakage error bound is not an existence proof.'))
    tt,YY=new.integrate((R(4,5),2,R(1,10),R(1,10)),60);ledger=2*YY[0]+YY[1]+3*YY[2]+2*YY[3];dump('release_trajectory_checks.json',dict(weighted_balance_error=float(max(abs(ledger-ledger[0]))),scope='For these chosen equal B-to-A and B-loss rates only, the weighted field balance cancels. This is not an asserted general conservation law under release.'));table('release_full_trajectory.csv',[dict(time_min=t,**{label:YY[j,i] for j,label in enumerate(new.species)}) for i,t in enumerate(tt)])
    return costs

def plot(out,sweep,eq,times,Y,V,release,extinction):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(11,4),layout='constrained');valid=[r for r in sweep if r['status']=='positive_equilibrium']
    axs[0].plot([r['load'] for r in valid],[r['b_pool'] for r in valid],label='Equilibrium multiplier pool b');axs[0].axhline(float(R(REQUIRED_B_FLOOR)),color='C1',ls='--',label='Required floor');axs[0].axvline(float(eq['floor_load_limit']),color='C1',ls=':');axs[0].axvline(float(eq['critical_load']),color='C3',ls=':');axs[0].set(xlabel='Specific load / per minute',ylabel='Concentration / micromolar',title='Supporting pool versus reactor load');axs[0].legend(fontsize=8)
    for i,name in enumerate(['a','b','c']):axs[1].plot(times,Y[i],label=name)
    axs[1].set(xlabel='Time / minutes',ylabel='Concentration / micromolar',title='Full loaded-reactor trajectory');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'operating_region_recovery.png',dpi=180);fig.savefig(out/'operating_region_recovery.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(11,4),layout='constrained');axs[0].loglog([float(R(r['b'])) for r in extinction],[float(R(r['b_residual'])) for r in extinction],'o-',label='Residual magnitude');axs[0].loglog([float(R(r['b'])) for r in extinction],[1]*len(extinction),'--',label='Concentration error stays 1');axs[0].set(xlabel='Multiplier pool b / micromolar',ylabel='Value (respective concentration/time units)',title='Rate residual and concentration error as B\nvanishes');axs[0].legend(fontsize=8)
    axs[1].loglog([r['release_rate'] for r in release],[r['unit_weight_storage'] for r in release],'o-',label='Intermediate concentration');axs[1].loglog([r['release_rate'] for r in release],[r['cargo_weighted_storage'] for r in release],'s-',label='Cargo-weighted inventory');axs[1].set(xlabel='Release rate / per minute',ylabel='Storage at input flux 1 micromolar/min',title='Intermediate storage versus product-release\nrate');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'conditioning_and_storage.png',dpi=180);fig.savefig(out/'conditioning_and_storage.svg');plt.close(fig)

if __name__=='__main__':main()