"""Shared NADPH: retained-state kinetics, exact thresholds, and a usable design model.

Run: python example.py --output outputs. Concentrations uM; time seconds.
The main inputs are immediately below. The exact paper audit stays fixed and is
reported separately from any edited preparation. See README for evidence scope.
"""
from __future__ import annotations
import argparse
import csv
from dataclasses import dataclass, replace, asdict
from fractions import Fraction as F
import hashlib
import json
import math
from pathlib import Path
import platform
from typing import Protocol
import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import brentq

def positive(**values):
    if any(not math.isfinite(float(v)) or v <= 0 for v in values.values()):
        raise ValueError('Finite positive values required: '+', '.join(values))


@dataclass(frozen=True)
class QuadraticResponse:
    """Au²+(B+C/x)u-Gamma=0; j=pu/(vu+w). B may be negative."""
    A: float
    B: float
    C: float
    Gamma: float
    p: float
    v: float
    w: float

    def __post_init__(self):
        positive(A=self.A, C=self.C, Gamma=self.Gamma, p=self.p, w=self.w)
        if not math.isfinite(float(self.B)) or not math.isfinite(float(self.v)) or self.v < 0:
            raise ValueError('B finite; v finite and nonnegative')

    def polynomial(self, x, u):
        return self.A*u*u+(self.B+self.C/x)*u-self.Gamma

    def carrier(self, x):
        if x == 0: return 0.
        positive(x=x)
        beta=float(self.B+self.C/x)
        disc=math.hypot(beta, 2*math.sqrt(float(self.A*self.Gamma)))
        # Avoid cancellation and also avoid overflow in beta+disc.
        return float(self.Gamma)/(beta/2+disc/2) if beta >= 0 else (-beta/2+disc/2)/float(self.A)

    def current_from_carrier(self, u):
        return self.p*u/(self.v*u+self.w)

    def current(self, x):
        return self.current_from_carrier(self.carrier(x))

    def invert(self, quota):
        positive(quota=quota)
        den=self.p-quota*self.v
        if den <= 0: raise ValueError('Quota reaches the kinetic asymptote')
        Y=quota*self.w/den
        delta=self.Gamma-self.A*Y*Y-self.B*Y
        if delta <= 0: raise ValueError('No finite positive NADPH floor')
        return Y, self.C*Y/delta

    def derivative(self, x):
        u=self.carrier(x); beta=self.B+self.C/x
        du=self.C*u/(x*x*(2*self.A*u+beta))
        return self.p*self.w*du/(self.v*u+self.w)**2

    def enclose(self, x, steps=100):
        """Fresh rational bisection, with no floating root used as evidence."""
        if not isinstance(x, F) or x <= 0: raise ValueError('Positive Fraction x required')
        lo,hi=F(0),F(1)
        while self.polynomial(x,hi) < 0: hi *= 2
        for _ in range(steps):
            mid=(lo+hi)/2
            if self.polynomial(x,mid) <= 0: lo=mid
            else: hi=mid
        return lo,hi


@dataclass(frozen=True)
class GlutathioneBranch:
    G: float; z0: float; E: float; a: float; b: float; c: float; k: float

    def __post_init__(self):
        positive(G=self.G,E=self.E,a=self.a,b=self.b,c=self.c,k=self.k)
        if not math.isfinite(float(self.z0)) or self.z0 < 0: raise ValueError('Invalid GSSG baseline')
        self.response  # enforce positive-pool condition at construction

    @property
    def rho(self): return self.a*(1/self.b+1/self.c)
    @property
    def S(self): return self.G-2*self.z0
    @property
    def response(self):
        return QuadraticResponse(1,self.rho-self.S,2*self.E*self.a/self.k,
            self.S*self.rho-self.E*self.a/self.c,self.E*self.a,1,self.rho)

    def reconstruct(self,x):
        positive(x=x); g=self.response.carrier(x); j=self.response.current_from_carrier(g)
        return dict(g=g,z=self.z0+j/(self.k*x),e0=j/self.a,e1=j/(self.b*g),e2=j/(self.c*g),jG=j)

    def redesign(self,x,quota):
        """Both algebraic candidates, with admissibility checked independently."""
        positive(x=x,quota=quota)
        B=self.S-2*quota/(self.k*x); C=quota/self.c; disc=B*B-4*C
        if B <= 0 or disc < 0: return []
        large=(B+math.sqrt(disc))/2
        roots=[large] if disc == 0 else [large,C/large]
        result=[]
        for g in roots:
            E=quota*(g+self.rho)/(self.a*g)
            try:
                branch=replace(self,E=E); state=branch.reconstruct(x)
                admissible=min(state.values()) > 0 and abs(state['jG']-quota) < 1e-8*quota
            except ValueError: admissible=False
            result.append(dict(g=g,E=E,admissible=admissible))
        return result


@dataclass(frozen=True)
class ThioredoxinBranch:
    T: float; z0: float; E: float; a: float; b: float; c: float; d: float; e: float; k: float

    def __post_init__(self):
        positive(T=self.T,E=self.E,a=self.a,c=self.c,d=self.d,e=self.e,k=self.k)
        if any(not math.isfinite(float(v)) or v < 0 for v in [self.z0,self.b]): raise ValueError('Invalid baseline/hyperoxidation rate')
        self.response

    @property
    def rho(self): return 1/self.a+1/self.d+self.b/(self.c*self.d)
    @property
    def S(self): return self.T-self.z0
    @property
    def response(self):
        v=self.rho*self.e
        return QuadraticResponse(v,1-self.S*v,self.E*self.e/self.k,self.S,self.E*self.e,v,1)

    def reconstruct(self,x):
        positive(x=x); y=self.response.carrier(x); j=self.response.current_from_carrier(y)
        return dict(y=y,zT=self.z0+j/(self.k*x),r=j/self.a,h=j/self.d,
            w=self.b*j/(self.c*self.d),v=j/(self.e*y),jT=j)


class DecreasingSource(Protocol):
    """Contract: strictly decreasing, positive below N, zero at N."""
    N: float
    def unit_current(self,x): ...


@dataclass(frozen=True)
class Regeneration:
    P: float; P0: float; K: float; V: float
    def __post_init__(self):
        positive(P=self.P,K=self.K,V=self.V)
        if not math.isfinite(float(self.P0)) or not 0 <= self.P0 < self.P: raise ValueError('Require 0 <= P0 < P')
    @property
    def N(self): return self.P-self.P0
    @property
    def D(self): return self.K+self.P
    def unit_current(self,x): return self.V*(self.N-x)/(self.D-x)
    def derivative(self,x): return -self.V*(self.D-self.N)/(self.D-x)**2


@dataclass(frozen=True)
class InhibitedSource:
    """Illustrative product inhibition of source; not inhibition of a branch."""
    base: Regeneration
    KI: float
    def __post_init__(self): positive(KI=self.KI)
    @property
    def N(self): return self.base.N
    def unit_current(self,x): return self.base.unit_current(x)/(1+x/self.KI)


@dataclass(frozen=True)
class ServiceDesign:
    responses: tuple[QuadraticResponse,...]
    source: DecreasingSource

    def __post_init__(self):
        if not self.responses: raise ValueError('At least one branch required')

    def demand(self,x): return sum(r.current(x) for r in self.responses)

    def equilibrium(self,scale):
        positive(scale=scale)
        return brentq(lambda x:scale*self.source.unit_current(x)-self.demand(x),0,float(self.source.N),xtol=1e-14)

    def decision(self,quotas):
        if len(quotas) != len(self.responses): raise ValueError('One positive quota per branch')
        for q in quotas: positive(quota=q)
        ceilings=[r.current(self.source.N) for r in self.responses]
        if any(q >= c for q,c in zip(quotas,ceilings)):
            return dict(repairable=False,ceilings=ceilings,reason='A quota reaches a finite-pool ceiling')
        inv=[r.invert(q) for r,q in zip(self.responses,quotas)]; L=max(x for _,x in inv)
        currents=[r.current(L) for r in self.responses]; R=self.source.unit_current(L)
        return dict(repairable=True,ceilings=ceilings,carrier_floors=[y for y,_ in inv],
            nadph_floors=[x for _,x in inv],L=L,currents=currents,
            minimum=sum(currents)/R,quota_only=sum(quotas)/R,
            overdelivery=[j-q for j,q in zip(currents,quotas)],
            isolated=[q/self.source.unit_current(x) for q,(_,x) in zip(quotas,inv)])

    @staticmethod
    def robust_command(scale,tolerance):
        if not 0 <= tolerance < 1: raise ValueError('Tolerance in [0,1) required')
        return scale/(1-tolerance)


@dataclass(frozen=True)
class MaintainedKinetics:
    """Literal eight-state ODE. Equilibrium response curves are NOT a transient closure."""
    gpx: GlutathioneBranch
    trx: ThioredoxinBranch
    source: Regeneration
    coordinates=('x','z','e1','e2','zT','h','w','v')

    @property
    def design(self): return ServiceDesign((self.gpx.response,self.trx.response),self.source)

    def full(self,u):
        x,z,e1,e2,zT,h,w,v=u
        return dict(zip(self.coordinates,u),g=self.gpx.G-2*z-e2,e0=self.gpx.E-e1-e2,
            y=self.trx.T-zT,r=self.trx.E-h-w-v,nadp=self.source.P-x)

    def state(self,x):
        s=dict(x=x,**self.gpx.reconstruct(x),**self.trx.reconstruct(x))
        return np.array([s[k] for k in self.coordinates])

    def currents(self,u,scale):
        s=self.full(u); G,T=self.gpx,self.trx
        return dict(source=scale*self.source.unit_current(s['x']),phi0=G.a*s['e0'],
            phi1=G.b*s['g']*s['e1'],phi2=G.c*s['g']*s['e2'],psiG=G.k*s['x']*(s['z']-G.z0),
            thetaR=T.a*s['r'],thetaW=T.b*s['h'],thetaH=T.c*s['w'],thetaV=T.d*s['h'],
            thetaY=T.e*s['y']*s['v'],psiT=T.k*s['x']*(s['zT']-T.z0))

    def rhs(self,t,u,scale):
        c=self.currents(u,scale)
        return np.array([c['source']-c['psiG']-c['psiT'],c['phi2']-c['psiG'],
            c['phi0']-c['phi1'],c['phi1']-c['phi2'],c['thetaY']-c['psiT'],
            c['thetaR']+c['thetaH']-c['thetaW']-c['thetaV'],c['thetaW']-c['thetaH'],c['thetaV']-c['thetaY']])

    def margins(self,u):
        s=self.full(u)
        return np.array([s[k] for k in self.coordinates]+[s[k] for k in ['g','e0','y','r']]+
            [self.source.N-s['x'],s['z']-self.gpx.z0,s['zT']-self.trx.z0])

    def jacobian(self,u,scale):
        return np.column_stack([self.rhs(0,np.asarray(u,dtype=complex)+1e-20j*np.eye(8)[i],scale).imag/1e-20 for i in range(8)])

    def integrate(self,initial,scale,times):
        positive(scale=scale)
        if self.margins(initial).min() < -1e-12: raise ValueError('Initial state outside conserved physical domain')
        times=np.asarray(times)
        if len(times)<2 or times[0]!=0 or not np.all(np.diff(times)>0): raise ValueError('Increasing times starting at zero required')
        sol=solve_ivp(lambda t,u:self.rhs(t,u,scale),(0,times[-1]),initial,t_eval=times,
            method='Radau',rtol=2e-9,atol=2e-11,jac=lambda t,u:self.jacobian(u,scale))
        if not sol.success: raise RuntimeError(sol.message)
        if min(self.margins(u).min() for u in sol.y.T) < -1e-7: raise RuntimeError('Numerical physical-domain violation')
        return sol.y.T


def nominal(exact=False):
    conv=lambda d:{k:F(str(v)) for k,v in d.items()} if exact else d.copy()
    # Literal paper parameters: deliberately independent of the editable block.
    G=GlutathioneBranch(**conv(dict(G=371.56,z0=1.78,E=50,a=.21,b=.04,c=10,k=3.2)))
    T=ThioredoxinBranch(**conv(dict(T=.505,z0=.075,E=19.096,a=.4,b=.00072,c=.003,d=15,e=2.1,k=20)))
    R=Regeneration(**conv(dict(P=30.3,P0=.3,K=57,V=375)))
    return MaintainedKinetics(G,T,R)