Example code — Exact common-activity compatibility of autocatalytic cores: finite cycle acceleration, monotone elimination, and fixed-parameter tractability in the longest path

Autocatalytic cores that share species must produce their surpluses at the same effective concentrations, or chemical activities. Each graph edge represents a two-reaction core. Both internal species must be produced by every core separately, using one common activity for each graph vertex. Summing production across neighboring cores answers a weaker question.

This example provides reusable core, activity-box, uncertainty and solver components. Its fixed-margin solver propagates forced lower bounds and jumps to the next admissible cycle root when ordinary iteration would approach a solution without reaching it exactly. Exact traces retain the root and simple path that produced each value; an independent algebraic query checks that no smaller feasible activity vector exists.

A forced cycle root replaces endless ordinary iteration; two cycle-root branches meet at the exact unit-core capacity of one forty-eighth. — The jump and least-state identities are exact. The convergence curve displays ordinary floating-point iteration for illustration; a small residual is not an exact stopping certificate.

Every edge of the weighted triangle produces both endpoints, while replication around a shared hub increases aggregate food demand. — The bars use exact nominal residuals. Common factor scaling keeps total food demand constant but lowers each individual production residual; longest-path complexity and resource cost are different quantities.

For the unit core at margin 0.01, one exact jump reaches the irrational least state. Further examples show why a tight selected subsystem can choose the wrong root, why zero-margin feasibility need not imply strict productivity, and why disconnected feasible sets require a next-admissible-value map.

The weighted triangle has a rational productive state, a sharp fixed-state factor tolerance of 19/301, and a joint activity/factor uncertainty certificate. Replicating it around one shared hub preserves the certificate at every size, but total food demand grows with the number of modules. Dividing every factor by that number controls aggregate demand while reducing each core's production rate.

The working tools use general exact nonlinear-real-arithmetic queries. They do not implement the paper's projected-CAD, positive-infinitesimal or complexity-guaranteed backend. A separate exact grid-elimination demonstration preserves disconnected sets; grid infeasibility makes no claim about real feasibility. These are static compatibility certificates, not steady-state or persistence guarantees. Lean is not rerun.

Python source

"""One shared activity per species; finite exact cycle jumps and reusable compatibility checks."""
# EDITABLE RATIONAL INPUTS. Activities are normalized; factors are illustrative, not measured.
MODULES = 3
PATH_FACTOR_A = '2'
PATH_FACTOR_B = '1'
SHORTCUT_FACTOR_A = '1'
HUB_BOX = ('0.1','0.9')
MIDDLE_BOX = ('0.01','0.9')
END_BOX = ('0.001','0.9')
FIXED_MARGIN = '0.001'
RATIO_RELATIVE_FACTOR_ERROR = '0.05'  # independent relative a,b errors -> worst-case ratio interval
NOMINAL_ACTIVITIES = ('0.1','0.065','0.043')
ACTIVITY_ERROR = '0.0001'
GRID_DENOMINATOR = 256
SOLVER_TIMEOUT_MS = 30000
MAX_CYCLE_JUMPS = 200
CAPACITY_BRACKET_BITS = 12
REFERENCE_CONCENTRATION_MOLAR = .001
REFERENCE_FLUX_MOLAR_PER_MINUTE = .001

from pathlib import Path
from dataclasses import replace
from fractions import Fraction as F
import argparse,csv,hashlib,json,math,platform
import numpy as np
import z3
from activity import ActivityBox,Core,ActivityProblem,windmill
from solver import CycleAccelerator,strict_decision,capacity_bracket,round_down,full_system,status,q,truth,exact,encode,numeric,BudgetUnknown
from elimination import GridElimination,ClosedUnionNext
from paper_checks import replay

MANUSCRIPT_SHA256='cdd3d14f5f78c6d309b35c7702b8379750fe5df5b5b08ec8e76bf80d3b9336ab'


def serialize(value):
    if isinstance(value,F):return str(value)
    if isinstance(value,z3.ExprRef):return encode(value)
    if isinstance(value,np.ndarray):return value.tolist()
    if isinstance(value,np.generic):return value.item()
    raise TypeError(type(value).__name__)


def unit_problem():
    return ActivityProblem({'u':ActivityBox(F(1,10),F(9,10)),'v':ActivityBox(F(1,100),F(9,10))},(Core('u','v'),))


def main():
    ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('--output',default='outputs');args=ap.parse_args();out=Path(args.output);out.mkdir(parents=True,exist_ok=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)
    dump('fresh_source_identities.json',replay())
    # Algorithm 1: exact forced-root jump, followed by an independent leastness query.
    unit=unit_problem();algorithm=CycleAccelerator(unit,F(1,100),SOLVER_TIMEOUT_MS,MAX_CYCLE_JUMPS);one=algorithm.solve()
    assert one['status']=='feasible' and one['jumps']==1
    alpha=exact((5-z3.Sqrt(13))/10)
    assert truth(one['state']['u']==alpha) and truth(one['state']['v']==alpha-q(F(7,100)))
    one['independent_leastness_check']=algorithm.audit_least(one['state']);one['rounded_half_margin_witness']=round_down(unit,one['state'],F(1,100))
    one['scope']='Exact fixed-margin Algorithm 1 with general QF_NRA algebraic roots. No projected-CAD, positive-infinitesimal or FPT implementation claim.'
    dump('unit_cycle_trace.json',one)
    iterations=[];x=.1;af=(5-math.sqrt(13))/10
    for i in range(40):iterations.append((i,x,af-x));x=(-1+math.sqrt(1+8*((x+2*x*x)/3+.02)))/2
    table('plain_iteration_diagnostic.csv',['iteration','activity','distance_to_exact_fixed_point'],iterations)
    bracket=capacity_bracket(unit,CAPACITY_BRACKET_BITS,SOLVER_TIMEOUT_MS)
    assert bracket['status']=='bracketed' and bracket['lower']<=F(1,48)<bracket['upper']
    dump('capacity_bracket.json',dict(unit=bracket,unit_exact='1/48',unit_optimizer={'u':'1/2','v':'17/48'},universal_ceiling='(3-2*sqrt(2))/8',scope='Capacity bisection is monotone in the margin. Binary search on a species coordinate is not generally valid. The universal ceiling is not attained by a rational ratio.'))
    # A tight policy can choose the wrong cycle root if it drops other constraints.
    policy=ActivityProblem({'x':ActivityBox(F(1,8),F(7,8)),'y':ActivityBox(F(1,16),F(7,8)),'z':ActivityBox(F(41,64),F(41,64))},(Core('x','y'),Core('x','z')))
    pa=CycleAccelerator(policy,F(1,64));pr=pa.solve();assert pr['status']=='feasible';assert all(truth(pr['state'][v]==q(value)) for v,value in {'x':F(3,4),'y':F(41,64),'z':F(41,64)}.items());pa.audit_least(pr['state'])
    pr['omitted_constraint_violation_at_selected_policy_least']=F(1,2);dump('policy_counterexample.json',pr)
    boxes={v:ActivityBox(F(1,100),F(9,10)) for v in 'ABCD'}
    diamond=ActivityProblem(boxes,(Core('A','B'),Core('B','C'),Core('A','D',F(2)),Core('D','C',F(2))))
    s,_=full_system(diamond,F(0));assert status(s)==z3.sat
    ds=strict_decision(diamond);assert ds['status']=='infeasible'
    dj=CycleAccelerator(diamond,F(1,10000)).solve();assert dj['status']=='infeasible'
    dump('zero_margin_diamond.json',dict(weak_zero_feasible=True,strict_decision=ds,positive_margin_trace=dj,identity='F(F(x)-t)-t - (F(F(x)+t)+t) = -t*(x*x+x+3)',scope='Zero-margin feasibility does not imply strict productivity. Repeatedly trying smaller margins would not prove infeasibility.'))
    # Exact finite-grid elimination preserves every gap; it does not solve the continuum FPT problem.
    nxt=ClosedUnionNext(((F(1,10),F(1,5)),(F(3,5),F(4,5))))
    assert nxt.value(F(3,10))==F(3,5) and nxt.value(F(9,10)) is None
    dump('disconnected_next_map.json',dict(admissible_intervals=nxt.intervals,lower_request='3/10',least_admissible=nxt.value(F(3,10)),outgoing_cap='2/5',extension_exists=False,
        scope='Abstract monotone-elimination example: with incoming lower bound .3 and outgoing identity capped at .4, no extension exists. Replacing the two intervals by their hull would invent one. This is not claimed to be a projection of the default chemical triangle.'))
    # Configured all-size assembly, each residual checked per core.
    problem=windmill(MODULES,F(PATH_FACTOR_A),F(PATH_FACTOR_B),F(SHORTCUT_FACTOR_A))
    problem.boxes={v:ActivityBox(*(F(x) for x in (HUB_BOX if v=='A' else MIDDLE_BOX if v.startswith('B') else END_BOX))) for v in problem.boxes}
    margin=F(FIXED_MARGIN);configured=CycleAccelerator(problem,margin,SOLVER_TIMEOUT_MS,MAX_CYCLE_JUMPS);cr=configured.solve()
    if cr['status']=='feasible':
        cr['independent_leastness_check']=configured.audit_least(cr['state'])
        if margin>0:cr['rounded_half_margin_witness']=round_down(problem,cr['state'],margin)
    dump('configured_cycle_closure.json',cr)
    strict=strict_decision(problem,SOLVER_TIMEOUT_MS);dump('configured_strict_decision.json',strict)
    grid=GridElimination(problem,margin,GRID_DENOMINATOR).solve();dump('configured_grid_elimination.json',grid)
    if grid['status']=='grid_feasible' and cr['status']=='feasible':
        assert all(truth(q(grid['state'][v])>=cr['state'][v]) for v in problem.boxes)
    error=F(RATIO_RELATIVE_FACTOR_ERROR)
    if not 0<=error<1:raise ValueError('factor uncertainty must lie in [0,1)')
    robust=ActivityProblem(problem.boxes,tuple(replace(e,ratio_lower=e.a/e.b*(1-error)/(1+error),ratio_upper=e.a/e.b*(1+error)/(1-error)) for e in problem.cores))
    rd=strict_decision(robust,SOLVER_TIMEOUT_MS);rd['scope']='One common state for every independent ratio realization, using exact worst-case band endpoints. This is not one separately chosen state per realization.';dump('configured_robust_decision.json',rd)
    nominal_values=list(map(F,NOMINAL_ACTIVITIES));state={v:nominal_values[0 if v=='A' else 1 if v.startswith('B') else 2] for v in problem.boxes}
    accounts=problem.total_accounts(state);accounts['per_core_strictly_productive']=problem.check_rational(state,strict=True)
    tolerance=[];mins=[]
    for e in problem.cores:
        x,y=state[e.tail],state[e.head];r=e.currents(x,y)
        if min(r['p'],r['q'],r['tail_production'],r['head_production'])<=0:continue
        tolerance.extend([r['tail_production']/(2*r['q']+r['p']),r['head_production']/(r['p']+r['q'])])
        lo,hi=e.uncertain_factor_activity_minima(x,y,error,F(ACTIVITY_ERROR));mins.append(dict(edge=[e.tail,e.head],tail_minimum=lo,head_minimum=hi))
    accounts['fixed_state_relative_factor_radius']=min(tolerance) if len(tolerance)==2*len(problem.cores) else None
    accounts['joint_box_minima']=mins;accounts['joint_box_certified']=len(mins)==len(problem.cores) and all(min(r['tail_minimum'],r['head_minimum'])>0 for r in mins)
    accounts['illustrative_units']=dict(concentration_reference_M=REFERENCE_CONCENTRATION_MOLAR,flux_reference_M_per_min=REFERENCE_FLUX_MOLAR_PER_MINUTE,activity_values_uM=[float(v)*REFERENCE_CONCENTRATION_MOLAR*1e6 for v in nominal_values],joint_worst_residual_uM_per_min=min(float(v) for r in mins for v in [r['tail_minimum'],r['head_minimum']])*REFERENCE_FLUX_MOLAR_PER_MINUTE*1e6 if mins else None)
    dump('configured_operating_accounts.json',accounts)
    # The canonical module is always independently reconstructed, irrespective of edited factors.
    reference=windmill(1);reference_state={'A':F(1,10),'B0':F(13,200),'C0':F(43,1000)};ref=reference.total_accounts(reference_state)
    assert ref['food']==F(5071,40000)
    shortcut=ActivityProblem(reference.boxes,tuple(replace(e,a=F(1),ratio_lower=F(1),ratio_upper=F(1)) for e in reference.cores))
    summed=shortcut.total_accounts(reference_state);summed['strict_per_core_decision']=strict_decision(shortcut)
    assert all(v>0 for v in summed['species_production'].values()) and summed['strict_per_core_decision']['status']=='infeasible'
    summed['scope']='All aggregate species productions are positive at this state, yet the unit shortcut triangle has no common per-core productive state. Summing residuals answers a weaker question.'
    dump('aggregate_production_counterexample.json',summed)
    radius=[];uncertain=[]
    for e in reference.cores:
        x,y=reference_state[e.tail],reference_state[e.head];r=e.currents(x,y)
        radius.extend([r['tail_production']/(2*r['q']+r['p']),r['head_production']/(r['p']+r['q'])]);uncertain.extend(e.uncertain_factor_activity_minima(x,y,F(1,20),F(1,10000)))
    assert min(radius)==F(19,301) and min(uncertain)==F(1175221,2000000000)
    ref.update(factor_radius=min(radius),joint_box_minimum=min(uncertain),normalized_margin=F(209,120000));dump('canonical_module_certificate.json',ref)
    scaling=[]
    for k in [1,2,3,5,10,100]:scaling.append((k,2*k+1,3*k,2 if k==1 else 4,str(k*ref['food']),str(ref['food']),str(F(209,40000*k))))
    table('windmill_resource_scaling.csv',['modules','species','cores','longest_path_edges','food_unscaled','food_after_common_factor_division','minimum_nominal_production_after_division'],scaling)
    plot(out,iterations,ref,cr,grid)
    print('One unit core reaches its irrational least state in one exact cycle jump; full-system audit confirms leastness.',flush=True)
    print(f'Configured {MODULES}-module assembly: fixed margin {cr["status"]}; strict {strict["status"]}; robust {rd["status"]}; finite grid {grid["status"]}.',flush=True)
    print('Canonical factor tolerance 19/301; joint activity/factor residual >= 1175221/2000000000; food 5071/40000 per module.',flush=True)
    print('QF_NRA decisions and exact fixed-margin control flow are implemented; the positive-infinitesimal projected-CAD and FPT backend are not. 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(),z3=z3.get_version_string(),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,iterations,ref,configured,grid):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,3,figsize=(12,4),layout='constrained')
    x=np.linspace(.05,.25,301);T=(-1+np.sqrt(1+8*((x+2*x*x)/3+.02)))/2;alpha=(5-math.sqrt(13))/10
    axs[0].plot(x,T,label='Cycle map');axs[0].plot(x,x,color='gray',ls=':',label='Identity');axs[0].scatter([.1,alpha],[.1,alpha]);axs[0].annotate('forced root',xy=(alpha,alpha),xytext=(.07,.2),arrowprops={'arrowstyle':'->'});axs[0].set(xlabel='Seed activity',ylabel='Cycle value',title='Cycle response and its exact fixed point');axs[0].legend(fontsize=7)
    axs[1].semilogy([r[0] for r in iterations],[r[2] for r in iterations]);axs[1].set(xlabel='Plain iteration',ylabel='Distance to exact fixed point',title='Iterative error relative to the exact fixed\npoint')
    t=np.linspace(0,1/48,300);root=np.sqrt(1-48*t);axs[2].plot(t,(1-root)/2);axs[2].plot(t,(1+root)/2);axs[2].set(xlabel='Normalized margin',ylabel='Cycle fixed point',title='Branches merge at capacity 1/48')
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'finite_cycle_jump.png',dpi=180);fig.savefig(out/'finite_cycle_jump.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    names=['AB','BC','AC'];ru=[float(r['tail_production']) for r in ref['cores']];rv=[float(r['head_production']) for r in ref['cores']];pos=np.arange(3)
    axs[0].bar(pos-.17,ru,.34,label='Tail production');axs[0].bar(pos+.17,rv,.34,label='Head production');axs[0].set(xticks=pos,xticklabels=names,ylabel='Production / reference flux',title='Species production within each triangle\ncore');axs[0].legend(fontsize=8)
    k=np.arange(1,21);axs[1].plot(k,k*float(ref['food']),label='Original factors');axs[1].plot(k,np.full(len(k),float(ref['food'])),label='All factors divided by module count');axs[1].set(xlabel='Modules sharing one hub',ylabel='Aggregate food / reference flux',title='Total food demand versus shared-hub module\ncount');axs[1].legend(fontsize=8)
    for ax in axs:ax.grid(alpha=.2)
    fig.savefig(out/'compatible_module_and_food.png',dpi=180);fig.savefig(out/'compatible_module_and_food.svg');plt.close(fig)


if __name__=='__main__':main()

Run output

One unit core reaches its irrational least state in one exact cycle jump; full-system audit confirms leastness.
Configured 3-module assembly: fixed margin feasible; strict feasible; robust feasible; finite grid grid_feasible.
Canonical factor tolerance 19/301; joint activity/factor residual >= 1175221/2000000000; food 5071/40000 per module.
QF_NRA decisions and exact fixed-margin control flow are implemented; the positive-infinitesimal projected-CAD and FPT backend are not. Lean not rerun.