Example code

"""Edit these inputs, then run: python example.py --output outputs.

Default run freshly certifies both support frontiers (~several minutes), searches
the cooperative class, and simulates the actual harvest/partition/refill protocol.
"""
INVENTORY = 32
QUOTA = 4
DEADLINE = 5                       # model time, fixed-volume count-rate convention
RELEASE_BAND = (10, 20)             # tied dissociation and productive release
ALLOCATION_THETA = '1/2'
RATE_BOX_HALF_WIDTH = '1/300'        # five groups: the two release channels stay tied
SUPPORT_TARGET = '999/1000'
FRONTIER_BUDGETS = range(15, 43)     # each lower budget is checked, never interpolated
COOPERATIVE_TARGETS = [('999/1000',1),('1999/2000',1),('1999/2000',8)]
COOPERATIVE_SEARCH_CAP = 800
REFERENCE_MAJORITY = (65,195,121)    # low, high, growth
REFERENCE_MINORITY = (8,64,30)
GATE_FORWARD_FLOOR = '20'
GATE_REVERSE_CEILING = '1/100000'
LINEAGE_CYCLES = 10
RANDOM_SEED = 740023
MANUSCRIPT_SHA256 = 'e95aab6b325f1d7d500ee7f13bf8d287e15f308bf38a561f02e5ebf3d819e32c'

import argparse,csv,hashlib,json,platform
from pathlib import Path
from fractions import Fraction as F
from math import log
import numpy as np
import scipy
from dataclasses import replace
from memory import (SupportReactor,SupportRates,ComplementaryAllocation,IntegerEnvelope,
                    SerialProtocol,Corridor,CooperativeGate,carrier_floor)
from ctmc_certificate import nominal_groups,SCALE
from frontier_proportion import frontier,verify_endpoint_lemma
from proportion_extras import spread,log10_k


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,value): (out/name).write_text(json.dumps(value,indent=2,default=lambda x:str(x) if isinstance(x,F) else 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)
    theta=ComplementaryAllocation(F(ALLOCATION_THETA));target=F(SUPPORT_TARGET)
    def certify(K,groups,allocation=theta):
        r=IntegerEnvelope(K,tuple(groups),DEADLINE,QUOTA,allocation).run()
        r['status']='certified' if F(r['lower'],SCALE)>=target else ('excluded' if F(r['upper'],SCALE)<target else 'undetermined')
        return r
    band=certify(INVENTORY,nominal_groups(*RELEASE_BAND));dump('support_band.json',band)
    print('Configured support band:',band['status'],band['lower'],'/',SCALE,flush=True)
    nominal=[];untied=[]
    for K in FRONTIER_BUDGETS:
        if K<=INVENTORY:
            r=certify(K,SupportRates().groups());nominal.append(r)
        r=certify(K,SupportRates(dissociation=F(20),release=F(10)).groups());untied.append(r)
        print('Fresh budget',K,'untied:',r['status'],flush=True)
    dump('support_frontiers.json',dict(nominal=nominal,untied=untied,
        analytic_exclusion_through=carrier_floor(QUOTA,1,1-target)-1,
        scope='Minimum asserted only with continuous coverage of all smaller integers by analytic or upper certificates. No monotonicity in inventory, deadline or individual rates.'))
    for name,rows in [('nominal',nominal),('untied',untied)]:
        successes=[r['K'] for r in rows if r['status']=='certified']
        if successes:
            k=min(successes);covered=set(range(2,carrier_floor(QUOTA,1,1-target)))|{r['K'] for r in rows if r['status']=='excluded'}
            print(name,'minimum' if set(range(2,k))<=covered else 'certified budget (coverage incomplete)',k,flush=True)
    robustness={
        'bias_039':certify(INVENTORY,nominal_groups(*RELEASE_BAND),ComplementaryAllocation(F(39,100))),
        'bias_038':certify(INVENTORY,nominal_groups(*RELEASE_BAND),ComplementaryAllocation(F(19,50))),
        'tied_rate_box':certify(INVENTORY,nominal_groups(*RELEASE_BAND,eps=F(RATE_BOX_HALF_WIDTH))),
        'wider_rate_box':certify(INVENTORY,nominal_groups(*RELEASE_BAND,eps=F(1,250))),
        'faster_product_release':certify(INVENTORY,SupportRates(dissociation=F(10),release=F(20)).groups())}
    dump('support_robustness.json',robustness)
    diagnostics=[]
    for kd,kr in [(20,10),(10,10),(10,20)]:
        reactor=SupportReactor(INVENTORY,SupportRates(dissociation=F(kd),release=F(kr)))
        states,values=reactor.deadline(DEADLINE,QUOTA,theta)
        assert np.max(np.abs(values[:,:3].sum(axis=1)-1))<1e-9
        for z in [(INVENTORY-1,0,0),(0,1,0)]:
            v=values[states.index(z)]
            diagnostics.append(dict(dissociation=kd,product_release=kr,start=str(z),joint_success=v[0],quota_failure=v[1],allocation_failure_after_quota=v[2],mean_carriers=v[3],mean_food=v[4],mean_product=v[5]))
    table('failure_mechanisms.csv',diagnostics)
    protocol=SerialProtocol(SupportReactor(INVENTORY),DEADLINE,QUOTA,theta)
    lineage=protocol.lineage((INVENTORY-1,0,0),LINEAGE_CYCLES,np.random.default_rng(RANDOM_SEED));dump('designated_lineage.json',lineage)
    lower=F(band['lower'],SCALE)
    dump('material_and_time_bounds.json',dict(support_carrier_floor=carrier_floor(QUOTA,1,1-target),two_type_carrier_floor=carrier_floor(QUOTA,2,1-target),
        lineage_joint_lower=lower**LINEAGE_CYCLES,initial_plus_refills_on_success_at_least=INVENTORY+LINEAGE_CYCLES*(INVENTORY+QUOTA),
        harvest_on_success_at_least=LINEAGE_CYCLES*QUOTA,
        minimum_tied_release_necessary_display=log(1/float(1-target))/(2*DEADLINE),
        scope='Lineage probability bound uses a uniform conditional guarantee and requires no cycle independence. Necessary release threshold is numerical display; material ceiling is exact. Not a whole binary-tree success guarantee.'))
    ref=CooperativeGate(Corridor(*REFERENCE_MAJORITY),Corridor(*REFERENCE_MINORITY),QUOTA,F(GATE_FORWARD_FLOOR),F(GATE_REVERSE_CEILING))
    def report_gate(g):
        X,Y=g.majority,g.minority;k,b=g.constants()
        return dict(core=g.core,initial_supply=g.core+1,molecularity=g.molecularity,majority=X.__dict__,minority=Y.__dict__,forward_constant=k,reverse_constant=b,
            partition=X.uniform_return()*Y.uniform_return(),joint_lower=g.joint_lower(),ten_cycle_lower=g.joint_lower()**10,
            worst_clock_and_partition_need_not_coincide=True,
            vertices=[dict(x=x,y=y,forward=g.pair_rates(x,y)[0],reverse=g.pair_rates(x,y)[1],joint_display=g.actual_joint(x,y)) for x in [X.low,X.high] for y in [Y.low,Y.high]])
    reference=report_gate(ref);reference['bias']=[dict(theta=str(t),partition=ref.majority.uniform_return(t)*ref.minority.uniform_return(t),joint_lower=ref.joint_lower(t)) for t in [F(1,2),F(49,100),F(12,25)]]
    reference['scope']='The 414-core reference and its bias tolerance are distinct from the optimized 230- or 253-core designs. No tolerance is transferred between them.'
    dump('cooperative_reference.json',reference)
    optimal=[]
    for rho,ly in COOPERATIVE_TARGETS:
        w=frontier(F(rho),q=QUOTA,ly_min=ly,cap=COOPERATIVE_SEARCH_CAP);w.pop('seconds',None)
        g=CooperativeGate(Corridor(w['lx'],w['ux'],w['gx']),Corridor(w['ly'],w['uy'],w['gy']),QUOTA,F(GATE_FORWARD_FLOOR),F(GATE_REVERSE_CEILING))
        w['realized_gate']=report_gate(g)
        w['configured_clock_meets_target']=g.joint_lower()>F(rho)
        optimal.append(w);print('Cooperative class',rho,'minority floor',ly,'minimum core',w['N'],flush=True)
    dump('cooperative_frontiers.json',optimal)
    scaling=[]
    for power in range(2,8):
        w=frontier(1-F(1,10**power),q=QUOTA,cap=COOPERATIVE_SEARCH_CAP);w.pop('seconds',None)
        scaling.append(dict(failure_target=str(F(1,10**power)),core=w['N'],molecularity=w['lx']+w['ly']+w['gx']+w['gy']+QUOTA+1,core_over_log_failure=w['N']/log(10**power)))
    table('cooperative_scaling.csv',scaling)
    w=dict(lx=ref.majority.low,ux=ref.majority.high,gx=ref.majority.growth,ly=ref.minority.low,uy=ref.minority.high,gy=ref.minority.growth,q=QUOTA)
    dump('cooperative_rate_spread.json',spread(w,ref.core,log10_k(w['lx'],w['ly'],ref.food_consumed,float(ref.forward_floor))))
    cases=[(l,u,g) for l in range(1,9) for u in range(l,25) for g in range(l,u+1)]
    dump('endpoint_reduction_replay.json',dict(cases=len(cases),violations=verify_endpoint_lemma(cases)))
    plot(out,nominal,untied,diagnostics,scaling)
    print('Joint quota/allocation certificates regenerated; full finite generator used. Floating diagnostics and one sampled lineage are not guarantees.',flush=True)
    print('Support rate box has five independent groups with tied exits. Cooperative gates permit extreme reaction order/rates; counts are not energy costs. 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(),numpy=np.__version__,scipy=scipy.__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,nominal,untied,diagnostics,scaling):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(11,4),layout='constrained')
    for rows,label in [(nominal,'Balanced exits'),(untied,'Dissociation favored')]:
        axs[0].semilogy([r['K'] for r in rows],[1-r['upper']/SCALE for r in rows],'.-',label=label+' (lower failure)')
        axs[0].semilogy([r['K'] for r in rows],[1-r['lower']/SCALE for r in rows],':',alpha=.7)
    axs[0].axhline(1-float(F(SUPPORT_TARGET)),color='black',ls='--');axs[0].set(xlabel='Working inventory K',ylabel='Worst-start failure probability',title='Certified failure bounds versus working\ninventory');axs[0].legend(fontsize=7)
    x=np.arange(len(diagnostics));axs[1].bar(x,[r['quota_failure'] for r in diagnostics],label='Quota missed');axs[1].bar(x,[r['allocation_failure_after_quota'] for r in diagnostics],bottom=[r['quota_failure'] for r in diagnostics],label='Quota met, daughter lost');axs[1].set(xticks=x,xticklabels=['20/10\nresident','20/10\nseed','10/10\nresident','10/10\nseed','10/20\nresident','10/20\nseed'],ylabel='Failure probability',title='Quota and daughter-allocation failures by\nreaction rates');axs[1].legend(fontsize=8)
    for a in axs:a.grid(alpha=.2)
    fig.savefig(out/'support_frontier_and_failure.png',dpi=180);fig.savefig(out/'support_frontier_and_failure.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4),layout='constrained')
    axs[0].bar(['Support\nbenchmark','Cooperative\nclass'],[32,231],color=['#537ba7','#bf8646']);axs[0].set(ylabel='Initial supplied units, including token',title='Quota 4, reliability 0.999');axs[0].text(.5,245,'Different architectural scopes',ha='center',fontsize=9);axs[0].set_ylim(0,280)
    x=[log(1/float(F(r['failure_target']))) for r in scaling];axs[1].plot(x,[r['core'] for r in scaling],'o-',label='Exact class frontier');axs[1].plot(x,[33.3*v for v in x],':',label='33.3 log(1 / failure): empirical');axs[1].set(xlabel='log(1 / failure target)',ylabel='Minimum core inventory',title='Minimum core inventory versus reliability\ntarget');axs[1].legend(fontsize=8)
    for a in axs:a.grid(alpha=.2)
    fig.savefig(out/'cooperative_material_cost.png',dpi=180);fig.savefig(out/'cooperative_material_cost.svg');plt.close(fig)


if __name__=='__main__': main()