Example code

"""Exact, composable activity models for collective core incompatibility.

Activities and kinetic factors are dimensionless; these are instantaneous
production constraints, not a steady-state or time-dependent reactor model.
"""
from __future__ import annotations

from dataclasses import dataclass
from fractions import Fraction as Q
from pathlib import Path
import argparse
import hashlib
import itertools
import json
import platform

# EDITABLE INPUTS: exact decimal strings or rational numbers, never binary floats.
PATH_LENGTH = 3                 # L >= 2; conflict order is L + 1
SHARED_FACTOR = Q(1)
SHARED_A_BOX = (Q('0.1'), Q('0.9'))
SHARED_B_BOX = (Q('0.1'), Q('0.9'))
# gain, B<->C factor, C<->mA factor, closed C activity box
BRANCH_INPUTS = ((3, Q(4), Q(2), (Q('.65'), Q('.71'))),
                 (4, Q('0.5'), Q(1), (Q('.65'), Q('.76'))),
                 (3, Q(1), Q(1), (Q('.65'), Q('.76'))))
MAX_ROOT_PRODUCT_DEGREE = 400   # exceeding this returns not_computed
MANUSCRIPT_SHA256 = '6b1d307cee4aa5f7cb65edbbc36b3071e2a62851bedd27741451321209cba972'


@dataclass(frozen=True)
class Box:
    lower: Q
    upper: Q

    def __post_init__(self):
        if not 0 < self.lower <= self.upper:
            raise ValueError('A positive, nonempty closed box is required.')

    def contains(self, value):
        return self.lower <= value <= self.upper

    def intersect_open(self, lower, upper):
        if not (lower < upper and self.lower < upper and lower < self.upper):
            return None
        return max(self.lower, min(self.upper, (lower + upper) / 2))


@dataclass(frozen=True)
class PairCore:
    source: int
    target: int

    def currents(self, activities, doubled=None):
        x, y = activities[self.source], activities[self.target]
        t = x*x if doubled is None else doubled[self.source]
        return x-y, y-t

    def residuals(self, activities, doubled=None):
        p, q = self.currents(activities, doubled)
        return 2*q-p, p-q


@dataclass(frozen=True)
class PairAssembly:
    species_count: int
    cores: tuple[PairCore, ...]

    def __post_init__(self):
        if self.species_count < 1 or any(not 0 <= v < self.species_count
                for e in self.cores for v in (e.source, e.target)):
            raise ValueError('Invalid species index.')

    def graded_state(self):
        """Return an exact constructive witness, or None (NOT infeasibility)."""
        adjacency = [[] for _ in range(self.species_count)]
        for edge in self.cores:
            adjacency[edge.source].append((edge.target, 1))
            adjacency[edge.target].append((edge.source, -1))
        ranks = {}
        for start in range(self.species_count):
            if start in ranks:
                continue
            ranks[start] = 0
            component, queue = [start], [start]
            while queue:
                u = queue.pop()
                for v, delta in adjacency[u]:
                    if v in ranks:
                        if ranks[v] != ranks[u] + delta:
                            return None
                    else:
                        ranks[v] = ranks[u] + delta
                        component.append(v)
                        queue.append(v)
            shift = min(ranks[v] for v in component)
            for v in component:
                ranks[v] -= shift
        height = max(ranks.values())
        state = tuple(1-Q(1,20)*Q(5,8)**(height-ranks[v])
                      for v in range(self.species_count))
        assert all(min(e.residuals(state)) > 0 for e in self.cores)
        return state


class ShortcutFamily:
    def __init__(self, length):
        if not isinstance(length, int) or length < 2:
            raise ValueError('Path length must be an integer >= 2.')
        self.length = length
        self.assembly = PairAssembly(length+1, tuple(
            [PairCore(i, i+1) for i in range(length)] + [PairCore(0, length)]))

    def report(self):
        L = self.length
        delta = Q(1,1000*(L-1))
        a = (Q(19,20),) + tuple(Q(941,1000)-(i-1)*delta for i in range(1,L+1))
        doubled = (Q(467,500),) + tuple(a[i]-Q(7,4)*delta for i in range(1,L+1))
        deletion = []
        for omitted in range(L+1):
            cores = self.assembly.cores[:omitted]+self.assembly.cores[omitted+1:]
            state = PairAssembly(L+1, cores).graded_state()
            deletion.append({'omitted': omitted, 'state': state,
                             'minimum_production': min(r for e in cores for r in e.residuals(state))})
        return {'order': L+1, 'physical_status': 'incompatible',
                'certificate': 'f(f(x))-g(x)=x*(x-1)*(3*x*x+9*x+2)/24 < 0 on (0,1)',
                'single_deletions': deletion, 'direction_state': a,
                'minimum_forward_current': min(j for e in self.assembly.cores for j in e.currents(a)),
                'independent_doubled_complexes': doubled,
                'minimum_relaxed_production': min(r for e in self.assembly.cores for r in e.residuals(a,doubled))}


@dataclass(frozen=True)
class Branch:
    gain: int
    forward_factor: Q
    return_factor: Q
    box: Box

    def __post_init__(self):
        if not isinstance(self.gain, int) or self.gain < 2 or min(self.forward_factor,self.return_factor) <= 0:
            raise ValueError('Integer gain >= 2 and positive factors required.')

    def responses(self, a, shared_factor):
        m, x, y, b0 = self.gain, self.forward_factor, self.return_factor, shared_factor
        t, ell, u = a**m, self.box.lower, self.box.upper
        lower = (((x+y)*b0*a+m*x*y*t)/(m*x*y+(x+y)*b0),
                 ((x+y)*ell-y*t)/x, (b0*a+m*y*(t-u))/b0)
        upper = (((x+y)*b0*a+x*y*t)/(x*y+(x+y)*b0), (x*u+b0*a)/(x+b0))
        return lower, upper

    def private_state(self, a, b, shared_factor):
        j, t = shared_factor*(a-b), a**self.gain
        lo = max(b-j/self.forward_factor, t+j/(self.gain*self.return_factor))
        hi = (self.forward_factor*b+self.return_factor*t)/(self.forward_factor+self.return_factor)
        return self.box.intersect_open(lo, hi)

    def residuals(self, a, b, c, shared_factor, monomial=None):
        j = shared_factor*(a-b)
        p = self.forward_factor*(b-c)
        q = self.return_factor*(c-(a**self.gain if monomial is None else monomial))
        return self.gain*q-j, j-p, p-q


@dataclass(frozen=True)
class Fan:
    shared_factor: Q
    a_box: Box
    b_box: Box
    branches: tuple[Branch, ...]

    def __post_init__(self):
        if self.shared_factor <= 0 or not self.branches:
            raise ValueError('Positive shared factor and at least one branch required.')

    def selected(self, selected):
        indices = tuple(range(len(self.branches))) if selected is None else tuple(selected)
        if len(set(indices)) != len(indices) or any(i < 0 or i >= len(self.branches) for i in indices):
            raise ValueError('Invalid or duplicate branch index.')
        return indices

    def witness_at(self, a, selected=None):
        indices = self.selected(selected)
        if not self.a_box.contains(a):
            return None
        if indices:
            responses = [self.branches[i].responses(a,self.shared_factor) for i in indices]
            b = self.b_box.intersect_open(max(v for lo,_ in responses for v in lo),
                                         min(v for _,hi in responses for v in hi))
            if b is None:
                return None
        else:
            b = (self.b_box.lower+self.b_box.upper)/2
        cs = tuple(self.branches[i].private_state(a,b,self.shared_factor) if i in indices
                   else (self.branches[i].box.lower+self.branches[i].box.upper)/2
                   for i in range(len(self.branches)))
        if any(c is None for c in cs):
            raise ArithmeticError('Response reconstruction failed.')
        state = (a,b,*cs)
        if not self.audit(state,indices):
            raise ArithmeticError('Reconstructed state failed independent physical audit.')
        return state

    def audit(self, state, selected=None):
        indices = self.selected(selected)
        if len(state) != len(self.branches)+2:
            return False
        a,b,*cs = state
        return (self.a_box.contains(a) and self.b_box.contains(b)
                and all(branch.box.contains(c) for branch,c in zip(self.branches,cs))
                and all(min(self.branches[i].residuals(a,b,cs[i],self.shared_factor)) > 0 for i in indices))

    def decide(self, selected=None, max_degree=MAX_ROOT_PRODUCT_DEGREE):
        """Exact rational sign-cell decision using SymPy real-root isolation.

        Resource guard returns not_computed. Root-isolation implementation is
        trusted, not formally verified. Every positive answer is audited again.
        """
        import sympy as sp
        indices = self.selected(selected)
        for a in (self.a_box.lower,self.a_box.upper):
            witness = self.witness_at(a,indices)
            if witness is not None:
                return {'status':'compatible','state':witness,'method':'box_endpoint'}
        if self.a_box.lower == self.a_box.upper:
            return {'status':'incompatible','method':'singleton_A_box'}
        a = sp.Symbol('a')
        responses = [self.branches[i].responses(a,self.shared_factor) for i in indices]
        lowers = [v for lo,_ in responses for v in lo]
        uppers = [v for _,hi in responses for v in hi]
        expressions = [u-l for l in lowers for u in uppers]
        expressions += [u-self.b_box.lower for u in uppers]
        expressions += [self.b_box.upper-l for l in lowers]
        polys = []
        for expression in expressions:
            p = sp.Poly(expression,a,domain=sp.QQ)
            if p.degree() <= 0:
                if p.eval(0) <= 0:
                    return {'status':'incompatible','method':'nonpositive_constant'}
            else:
                polys.append(p)
        # Include box endpoints so every in-domain sign cell is bounded.
        product = sp.Poly((a-self.a_box.lower)*(a-self.a_box.upper),a,domain=sp.QQ)
        for p in polys:
            product = sp.lcm(product,p)
            if product.degree() > max_degree:
                return {'status':'not_computed','method':'degree_limit','degree':product.degree()}
        roots = product.sqf_part().intervals(eps=sp.Rational(1,10**12))
        samples = []
        for ((_,left),_),((right,_),_) in zip(roots,roots[1:]):
            sample = Q((left+right)/2)
            if self.a_box.lower < sample < self.a_box.upper:
                samples.append(sample)
                witness = self.witness_at(sample,indices)
                if witness is not None:
                    return {'status':'compatible','state':witness,'method':'exact_sign_cells',
                            'degree':product.degree(),'cells_checked':len(samples)}
        return {'status':'incompatible','method':'exact_sign_cells',
                'degree':product.degree(),'cells_checked':len(samples)}


def reference_fan():
    return Fan(SHARED_FACTOR,Box(*SHARED_A_BOX),Box(*SHARED_B_BOX),
               tuple(Branch(m,x,y,Box(*box)) for m,x,y,box in BRANCH_INPUTS))


def robustness_certificate(rho=Q(1,10**6)):
    """Fixed reference theorem, separate from editable fan inputs.

    Certifies a whole cube of independent factor/endpoint perturbations through
    uniform inequalities, not through sampling its corners.
    """
    if not 0 <= rho <= Q(1,10**6):
        return {'status':'not_certified','reason':'outside implemented theorem radius'}
    fan = Fan(Q(1),Box(Q(1,10),Q(9,10)),Box(Q(1,10),Q(9,10)),(
        Branch(3,Q(4),Q(2),Box(Q(13,20),Q(71,100))),
        Branch(4,Q(1,2),Q(1),Box(Q(13,20),Q(19,25))),
        Branch(3,Q(1),Q(1),Box(Q(13,20),Q(19,25)))))
    eta = Q(1,10**5)
    states = (((0,1),(Q(141,160),Q(93,125),Q(71,100)-eta,Q(13,20)+eta,Q(7,10))),
              ((0,2),(Q(169,200),Q(88,125),Q(67,100),Q(7,10),Q(653,1000))),
              ((1,2),(Q(22,25),Q(157,200),Q(7,10),Q(13,20)+eta,Q(18,25))))
    margins = []
    boxes = (fan.a_box,fan.b_box,*(b.box for b in fan.branches))
    for selected,state in states:
        assert fan.audit(state,selected)
        assert min(min(v-box.lower,box.upper-v) for v,box in zip(state,boxes)) >= eta
        margins.append(min(r for i in selected for r in
            fan.branches[i].residuals(state[0],state[1],state[i+2],fan.shared_factor)))
    a, u, ell = Q(7,8),Q(89,125),Q(81,125)
    split_bounds = ((a+4*u+10*a**4-15*ell)/11,(4*u-a-3*a**3)/6)
    assert max(split_bounds) < -Q(3,500)
    relaxed = (Q(9,10)-eta,Q(37,50),Q(701,1000),Q(651,1000),Q(341,500))
    complexes = {3:Q(157,250),4:Q(61,100)}
    relaxed_margin = min(r for i,b in enumerate(fan.branches) for r in
        b.residuals(relaxed[0],relaxed[1],relaxed[i+2],Q(1),complexes[b.gain]))
    assert all((Q(1,10)+rho)**m < t < (Q(9,10)-rho)**m for m,t in complexes.items())
    assert min(min(v-box.lower,box.upper-v) for v,box in zip(relaxed,boxes)) > rho
    assert min(margins)-5*rho > 0 and relaxed_margin-5*rho > 0
    return {'status':'certified','radius':rho,'negative_split_bounds':split_bounds,
            'full_family_production_upper_bound':-Q(3,500)+5*rho,
            'deletion_reference_margins':margins,
            'deletion_uniform_lower_bound':min(margins)-5*rho,
            'independent_complex_reference_margin':relaxed_margin,
            'independent_complex_uniform_lower_bound':relaxed_margin-5*rho}


def rational_json(value):
    if isinstance(value,Q):
        return str(value)
    if isinstance(value,dict):
        return {k:rational_json(v) for k,v in value.items()}
    if isinstance(value,(list,tuple)):
        return [rational_json(v) for v in value]
    return value


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--output',type=Path,default=Path('outputs'))
    args = parser.parse_args()
    args.output.mkdir(parents=True,exist_ok=True)
    shortcut, fan = ShortcutFamily(PATH_LENGTH).report(), reference_fan()
    subsets = [{'selected':selected,**fan.decide(selected)}
               for k in range(len(fan.branches)+1)
               for selected in itertools.combinations(range(len(fan.branches)),k)]
    result = {'shortcut':shortcut,'fan_subfamilies':subsets,
              'fixed_reference_robustness':robustness_certificate()}
    (args.output/'results.json').write_text(json.dumps(rational_json(result),indent=2)+'\n')
    lines = [f'Path with shortcut: {shortcut["order"]} cores; full family incompatible.',
             'All single deletions have exact productive witnesses.',
             f'Minimum direction current: {shortcut["minimum_forward_current"]}',
             f'Minimum independent-complex production: {shortcut["minimum_relaxed_production"]}']
    lines += [f'Fan branches {r["selected"]}: {r["status"]} ({r["method"]})' for r in subsets]
    (args.output/'console.txt').write_text('\n'.join(lines)+'\n')
    print('\n'.join(lines))
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    xs = [i/500 for i in range(501)]
    f = lambda x:(x+x*x)/2
    g = lambda x:(x+2*x*x)/3
    fig,axes = plt.subplots(1,2,figsize=(11,4.4),layout='constrained')
    axes[0].plot(xs,[f(x) for x in xs],label='One-step upper response f(x)')
    axes[0].plot(xs,[g(x) for x in xs],label='Shortcut lower response g(x)')
    axes[0].plot(xs,[f(f(x)) for x in xs],label='Two-step upper response f(f(x))')
    axes[0].set(xlabel='Source activity x',ylabel='Target activity bound',title='Activity bounds for a two-step path and\nshortcut')
    axes[0].legend(fontsize=8)
    aa = [Q(80,100)+Q(i,10000) for i in range(1001)]
    lower,upper = [],[]
    for a in aa:
        rr = [b.responses(a,fan.shared_factor) for b in fan.branches]
        lower.append(float(max(v for lo,_ in rr for v in lo)))
        upper.append(float(min(v for _,hi in rr for v in hi)))
    axes[1].plot([float(a) for a in aa],lower,label='Largest strict lower bound on B')
    axes[1].plot([float(a) for a in aa],upper,label='Smallest strict upper bound on B')
    axes[1].set(xlabel='Shared activity A',ylabel='Shared activity B bound',title='Shared-activity bounds for three branches')
    axes[1].legend(fontsize=8)
    for ax in axes:
        ax.grid(alpha=.2)
    fig.savefig(args.output/'compatibility.png',dpi=180)
    fig.savefig(args.output/'compatibility.svg')
    plt.close(fig)
    digest = lambda p:hashlib.sha256(p.read_bytes()).hexdigest()
    metadata = {'manuscript_sha256':MANUSCRIPT_SHA256,'source_sha256':digest(Path(__file__)),
                'python':platform.python_version(),'output_sha256':{
                    p.name:digest(p) for p in sorted(args.output.iterdir())
                    if p.is_file() and p.name != 'run_metadata.json'}}
    (args.output/'run_metadata.json').write_text(json.dumps(metadata,indent=2)+'\n')


if __name__ == '__main__':
    main()