"""Exact boundary reduction and rational operating states for junction paths.

python example.py --output outputs
python example.py --input assembly.json --output custom_outputs
Compatibility is static; it does not imply a sustained operating trajectory.
"""
# EDITABLE INPUTS: exact rational model conductances and dimensionless activities.
LOWER_ACTIVITY = '9/10'
PATH_FACTORS = (('9','1'),('9','1'))
SHORTCUT_FACTORS = ('2','1')
REFERENCE_ACTIVITIES = {'A':'99/100','B':'593/600','C':'3551/3600'}
RELATIVE_FACTOR_RADIUS = '1/25'
ACTIVITY_RADIUS = '1/1000000'
FOOD_INTERVAL = ('1','1')
MAX_JUNCTIONS = 6
MAX_PATH_LENGTH = 8
QUERY_TIMEOUT_MS = 3000
MAX_QUERIES = 400

from dataclasses import dataclass
from fractions import Fraction as F
from itertools import product
import argparse
import csv
import hashlib
import json
from pathlib import Path
import platform
import numpy as np
import z3

MANUSCRIPT_SHA256='c71bb73dd23e95817ca599f6a14b3fcb7ca4ec611d6549e411292a3aa23b80c0'


def rational(value):
    if isinstance(value,F):return value
    if type(value) is int or isinstance(value,str):return F(value)
    raise ValueError('Use integer or exact rational strings, not floating-point inputs.')


class Unresolved(RuntimeError):pass


@dataclass(frozen=True)
class Core:
    source:str
    target:str
    a:F
    b:F

    def __post_init__(self):
        object.__setattr__(self,'a',rational(self.a));object.__setattr__(self,'b',rational(self.b))
        if min(self.a,self.b)<=0:raise ValueError('Positive barrier factors required.')
    def response(self,x,c):return (self.a*x+c*self.b*x*x)/(self.a+c*self.b)
    def lower(self,x):return self.response(x,2)
    def upper(self,x):return self.response(x,1)
    def currents(self,x,y,food=F(1)):
        p=self.a*(x-y);q=self.b*(food*y-x*x)
        return p,q,2*q-p,p-q


@dataclass(frozen=True)
class PolynomialResponse:
    """Nonnegative integer coefficients / positive denominator, un-reduced."""
    coefficients:tuple[int,...]=(0,1)
    denominator:int=1

    def append(self,core,c):
        if c not in (1,2):raise ValueError('Response selector must be 1 or 2.')
        a,b=core.a,core.b;alpha=a.numerator*b.denominator;beta=b.numerator*a.denominator
        n=self.coefficients;Q=self.denominator;values=[0]*(2*len(n)-1)
        for i,u in enumerate(n):
            for j,v in enumerate(n):values[i+j]+=c*beta*u*v
        for i,u in enumerate(n):values[i]+=alpha*Q*u
        denominator=(alpha+c*beta)*Q*Q
        if min(values)<0 or sum(values)!=denominator:raise ArithmeticError('Positive representation invariant failed.')
        return PolynomialResponse(tuple(values),denominator)
    def evaluate(self,x):
        value=0
        for c in reversed(self.coefficients):value=value*x+c
        return value/self.denominator
    def z3(self,x):
        value=z3.RealVal(0)
        for c in reversed(self.coefficients):value=value*x+c
        return value/self.denominator


@dataclass(frozen=True)
class DirectedPath:
    cores:tuple[Core,...]

    @property
    def vertices(self):return (self.cores[0].source,)+tuple(e.target for e in self.cores)
    def response(self,x,c):
        for core in self.cores:x=core.response(x,c)
        return x
    def polynomial(self,c):
        value=PolynomialResponse()
        for core in self.cores:value=value.append(core,c)
        return value
    def reconstruct(self,x,y,ell):
        """Keep the final endpoint exact; only the private prefix is interpolated."""
        x,y,ell=map(rational,(x,y,ell));lo,hi=self.response(x,2),self.response(x,1)
        if not 0<ell<=y<x<1 or not lo<y<hi:raise ValueError('Strict response band and common positive box required.')
        if len(self.cores)==1:return (x,y),{'iterations':0,'proved_budget':0}
        r=len(self.cores);C=2**r-1;Cp=2**(r-1)-1;last=self.cores[-1]
        delta=min(y-lo,hi-y);gap=last.a*last.b*ell*(1-x)/((last.a+last.b)*(last.a+2*last.b));radius=delta*gap/(4*C*Cp)
        width=F(1);budget=1
        while width>=radius:width/=2;budget+=1
        left,right=F(0),F(1)
        for iteration in range(1,budget+1):
            lam=(left+right)/2;values=[x]
            for edge in self.cores[:-1]:values.append((1-lam)*edge.lower(values[-1])+lam*edge.upper(values[-1]))
            t=values[-1]
            if last.lower(t)<y<last.upper(t):
                values.append(y)
                if not all(ell<=v<=1 for v in values) or not all(min(edge.currents(u,v)[2:])>0 for edge,u,v in zip(self.cores,values,values[1:])):raise ArithmeticError('Private reconstruction failed exact check.')
                return tuple(values),{'iterations':iteration,'proved_budget':budget,'interpolation_parameter':str(lam),'certified_parameter_radius':str(radius)}
            end=(1-lam)*last.lower(t)+lam*last.upper(t)
            if end<y:left=lam
            else:right=lam
        raise ArithmeticError('Proved rational reconstruction budget exhausted.')


class Assembly:
    """Validates a supplied decomposition; never invents private/junction copies."""
    def __init__(self,data):
        self.species=tuple(data['species']);self.junctions=tuple(data['junctions']);self.ell=rational(data['ell'])
        if any(not isinstance(s,str) or not s for s in self.species) or len(set(self.species))!=len(self.species) or len(set(self.junctions))!=len(self.junctions):raise ValueError('Distinct named species/junctions required.')
        if not set(self.junctions)<=set(self.species) or not 0<self.ell<1:raise ValueError('Invalid junctions or lower activity.')
        raw=data.get('boxes',{})
        if set(raw)-set(self.junctions):raise ValueError('Independent private boxes are outside the two-response reduction; promote those species to junctions.')
        self.boxes={s:tuple(map(rational,raw.get(s,(self.ell,1)))) for s in self.junctions}
        if any(len(b)!=2 or not self.ell<=b[0]<=b[1]<=1 for b in self.boxes.values()):raise ValueError('Invalid closed junction box.')
        self.cores=tuple(Core(e['u'],e['v'],e['a'],e['b']) for e in data['edges']);pairs=set()
        for e in self.cores:
            pair=frozenset((e.source,e.target))
            if e.source not in self.species or e.target not in self.species or e.source==e.target or pair in pairs:raise ValueError('Graph must have no loops, duplicates or antiparallel edges.')
            pairs.add(pair)
        used=[];private=set();paths=[]
        for ids in data['paths']:
            if not ids or any(type(i) is not int or not 0<=i<len(self.cores) for i in ids):raise ValueError('Invalid path indices.')
            path=DirectedPath(tuple(self.cores[i] for i in ids));vs=path.vertices
            if any(a.target!=b.source for a,b in zip(path.cores,path.cores[1:])):raise ValueError('Path orientation/discontinuity.')
            if vs[0] not in self.junctions or vs[-1] not in self.junctions or len(set(vs))!=len(vs):raise ValueError('Path endpoints must be junctions, with no repeated vertices.')
            for s in vs[1:-1]:
                if s in self.junctions or s in private:raise ValueError('Interior species must be private to exactly one path.')
                private.add(s)
            used.extend(ids);paths.append(path)
        if sorted(used)!=list(range(len(self.cores))):raise ValueError('Paths must partition the selected edges exactly.')
        self.paths=tuple(paths)

    def check(self,activities):
        if set(activities)!=set(self.species):raise ValueError('Exactly one activity per global species required.')
        x={s:rational(v) for s,v in activities.items()}
        if any(not self.ell<=v<=1 for v in x.values()) or any(not lo<=x[s]<=hi for s,(lo,hi) in self.boxes.items()):raise ValueError('Activity violates a box.')
        currents=[e.currents(x[e.source],x[e.target]) for e in self.cores]
        if any(min(row[2:])<=0 for row in currents):raise ValueError('Nonpositive literal production residual.')
        return currents


def zq(v):
    v=rational(v);return z3.RealVal(f'{v.numerator}/{v.denominator}')


class BoundarySolver:
    """Z3 is an exact NRA backend, not the paper's complexity-bound routine.

    Every SAT witness is reconstructed and checked with rational arithmetic.
    Timeout/query limits return UNKNOWN, never INCOMPATIBLE.
    """
    def __init__(self,max_junctions=MAX_JUNCTIONS,max_length=MAX_PATH_LENGTH,timeout_ms=QUERY_TIMEOUT_MS,max_queries=MAX_QUERIES):
        self.max_junctions,self.max_length,self.timeout_ms,self.max_queries=max_junctions,max_length,timeout_ms,max_queries
        if timeout_ms<=0 or max_queries<0:raise ValueError('Positive per-query timeout and nonnegative query budget required.')
    def solve(self,assembly):
        self.queries=0
        if len(assembly.junctions)>self.max_junctions or any(len(p.cores)>self.max_length for p in assembly.paths):return {'status':'OUTSIDE_IMPLEMENTATION_LIMITS'}
        if not assembly.cores:
            values={s:(sum(assembly.boxes[s])/2 if s in assembly.boxes else assembly.ell) for s in assembly.species};assembly.check(values)
            return {'status':'SAT','activities':{s:str(v) for s,v in values.items()},'queries':0,'reconstruction':[]}
        solver=z3.SolverFor('QF_NRA');solver.set(timeout=self.timeout_ms);variables={s:z3.Real(f'junction_{i}') for i,s in enumerate(assembly.junctions)}
        for s,(lo,hi) in assembly.boxes.items():solver.add(variables[s]>=zq(lo),variables[s]<=zq(hi))
        residuals=[]
        for p in assembly.paths:
            x,y=variables[p.vertices[0]],variables[p.vertices[-1]]
            residuals.extend((y-p.polynomial(2).z3(x),p.polynomial(1).z3(x)-y))
        def query(constraints):
            if self.queries>=self.max_queries:raise Unresolved('query budget exhausted')
            self.queries+=1;solver.push();solver.add(*constraints);answer=solver.check();reason=solver.reason_unknown();solver.pop()
            if answer==z3.unknown:raise Unresolved(reason)
            return answer==z3.sat
        try:
            if not query([f>0 for f in residuals]):return {'status':'UNSAT','queries':self.queries,'scope':'Exact nonlinear real-arithmetic backend verdict for the validated reduction; no independent UNSAT proof export.'}
            slack=F(1)
            while not query([f>=zq(slack) for f in residuals]):slack/=2
            solver.add(*[f>=zq(slack) for f in residuals]);D=max(2**len(p.cores) for p in assembly.paths);eta=slack/(4*(D+1))
            brackets=dict(assembly.boxes)
            for s in assembly.junctions:
                lo,hi=brackets[s]
                while hi-lo>eta:
                    mid=(lo+hi)/2;constraints=[c for t,(l,u) in brackets.items() for c in (variables[t]>=zq(l),variables[t]<=zq(u))]
                    if query(constraints+[variables[s]<=zq(mid)]):hi=mid
                    else:lo=mid
                    brackets[s]=(lo,hi)
            values={s:assembly.ell for s in assembly.species}
            for s,(lo,hi) in brackets.items():values[s]=(lo+hi)/2
            reconstruction=[]
            for p in assembly.paths:
                local,stats=p.reconstruct(values[p.vertices[0]],values[p.vertices[-1]],assembly.ell)
                for s,v in zip(p.vertices,local):
                    if s in assembly.junctions and values[s]!=v:raise ArithmeticError('Conflicting junction activity.')
                    values[s]=v
                reconstruction.append(stats)
            rows=assembly.check(values)
            return {'status':'SAT','activities':{s:str(v) for s,v in values.items()},'currents_and_residuals':[list(map(str,row)) for row in rows],
                'queries':self.queries,'normalized_response_slack':str(slack),'brackets':{s:list(map(str,b)) for s,b in brackets.items()},'reconstruction':reconstruction,
                'scope':'Exact rational literal-current witness; Z3 QF_NRA used for decisions. No claim of the manuscript backend complexity or rerun Lean verification.'}
        except Unresolved as error:return {'status':'UNKNOWN','queries':self.queries,'reason':str(error)}


class OperatingCertificate:
    def __init__(self,assembly,activities):
        self.assembly=assembly;self.activities={s:rational(v) for s,v in activities.items()};self.rows=assembly.check(self.activities)
    def relative_radius(self):
        if not self.rows:raise ValueError('Nonempty family required for a finite tolerance radius.')
        radii=[min(ru/(2*q+p),rv/(p+q)) for p,q,ru,rv in self.rows]
        return min(radii),radii
    def relative_margin(self,rho):
        rho=rational(rho)
        if not 0<=rho<1:raise ValueError('Relative radius must lie in [0,1).')
        return min(min(ru-rho*(2*q+p),rv-rho*(p+q)) for p,q,ru,rv in self.rows)
    def joint_rectangle(self,activities,factors,food=(F(1),F(1))):
        """Eight checks per core; all bounds are simultaneous and closed."""
        def box(b):
            b=tuple(map(rational,b))
            if len(b)!=2 or not 0<b[0]<=b[1]:raise ValueError('Positive ordered interval required.')
            return b
        activity={s:box(activities[s]) for s in self.assembly.species};food=box(food)
        if len(factors)!=len(self.rows):raise ValueError('One factor rectangle per edge required.')
        checks=[]
        for i,(edge,pair) in enumerate(zip(self.assembly.cores,factors)):
            a_box,b_box=map(box,pair)
            for a,b in product(a_box,b_box):
                for side in (0,1):
                    x=activity[edge.source][1-side];y=activity[edge.target][side];f=food[side]
                    core=Core(edge.source,edge.target,a,b);residual=core.currents(x,y,f)[2+side]
                    checks.append({'edge':i,'side':'source' if side==0 else 'target','a':str(a),'b':str(b),'x':str(x),'y':str(y),'food':str(f),'residual':str(residual)})
        minimum=min(F(row['residual']) for row in checks)
        return {'productive_throughout':minimum>0,'minimum_residual':str(minimum),'corners':checks}
    def food_window(self):
        intervals=[]
        for edge,row in zip(self.assembly.cores,self.rows):
            x,y=self.activities[edge.source],self.activities[edge.target];p=row[0]
            intervals.append(((x*x+p/(2*edge.b))/y,(x*x+p/edge.b)/y))
        return max(b[0] for b in intervals),min(b[1] for b in intervals)
    def instantaneous_budget(self,degradation=None):
        d={s:F(0) for s in self.assembly.species} if degradation is None else {s:rational(degradation[s]) for s in self.assembly.species}
        if min(d.values())<0:raise ValueError('Nonnegative degradation required.')
        production={s:F(0) for s in self.assembly.species}
        for edge,row in zip(self.assembly.cores,self.rows):production[edge.source]+=row[2];production[edge.target]+=row[3]
        return {'production':{s:str(v) for s,v in production.items()},'food_consumption':str(sum(row[1] for row in self.rows)),
            'strict_dilution_upper':str(min(production[s]/self.activities[s]-d[s] for s in production)),
            'scope':'Only the instantaneous field P_s-(D+d_s)*x_s at this prepared state. No persistence, invariant-region or steady-state conclusion.'}


class DeficitConstruction:
    """Linear sufficient construction; failure does not decide compatibility."""
    @staticmethod
    def from_deficits(assembly,deficits,epsilon=None):
        d={s:rational(deficits[s]) for s in assembly.species}
        if min(d.values())<=0:raise ValueError('Positive deficits required.')
        strict_upper=[]
        for e in assembly.cores:
            U=(e.a+4*e.b)*d[e.source]-(e.a+2*e.b)*d[e.target];V=(e.a+e.b)*d[e.target]-(e.a+2*e.b)*d[e.source]
            if min(U,V)<=0:return {'status':'NO_SUFFICIENT_CERTIFICATE','reason':'strict linear deficits fail'}
            strict_upper.append(U/(2*e.b*d[e.source]**2))
        lo=F(0);hi=(1-assembly.ell)/max(d.values())
        for s,(lower,upper) in assembly.boxes.items():lo=max(lo,(1-upper)/d[s]);hi=min(hi,(1-lower)/d[s])
        bound=min(strict_upper,default=hi+1);cap=min(hi,bound)
        if epsilon is None:
            epsilon=(lo+cap)/2 if lo<cap else lo
        else:epsilon=rational(epsilon)
        if not 0<epsilon<bound or not lo<=epsilon<=hi:return {'status':'NO_SUFFICIENT_CERTIFICATE','reason':'epsilon interval empty or supplied epsilon outside it'}
        values={s:1-epsilon*d[s] for s in assembly.species};rows=assembly.check(values)
        return {'status':'SAT','activities':{s:str(v) for s,v in values.items()},'deficits':{s:str(v) for s,v in d.items()},'epsilon':str(epsilon),'minimum_residual':str(min((min(r[2:]) for r in rows),default=F(0)))}
    @staticmethod
    def linear_program(assembly,timeout_ms=QUERY_TIMEOUT_MS):
        solver=z3.Optimize();solver.set(timeout=timeout_ms);d={s:z3.Real(f'deficit_{i}') for i,s in enumerate(assembly.species)};t=z3.Real('margin')
        solver.add(sum(d.values())==1,*[v>=t for v in d.values()])
        for e in assembly.cores:
            alpha=(e.a+2*e.b)/(e.a+e.b);beta=(e.a+4*e.b)/(e.a+2*e.b)
            solver.add(d[e.target]-zq(alpha)*d[e.source]>=t,zq(beta)*d[e.source]-d[e.target]>=t)
        solver.maximize(t);answer=solver.check()
        if answer==z3.unknown:return {'status':'UNKNOWN','reason':solver.reason_unknown()}
        if answer!=z3.sat:return {'status':'NO_SUFFICIENT_CERTIFICATE'}
        model=solver.model()
        def fraction(value):
            value=model.eval(value)
            if not z3.is_rational_value(value):raise Unresolved('Nonrational linear-program output.')
            return F(value.numerator_as_long(),value.denominator_as_long())
        margin=fraction(t)
        if margin<=0:return {'status':'NO_SUFFICIENT_CERTIFICATE','linear_margin':str(margin)}
        result=DeficitConstruction.from_deficits(assembly,{s:fraction(v) for s,v in d.items()});result['linear_margin']=str(margin);return result


def shortcut_data(private=True,unit=False,singleton=False):
    factors=(('1','1'),)*3 if unit else (*PATH_FACTORS,SHORTCUT_FACTORS)
    return {'species':['A','B','C'],'junctions':['A','C'] if private else ['A','B','C'],'ell':LOWER_ACTIVITY,
        'edges':[dict(u=u,v=v,a=a,b=b) for (u,v),(a,b) in zip([('A','B'),('B','C'),('A','C')],factors)],
        'paths':[[0,1],[2]] if private else [[0],[1],[2]],
        'boxes':{s:[REFERENCE_ACTIVITIES[s]]*2 for s in (['A','C'] if private else ['A','B','C'])} if singleton else {}}


def nongraded_data():
    p4=['A','u1','u2','u3','Z'];p5=['A','v1','v2','v3','v4','Z'];edges=[];paths=[]
    for vertices in (p4,p5):
        indices=[]
        for u,v in zip(vertices,vertices[1:]):indices.append(len(edges));edges.append(dict(u=u,v=v,a='1',b='1'))
        paths.append(indices)
    deficits=dict(zip(p4,map(F,['1','83/50','69/25','459/100','191/25'])));deficits.update(zip(p5,map(F,['1','1501/1000','1127/500','3383/1000','1269/250','191/25'])))
    return {'species':list(deficits),'junctions':['A','Z'],'ell':'9/10','edges':edges,'paths':paths},deficits


def main():
    parser=argparse.ArgumentParser(description=__doc__);parser.add_argument('--output',type=Path,default=Path('outputs'));parser.add_argument('--input',type=Path);args=parser.parse_args();out=args.output;out.mkdir(parents=True,exist_ok=True)
    if args.input:
        data=json.loads(args.input.read_text());assembly=Assembly(data);result=BoundarySolver().solve(assembly)
        (out/'custom_result.json').write_text(json.dumps(result,indent=2)+'\n');print(json.dumps(result,indent=2));return
    fixtures={'weighted_private':shortcut_data(),'weighted_junctions':shortcut_data(False),'unit_incompatible':shortcut_data(unit=True),'singleton_junctions':shortcut_data(singleton=True)}
    fixtures['cycle']={'species':['A','B','C'],'junctions':['A','B','C'],'ell':'1/10','edges':[dict(u=u,v=v,a=1,b=1) for u,v in [('A','B'),('B','C'),('C','A')]],'paths':[[0],[1],[2]]}
    decisions={name:BoundarySolver().solve(Assembly(data)) for name,data in fixtures.items()}
    (out/'fixture_inputs.json').write_text(json.dumps(fixtures,indent=2)+'\n')
    a=Assembly(shortcut_data());state={s:F(v) for s,v in REFERENCE_ACTIVITIES.items()};certificate=OperatingCertificate(a,state);radius,radii=certificate.relative_radius();rho=F(RELATIVE_FACTOR_RADIUS);delta=F(ACTIVITY_RADIUS)
    boxes={s:(v-delta,v+delta) for s,v in state.items()};factors=[((e.a*(1-rho),e.a*(1+rho)),(e.b*(1-rho),e.b*(1+rho))) for e in a.cores]
    joint=certificate.joint_rectangle(boxes,factors,tuple(map(F,FOOD_INTERVAL)));food=certificate.food_window();ng,deficits=nongraded_data();nongraded=Assembly(ng)
    fixed_deficits=DeficitConstruction.from_deficits(nongraded,deficits,F(1,10000));lp=DeficitConstruction.linear_program(nongraded)
    (out/'nongraded_input.json').write_text(json.dumps(ng,indent=2)+'\n')
    result={'decisions':decisions,'reference_state':REFERENCE_ACTIVITIES,'relative_radius':str(radius),'edge_relative_radii':list(map(str,radii)),
        'margin_at_configured_relative_radius':str(certificate.relative_margin(rho)),'joint_rectangle':joint,'food_window_open':list(map(str,food)),
        'instantaneous_loss_budget':certificate.instantaneous_budget(),'nongraded_four_five_explicit':fixed_deficits,'nongraded_four_five_linear_program':lp,
        'scope':'Static common-activity compatibility and fixed-state tolerances. No claim of global dynamics, automatic decomposition, or Z3 complexity bound; Lean not rerun.'}
    (out/'results.json').write_text(json.dumps(result,indent=2)+'\n')
    def write(name,header,rows):
        with (out/name).open('w',newline='') as f:w=csv.writer(f);w.writerow(header);w.writerows(rows)
    write('exact_currents.csv',['edge','p','q','source_production','target_production','relative_radius'],[[e.source+'->'+e.target,*map(str,row),str(r)] for e,row,r in zip(a.cores,certificate.rows,radii)])
    path,shortcut=a.paths;xs=np.linspace(.97,1,251);bands=[[x,float(path.response(x,2)),float(path.response(x,1)),float(shortcut.response(x,2)),float(shortcut.response(x,1))] for x in xs]
    write('response_bands.csv',['xA','path_lower','path_upper','shortcut_lower','shortcut_upper'],bands)
    rhos=[F(i,10000) for i in range(501)];tolerances=[[float(r),*[float(min(ru-r*(2*q+p),rv-r*(p+q))) for p,q,ru,rv in certificate.rows]] for r in rhos]
    write('factor_tolerances.csv',['rho','AB_minimum','BC_minimum','AC_minimum'],tolerances)
    foods=[F(998,1000)+F(i,100000) for i in range(501)];foodrows=[[float(f),*[float(min(e.currents(state[e.source],state[e.target],f)[2:])) for e in a.cores]] for f in foods]
    write('food_response.csv',['food','AB_minimum','BC_minimum','AC_minimum'],foodrows)
    representations=[]
    for h in (1,2,4,8):
        p=DirectedPath(tuple(Core(str(i),str(i+1),F(1),F(1)) for i in range(h)));response=p.polynomial(1)
        representations.append({'length':h,'expanded_degree':len(response.coefficients)-1,'denominator_bits':response.denominator.bit_length(),'coefficient_sum_equals_denominator':sum(response.coefficients)==response.denominator})
    (out/'representation_sizes.json').write_text(json.dumps(representations,indent=2)+'\n')
    lines=[f'Weighted private-path decision: {decisions["weighted_private"]["status"]}; unit-factor shortcut: {decisions["unit_incompatible"]["status"]}; directed cycle: {decisions["cycle"]["status"]}.',
        f'Fixed-state sharp relative radius: {radius} = {float(radius):.8g}; strict inequality required.',
        f'Joint preparation/factor/food corner minimum: {joint["minimum_residual"]}.',
        f'Fixed-state food window: ({food[0]}, {food[1]}).',
        f'Nongraded four/five construction minimum production: {fixed_deficits.get("minimum_residual")}; exact LP candidate: {lp["status"]}.',
        'All SAT states are checked against literal currents and boxes. UNKNOWN does not mean incompatible.']
    (out/'console.txt').write_text('\n'.join(lines)+'\n');print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
    for ax,domain in zip(axs,(np.linspace(.97,1,251),np.linspace(.989,.991,251))):
        for p,label,color in [(path,'Two-edge path','#175477'),(shortcut,'Direct shortcut','#bd5a24')]:
            lower=[float(p.response(x,2)) for x in domain];upper=[float(p.response(x,1)) for x in domain]
            ax.fill_between(domain,lower,upper,alpha=.25,color=color,label=label);ax.plot(domain,lower,color=color,lw=.8);ax.plot(domain,upper,color=color,lw=.8)
        ax.plot(float(state['A']),float(state['C']),'ko',ms=4,label='Exact prepared endpoints');ax.set(xlabel='Activity of A',ylabel='Activity of C');ax.set_xticks(np.linspace(domain[0],domain[-1],5));ax.grid(alpha=.2)
    axs[0].set_title('Compatible path-response bands');axs[1].set_title('Shared endpoint window, enlarged');axs[0].legend(fontsize=8)
    fig.savefig(out/'response_bands.png',dpi=180);fig.savefig(out/'response_bands.svg');plt.close(fig)
    fig,axs=plt.subplots(1,2,figsize=(10,4.2),layout='constrained')
    rr=np.array(tolerances);ff=np.array(foodrows)
    for i,label in enumerate(('AB','BC','AC'),1):axs[0].plot(100*rr[:,0],rr[:,i],label=label);axs[1].plot(ff[:,0],ff[:,i],label=label)
    axs[0].set(xlabel='Independent relative factor radius (%)',ylabel='Worst core production (flux units)',title='Fixed prepared state and unit food')
    axs[1].set(xlabel='Food activity at fixed state and factors',ylabel='Minimum core production (flux units)',title='A finite food window for this state')
    for ax in axs:ax.axhline(0,color='gray',lw=.8);ax.legend(fontsize=8);ax.grid(alpha=.2)
    fig.savefig(out/'operating_windows.png',dpi=180);fig.savefig(out/'operating_windows.svg');plt.close(fig)
    digest=lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    (out/'run_metadata.json').write_text(json.dumps({'manuscript_sha256':MANUSCRIPT_SHA256,'source_sha256':digest(Path(__file__)),'python':platform.python_version(),'z3':z3.get_version_string(),
        'output_sha256':{p.name:digest(p) for p in sorted(out.iterdir()) if p.is_file() and p.name!='run_metadata.json'}},indent=2)+'\n')


if __name__=='__main__':main()