"""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()