Example code

"""Four stable labels are not yet an economical memory. Run: python example.py"""
# EDITABLE INPUTS. Exact paper certificates remain separate from the exploratory reactor.
SITES = 3
# One six-rate tuple per sequential level: E binding, release, catalytic; F binding, release, catalytic.
# These moderate, illustrative values define a symmetric cube lift, not the four-sink source.
CHAIN_RATES = [('1','1','1','1','1','1')]*3
TOTAL_KINASE = '2'
TOTAL_PHOSPHATASE = '2'
TOTAL_SUBSTRATE = '20'
INITIAL_PHOSPHOFORM = 0
HORIZON = 50.
REVERSE_ACTIVITY = '0'           # ADP = phosphate activity; zero is irreversible reference
CLOCK_MULTIPLIER = '1'           # applied consistently to all elementary rates
RUN_SSA = True
SYSTEM_SIZE = 20                # count = system size * concentration
SSA_SEED = 72
SSA_EVENT_BUDGET = 200000
# Resource exploration uses the verified four-sink source and its own model clock.
CERTIFIED_SIZE_MULTIPLIER = '1'
CERTIFIED_CLOCK_MULTIPLIER = '1'
LOADING_KERNEL_SHIFT = '0'      # separate family exploration, preserves fixed target/effective rates
OBSERVATION_ERROR_FRACTION = '0.01'

from pathlib import Path
from fractions import Fraction as F
import argparse,csv,hashlib,json,math,platform,sys
import numpy as np
import mpmath as mp
import sympy as sp
from reactor import EdgeKinetics,PhosphorylationNetwork
from geometry import EquilibriumGeometry,capacity_bounds
from certificates import ExactPaperReplay,PublishedSource,PublishedOperatingContract,driving_ratios,load
import paper_checks as cp

MANUSCRIPT_SHA256='7970bdb2f32090091d3a3a5dc2bc73e835db7b7d23165410ff053eb5aa2175ea'


def serialize(value):
    if isinstance(value,F):return str(value)
    if isinstance(value,np.ndarray):return value.tolist()
    if isinstance(value,np.generic):return value.item()
    if isinstance(value,mp.mpf):return mp.nstr(value,60)
    raise TypeError(type(value).__name__)


def main():
    sys.set_int_max_str_digits(0)
    ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('--output',default='outputs');args=ap.parse_args();out=Path(args.output);out.mkdir(exist_ok=True,parents=True)
    def dump(name,value):(out/name).write_text(json.dumps(value,indent=2,default=serialize)+'\n',encoding='utf-8')
    def table(name,headers,rows):
        with (out/name).open('w',newline='',encoding='utf-8') as f:
            w=csv.writer(f);w.writerow(headers);w.writerows(rows)
    print('Replaying the original and final source: roots, full-class stability, operating contracts and integer preparations...',flush=True)
    replay=ExactPaperReplay().run();dump('fresh_exact_replay.json',replay)
    assert replay['admissible_equilibria']==7 and len(replay['final_sinks'])==4
    assert replay['final_worst_retention_bound']<.00105 and replay['final_worst_recovery_bound']<.008987
    source=PublishedSource();net=source.network
    tw=load('tree_weights.json')['weights_descending_coefficients'];u=sp.Symbol('u')
    tau=[sp.Poly.from_list(row,u).as_expr() for row in tw]
    geometry=EquilibriumGeometry(net,tau);roots=geometry.census(source.totals)
    assert len(roots)==8 and sum(r['admissible'] for r in roots)==7
    assert sp.expand(geometry.B-10*geometry.D-sp.prod(u-j for j in range(1,7)))==0
    dump('independent_root_census.json',dict(eliminant_coefficients=[str(x) for x in geometry.eliminant(source.totals).all_coeffs()],roots=roots,scope='All roots isolated exactly; L and M have no roots in each regular isolating interval. u=r is handled separately by the reusable geometry. Simple roots plus four interval-Hurwitz certificates and the paper degree theorem imply the other three have index -1.'))
    states=[];sink_intervals=[list(map(F,c['root_interval'])) for c in load('postproof_candidate_sinks.json')['sink_certificates']]
    with mp.workdps(80):
        for row in roots:
            if not row['admissible']:continue
            ratio=sum(row['interval'])/2;x=geometry.state_on_curve(ratio,source.totals[:2],digits=80)
            idx=len(states);sink=any(max(row['interval'][0],lo)<=min(row['interval'][1],hi) for lo,hi in sink_intervals)
            read=x[7]+sum(x[22+i] for i,(a,b) in enumerate(net.edges) if b==7)
            totals=[sum(mp.mpf(int(v))*x[j] for j,v in enumerate(r)) for r in net.inventory]
            states.append(dict(index=idx,ratio=float(ratio),sink_certified=sink,normalized_readout=float(read/cp.ST),state=[mp.nstr(v,65) for v in x],total_residual=[mp.nstr(t-mp.mpf(str(expected)),10) for t,expected in zip(totals,source.totals)]))
    dump('equilibrium_states.json',dict(states=states,scope='80-digit midpoint reconstructions for reuse and plotting. Stability comes from interval certificates, not floating-point spectra. Bound substrate is retained in the total.'))
    squares=net.square_ratios();assert any(r[-1]!=1 for r in squares)
    dump('cycle_structure.json',dict(squares=[dict(tail=v,sites=[i,j],path_ratio=r,log_effective_affinity=math.log(float(r))) for v,i,j,r in squares],
        detailed_balance_ratios=driving_ratios(net),scope='Effective pathway affinity differs from the common ATP-hydrolysis potential. Nonzero square affinity is necessary, not sufficient, for more than n sinks.'))
    # A fixed-target family direction: compare states at the same free-enzyme ratio, not at the same total.
    shifted=PublishedSource(F(LOADING_KERNEL_SHIFT));sg=EquilibriumGeometry(shifted.network,tau)
    assert sp.expand(sg.B-10*sg.D-(geometry.B-10*geometry.D))==0
    position=F(83,20);x0=geometry.state_on_curve(position,source.totals[:2]);x1=sg.state_on_curve(position,source.totals[:2])
    lower=F(2260898025826008319,15168384)
    assert F(10)/lower<F(68,10**12)
    dump('loading_family_obstruction.json',dict(shift=F(LOADING_KERNEL_SHIFT),ratio=position,free_substrate_reference=str(sum(x0[:8])),free_substrate_shifted=str(sum(x1[:8])),exact_free_substrate_lower_on_4_1_to_4_2=lower,kinase_to_substrate_upper=F(10)/lower,
        scope='Fixed effective rates, target L, enzyme totals and retained ratio interval. A kernel shift preserves free substrate at a fixed ratio but generally changes the compatibility-class total and its equilibrium locations. No stability certificate is transferred to the shifted source.'))
    operating=PublishedOperatingContract().calculate(F(CERTIFIED_SIZE_MULTIPLIER),F(CERTIFIED_CLOCK_MULTIPLIER));dump('configured_operating_contract.json',operating)
    baseline=PublishedOperatingContract().calculate();clockcheck=PublishedOperatingContract().calculate(clock_multiplier=F(1,100))
    assert baseline['required_size']==clockcheck['required_size'] and baseline['rows']==clockcheck['rows']
    dump('common_clock_identity.json',dict(clock_multiplier='1/100',budget_unchanged=True,failure_bounds_unchanged=True,time_multiplier=100,scope='Matched recovery and storage horizons. A fixed absolute horizon is a different contract.'))
    caps=[capacity_bounds(n) for n in range(2,13)];dump('capacity_bounds.json',caps)
    gap=min(states[j+2]['normalized_readout']-states[j]['normalized_readout'] for j in [0,2,4]);obs=F(OBSERVATION_ERROR_FRACTION)
    if obs<=0:raise ValueError('positive observation tolerance required')
    dump('observation_and_crowding.json',dict(normalized_minimum_sink_gap=gap,configured_error=obs,ideal_center_error_boxes_disjoint=2*float(obs)<gap,scalar_packing_upper=1+math.floor(1/(2*obs)),
        source_guaranteed_driving_activity='1/10^102',common_driving_kBT=204*math.log(10),relative_rate_error_certificate='1/10^99',
        scope='Ideal-center separation alone is not an operational guarantee: safe-set variation also belongs in decoder boxes. The fresh paper replay checks the complete source decoder. Driving preserves four deterministic sinks only; its finite-molecule contract is not recomputed. Crowding scaling applies to a different sequential split family and an iterated limit.'))
    crowd=[]
    for delta in [.5,.25,.125,.0625]:crowd.append((delta,delta**6,delta**-4,delta))
    table('crowding_scaling_reference.csv',['separation_delta','relative_barrier_delta6','relative_recovery_delta_minus4','relative_readout_gap_delta'],crowd)
    # A tractable full-network example for editing and simulation; no four-label claim.
    if len(CHAIN_RATES)!=SITES:raise ValueError('one six-rate tuple per site is required')
    chain_rates=[EdgeKinetics(*(F(v) for v in row)).scaled(F(CLOCK_MULTIPLIER)) for row in CHAIN_RATES]
    model=PhosphorylationNetwork.symmetric_lift(chain_rates,F(REVERSE_ACTIVITY));chain=PhosphorylationNetwork.chain(chain_rates,F(REVERSE_ACTIVITY));A=model.cube_aggregation()
    if type(INITIAL_PHOSPHOFORM) is not int or not 0<=INITIAL_PHOSPHOFORM<model.q:raise ValueError('invalid initial phosphoform')
    totals=list(map(F,[TOTAL_KINASE,TOTAL_PHOSPHATASE,TOTAL_SUBSTRATE]))
    if min(totals)<=0:raise ValueError('positive inventories required')
    initial=np.zeros(model.size);initial[model.E]=float(totals[0]);initial[model.F]=float(totals[1]);initial[INITIAL_PHOSPHOFORM]=float(totals[2])
    times=np.linspace(0,HORIZON,251);trajectory=model.simulate(initial,times);[email protected]
    # The paper's exact same-rate chain reduction is for the irreversible lift.
    # Adding the thermodynamic reverse channels changes their aggregation multiplicities.
    reference=chain.simulate(A@initial,times) if not F(REVERSE_ACTIVITY) else None
    residual=float(np.max(abs(aggregated-reference))) if reference is not None else None
    ledger=float(np.max(abs([email protected](totals,float))))
    table('configured_full_trajectory.csv',['time']+model.names,[(t,*x) for t,x in zip(times,trajectory)])
    table('configured_reactions.csv',['reaction','rate','reactants','products'],[(r.name,str(r.rate),'+'.join(model.names[i] for i in r.reactants),'+'.join(model.names[i] for i in r.products)) for r in model.reactions])
    ssa=None
    if RUN_SSA:
        count_totals=[t*SYSTEM_SIZE for t in totals]
        if any(t.denominator!=1 for t in count_totals):raise ValueError('SSA size must make all conserved totals integral')
        counts=[0]*model.size;counts[model.E]=int(count_totals[0]);counts[model.F]=int(count_totals[1]);counts[INITIAL_PHOSPHOFORM]=int(count_totals[2])
        st,sx,events=model.ssa(counts,SYSTEM_SIZE,HORIZON,SSA_SEED,SSA_EVENT_BUDGET)
        for x in sx:assert all(a==b for a,b in zip(model.inventory@x,count_totals))
        stride=max(1,len(st)//1000);selection=sorted(set(range(0,len(st),stride))|{len(st)-1})
        table('configured_ssa_sampled.csv',['time']+model.names,[(st[i],*sx[i]) for i in selection])
        ssa=dict(events=len(events),seed=SSA_SEED,system_size=SYSTEM_SIZE,exact_inventories_preserved=True,completed_horizon=HORIZON);dump('configured_ssa.json',ssa)
    dump('configured_reactor.json',dict(sites=SITES,species=model.size,reactions=len(model.reactions),totals=totals,reverse_activity=F(REVERSE_ACTIVITY),aggregation_max_residual=residual,inventory_max_residual=ledger,
        scope='Moderate illustrative symmetric lift, not the four-label source. Full retained-complex ODE and literal finite-count CTMC. Same-rate chain comparison applies only at zero reverse activity; its residual is null otherwise. One path does not estimate a rare-event retention probability.'))
    curve=[]
    for left,right in [(1e-9,.999999),(2.000001,2.999999),(4.000001,4.999999),(6.000001,9.999999)]:
        grid=np.geomspace(left,right,160) if left<1 else np.linspace(left,right,160)
        for z in grid:
            x=geometry.state_on_curve(str(z),source.totals[:2],digits=50);total=sum(x[:8])+sum(x[10:]);curve.append((left,z,float(total)))
    table('equilibrium_curve.csv',['branch_start','free_enzyme_ratio','retained_substrate_total'],curve)
    plot(out,states,curve,caps,replay,times,trajectory,model)
    print(f'Seven positive equilibria, four exact full-class sinks; minimum normalized readout gap {gap:.6f}.',flush=True)
    print(f'Final sufficient log10(substrate molecules) {replay["physical_final_local"]["log10_substrate"]:.3f}; log10(recovery seconds) {replay["physical_final_local"]["log10_trec_s"]:.3f}.',flush=True)
    print(f'Illustrative full-network aggregation residual {residual}; inventory residual {ledger:.3g}; SSA {ssa}.',flush=True)
    print('Sufficient resources are not minimum requirements. No writing protocol or finite-count reversible guarantee is claimed; 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(),module_sha256={p.name:digest(p) for p in sorted(here.glob('*.py')) if p.name not in ['example.py','test_example.py']},input_sha256={p.name:digest(p) for p in sorted(here.glob('*.json')) if p.name!='release.json'},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,states,curve,caps,replay,times,trajectory,model):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,3,figsize=(12,4),layout='constrained')
    for left in sorted(set(r[0] for r in curve)):
        rows=[r for r in curve if r[0]==left];axs[0].semilogy([r[1] for r in rows],[r[2] for r in rows],color='C0')
    for row in states:axs[0].plot(row['ratio'],cp.ST,'o' if row['sink_certified'] else 'x',color='C1' if row['sink_certified'] else 'black')
    axs[0].axhline(cp.ST,color='gray',ls=':');axs[0].set(xlabel='Free-enzyme ratio E/F',ylabel='Retained substrate total',title='Seven equilibria, four sinks',ylim=(1e6,1e16))
    sinks=[r for r in states if r['sink_certified']]
    axs[1].errorbar([r['normalized_readout'] for r in sinks],range(1,5),xerr=replay['final_min_normalized_gap']/8,fmt='o',capsize=4)
    axs[1].set(xlabel='Fully phosphorylated retained fraction',ylabel='Stable label',title='Phosphorylation readouts of the four stable\nlabels',yticks=[1,2,3,4],xlim=(-.04,.75))
    ns=[r['sites'] for r in caps];axs[2].fill_between(ns,[r['constructive_lower'] for r in caps],[r['universal_upper'] for r in caps],alpha=.2,label='Proved interval')
    axs[2].plot(ns,[r['square_balanced_capacity'] for r in caps],label='Square balanced: exact n');axs[2].set(yscale='log',xlabel='Sites',ylabel='Stable-state capacity',title='Constructive and upper bounds on\nstable-state count');axs[2].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'structure_and_labels.png',dpi=180);fig.savefig(out/'structure_and_labels.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    for key,label in [('physical_original_global','Original / global noise'),('physical_original_local','Original / local noise'),('physical_final_local','Final / local noise')]:
        r=replay[key];axs[0].scatter(r['log10_trec_s'],r['log10_substrate'],label=label)
    axs[0].set(xlabel='log10 certified recovery time (s)',ylabel='log10 sufficient substrate molecules',title='Sufficient molecule counts and recovery\ntimes');axs[0].legend(fontsize=8)
    for level in range(SITES+1):
        cols=[v for v in range(model.q) if v.bit_count()==level];axs[1].plot(times,trajectory[:,cols].sum(axis=1),label=f'{level} modified sites')
    axs[1].set(xlabel='Time (illustrative units)',ylabel='Free substrate by level',title='Separate editable symmetric reactor');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'resources_and_reactor.png',dpi=180);fig.savefig(out/'resources_and_reactor.svg');plt.close(fig)


if __name__=='__main__':main()