"""Energy service, peroxide fates, finite carrier kinetics, and paired signals. Run python example.py --output outputs. All S7 certificates are checked afresh. Illustrative kinetic coefficients are NOT fitted red-cell parameters. """ from __future__ import annotations import argparse import csv from dataclasses import dataclass from fractions import Fraction as Q import hashlib import json import math from pathlib import Path import platform import numpy as np from scipy.integrate import solve_ivp import sympy as sp import mpmath as mp from metabolic_model import SourceModel,RationalAudit,LinearRelaxation,OBJECTIVES # EDITABLE INPUTS ----------------------------------------------------------- HORIZON_HOURS = 1 REQUIRED_SERVICE = 1.0 # mmol/gDW over horizon; source Na/K ATPase MEDIUM = 'restricted' # original, sulfur_closed, restricted, internal_blocked INITIAL_CHEMICAL_AMOUNTS = {} # zero is the paper's algebraic diagnostic, not a viable cell TERMINAL_CHEMICAL_FLOORS = {} OXIDATION_CAP = 2.0 # a, per hour; illustrative carrier module REGENERATION_CAP = 1.0 # b, per hour; illustrative TOTAL_CARRIER = 1.0 # arbitrary unit; results per unit, no invented physiological pool INITIAL_REDUCED_FRACTIONS = (0., 1/3, 1.) KINETIC_HORIZON = 3.0 PAIRED_PEAK = 1.0 # synthetic raw signal units RESIDUAL_INPUT = .03 # one constant r at every time SIGNAL_HORIZON = 4.0 RAW_OBSERVATIONS = ((Q('0.8'),Q('0.5')),(Q('0.7'),Q('0.6'))) # (untreated, treated) RAW_ERROR = Q('0.01') RESIDUAL_INPUT_INTERVAL = (Q(0),Q('0.05')) LP_TOLERANCE = 1e-8 # numerical objective-fixing tolerance, not certification RTOL, ATOL = 1e-10, 1e-12 MANUSCRIPT_SHA256 = 'a8cd8efce97007780e7b88fc7f19dc94a31c3e7da20f22b66dbc3f1a0db97774' # -------------------------------------------------------------------------- class PlateauCertificate: """One exact feasible witness plus two global bounds covers an entire interval.""" def __init__(self,primal,turnover_upper,service_upper): self.M=primal['M'];self.L=primal['Q'];self.U=turnover_upper['bound_without_reverse'];self.UM=service_upper['bound_without_reverse'] if not 0<=self.M<=self.UM or not 0<=self.L<=self.U:raise ValueError('Inconsistent enclosure inputs.') def query(self,required): required=Q(required) if required<0:raise ValueError('Nonnegative required service expected.') if required>self.UM:return dict(status='infeasible_by_exact_service_bound') if required<=self.M:return dict(status='uniformly_enclosed',lower=self.L,upper=self.U) return dict(status='unresolved_feasibility',global_turnover_upper=self.U) @dataclass(frozen=True) class CarrierPool: oxidation: float = OXIDATION_CAP regeneration: float = REGENERATION_CAP total: float = TOTAL_CARRIER def __post_init__(self): if not all(np.isfinite(v) for v in (self.oxidation,self.regeneration,self.total)) or min(self.oxidation,self.regeneration)<=0 or self.total<0:raise ValueError('Positive finite a,b and nonnegative finite total required.') def validate(self,g0,time): if not np.isfinite(g0) or not 0<=g0<=self.total or not np.isfinite(time) or time<0:raise ValueError('Initial reduced pool and nonnegative horizon out of range.') def bound(self,time,g0): self.validate(g0,time);a,b,C=self.oxidation,self.regeneration,self.total;k=a+b;s=k*time # T - (1-exp(-kT))/k loses its leading quadratic term at short times. if abs(s)<1e-3: startup=sum((-1)**n*s**n/math.factorial(n) for n in range(2,13))/k else:startup=time+math.expm1(-s)/k return a/k*(-math.expm1(-s))*g0+a*b*C/k*startup def saturated_state(self,time,g0): self.validate(g0,time);a,b,C=self.oxidation,self.regeneration,self.total return b*C/(a+b)+(g0-b*C/(a+b))*math.exp(-(a+b)*time) @property def harmonic_rate(self):return self.total/(1/self.oxidation+1/self.regeneration) def required_total(self,turnover,time,theta): if not np.isfinite(turnover) or turnover<0 or not np.isfinite(theta) or not 0<=theta<=1:raise ValueError('Nonnegative turnover and reduced fraction in [0,1] required.') per_unit=CarrierPool(self.oxidation,self.regeneration,1).bound(time,theta) if per_unit==0:return 0. if turnover==0 else None return turnover/per_unit def evolve(self,g0,duration,control=None): self.validate(g0,duration) if duration<=0:raise ValueError('Use bound() for the zero-time endpoint.') def rhs(t,v): oxidation,regeneration=(1.,1.) if control is None else control(t,float(v[0])) if not 0<=oxidation<=1 or not 0<=regeneration<=1:raise ValueError('Control fractions outside rate envelopes.') g=v[0];vout=oxidation*self.oxidation*g;uin=regeneration*self.regeneration*(self.total-g) return [uin-vout,vout,uin] t=np.linspace(0,duration,301);sol=solve_ivp(rhs,(0,duration),[g0,0,0],t_eval=t,rtol=RTOL,atol=ATOL,method='DOP853') if not sol.success:raise RuntimeError(sol.message) return t,sol.y.T class PairedRecovery: """Same initial peak, one residual-input fraction, common effective reduction.""" def __init__(self,peak=PAIRED_PEAK,residual=RESIDUAL_INPUT): if not np.isfinite(peak) or peak<=0 or not np.isfinite(residual) or not 0<=residual<1:raise ValueError('Positive common peak and r in [0,1) required.') self.peak=peak;self.r=residual def transform(self,untreated,treated):return (np.asarray(treated)-self.r*np.asarray(untreated))/((1-self.r)*self.peak) def evolve(self,duration,input_rate,reduction,discrepancy=None): def rhs(t,v): q=input_rate(t);k=reduction(t);extra=0 if discrepancy is None else discrepancy(t,v[1]) if q<0 or k<0:raise ValueError('Synthetic input and reduction rates must be nonnegative.') return [q-k*v[0],self.r*q-k*v[1]+extra,k,abs(extra)/((1-self.r)*self.peak)] t=np.linspace(0,duration,401);sol=solve_ivp(rhs,(0,duration),[self.peak,self.peak,0,0],t_eval=t,rtol=RTOL,atol=ATOL,method='DOP853') if not sol.success:raise RuntimeError(sol.message) return t,sol.y.T @staticmethod def transformed_error(r,peak,error_C,error_D,error_peak,Zmax=1): r,peak,ec,ed,ep,z=map(Q,(r,peak,error_C,error_D,error_peak,Zmax)) if not 0<=r<1 or peak<=ep or min(ec,ed,ep,z)<0:raise ValueError('Invalid raw-error assumptions.') return (ec+r*ed)/((1-r)*(peak-ep))+z*ep/(peak-ep) @staticmethod def infer_integrated(z_observed,total_error): z,e=Q(z_observed),Q(total_error) if e<0:raise ValueError('Nonnegative discrepancy plus observation allowance required.') lo=z-e;hi=min(Q(1),z+e) if hi<=0 or lo>hi:return dict(status='incompatible') mp.iv.dps=45 iv=lambda v:mp.iv.mpf(v.numerator)/v.denominator lower=math.nextafter(float((-mp.iv.log(iv(hi))).a),-math.inf) if lo<=0:return dict(status='no_finite_upper_identification',K_lower=max(0.,lower),K_upper=None) upper=math.nextafter(float((-mp.iv.log(iv(lo))).b),math.inf) return dict(status='bounded',K_lower=max(0.,lower),K_upper=upper,z_lower=str(lo),z_upper=str(hi)) @staticmethod def reject_increase(first,second,error,r_interval,peak=1,discrepancy_sum=0): """A sufficient exact rejection with ONE common r and peak. discrepancy_sum bounds the sum of the two normalized intervention errors. It is subtracted before declaring rejection; passing is only consistency. """ d1,c1=map(Q,first);d2,c2=map(Q,second);e=Q(error);rlo,rhi=map(Q,r_interval);P=Q(peak);gamma=Q(discrepancy_sum) if e<0 or not 0<=rlo<=rhi<1 or P<=0 or gamma<0:raise ValueError('Invalid common-r observation assumptions.') margins=[c2-c1-2*e+r*(d1-d2-2*e)-(1-r)*P*gamma for r in (rlo,rhi)] return dict(rejected=min(margins)>0,minimum_numerator_margin=min(margins),common_r=list(map(str,(rlo,rhi)))) def symbolic_checks(): a,b,C,T,t,g0=sp.symbols('a b C T t g0',positive=True);k=a+b weight=a/k*(1-sp.exp(-k*(T-t)));g=b*C/k+(g0-b*C/k)*sp.exp(-k*t) bound=a*b*C*T/k+a/k*(g0-b*C/k)*(1-sp.exp(-k*T)) assert sp.simplify(sp.diff(weight,t)-(k*weight-a))==0 and weight.subs(t,T)==0 assert sp.simplify(sp.diff(g,t)-b*(C-g)+a*g)==0 assert sp.simplify(sp.integrate(a*g,(t,0,T))-bound)==0 assert sp.simplify(sp.limit(bound/T,T,sp.oo)-a*b*C/k)==0 quadratic=sp.expand(sp.series(bound.subs(g0,0),T,0,3).removeO()) assert sp.simplify(quadratic-a*b*C*T*T/2)==0 D,Cs,q,r,kk=sp.symbols('D Cs q r kk') assert sp.expand((r*q-kk*Cs)-r*(q-kk*D)+kk*(Cs-r*D))==0 return dict(adjoint=str(weight),attaining_state=str(g),finite_time_bound=str(bound),oxidized_start=str(quadratic),paired_cancellation='(Cs-rD)prime = -k(Cs-rD) exactly') def run_lp_diagnostics(model): lp=LinearRelaxation(model,MEDIUM,HORIZON_HOURS,INITIAL_CHEMICAL_AMOUNTS,TERMINAL_CHEMICAL_FLOORS) baseline=lp.solve(required_service=REQUIRED_SERVICE) if not baseline['finished']:return dict(baseline=baseline,scope='No successful optimum; later cleanup omitted.') fixed={name:(baseline[key]-LP_TOLERANCE,baseline[key]+LP_TOLERANCE) for name,key in [('service','M'),('turnover','Q')]} cleaned=lp.solve('peroxide',REQUIRED_SERVICE,fixed) parsimonious=None if cleaned['finished']: fixed['peroxide']=(cleaned['terminal_peroxide']-lp.initial[model.ri['M_h2o2_c']]-LP_TOLERANCE,cleaned['terminal_peroxide']-lp.initial[model.ri['M_h2o2_c']]+LP_TOLERANCE) parsimonious=lp.solve('turnover',REQUIRED_SERVICE,fixed,parsimony=True) sensitivity=[dict(fraction=f,result=lp.solve(required_service=REQUIRED_SERVICE,capacity_fractions={'R_GTHOx':f})) for f in (0.,.5,1.)] return dict(baseline=baseline,minimum_terminal=cleaned,parsimonious=parsimonious,GTHOx_sensitivity=sensitivity, scope='Fresh floating-point LPs on the configured medium; intermediate-burden and damage constraints absent. Absolute-flux parsimony mixes source column scales and is a declared selection rule.') def main(): ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('--output',type=Path,default=Path('outputs'));ap.add_argument('--skip-lp',action='store_true');args=ap.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True) source=Path(__file__).parent;model,certs,manifest=SourceModel.load(source);audit=RationalAudit(model) primal=audit.primal(certs['primal']);duals={name:audit.upper(certs[name],OBJECTIVES[name],medium='original' if name=='joint' else 'restricted') for name in ('joint','turnover','service')} assert duals['joint']['interval_cost']==Q(1283461673726143619667945455169487587,50452437500000000000000000000000000)=endpoint and plateau.L>=Q('6.50048030') and plateau.U<=Q('6.50048050') assert plateau.U-plateau.L