# Same reaction skeleton, different mass-action stability

This example implements the constructions and exact certificates in *Kinetic
order is invisible to D-cores: a support-preserving mass-action lift and an exact
counterexample*. It compares the original four-species network, the manuscript's
exact-eigenpair padding, and its smaller maximum-exponent-200 padding. It also
provides a general rational lifting routine and a reusable mass-action model.
The manuscript hash is recorded in the source. This is original Python from the
paper's definitions, not extracted Lean code.

## Model and inputs

Top-level inputs specify the reactant and product matrices (species rows,
reaction columns), stationary flux, equilibrium concentrations, padding orders,
and integer-order sweep. Values are abstract nondimensional mathematical
witnesses, not realistic calibrated reaction mechanisms. Large orders are
essential to these particular examples; no low-molecularity claim is made.

The skeleton consists of S=P-Y and the pattern of positive entries in Y.
Adding the same nonnegative integer padding C to Y and P, only at existing
reactant incidences, preserves that skeleton. It preserves every child-selection
matrix: choose distinct reactions for selected species that occur as reactants,
then take the corresponding square submatrix of S. It changes the mass-action
Jacobian because reactant multiplicities are the kinetic exponents.

`Network` validates integer complexes, applies padding, and enumerates children.
`EquilibriumModel` takes a network, positive state xbar and positive stationary
flux v satisfying S v=0. It defines k[j]=v[j]/product(xbar[i]**Y[i,j]), making
xbar an equilibrium, and computes the exact rational reactivity
R[j,i]=v[j]Y[i,j]/xbar[i] and Jacobian S R. `RationalLift` implements the general
common-denominator construction to reproduce any admissible rational target R
on a given skeleton and stationary positive flux.

Rate constants can be unimaginably large or small. The model therefore stores
their reconstruction formula and reports log(k), while evaluating the equivalent
rates r(x)=v*exp(Y.T@log(x/xbar)). It never expands the huge exact powers.
Numerical rate evaluation is for positive states; overflow raises an exception
rather than silently clipping rates. Extreme orders still make off-equilibrium
evaluation sensitive. Exact equilibrium certificates do not depend on floating
point or on successful nonlinear integration.

## Run

Python 3.11 or newer. Direct dependencies are pinned in requirements.txt;
transitive dependencies are not fully locked.

Windows PowerShell:

```powershell
py -m venv .venv
.\.venv\Scripts\python.exe -m pip install -r requirements.txt
.\.venv\Scripts\python.exe example.py --output outputs
.\.venv\Scripts\python.exe -m unittest -v
```

POSIX:

```sh
python3 -m venv .venv
.venv/bin/python -m pip install -r requirements.txt
.venv/bin/python example.py --output outputs
.venv/bin/python -m unittest -v
```

The default entry command checks the fixed paper certificates, so changing the
principal witness will deliberately fail its exact-eigenpair assertion. Use the
imported model for independent networks and target reactivities, or edit the
compact exponent and sweep to explore that controlled family.

```python
from example import Network, EquilibriumModel, RationalLift
from fractions import Fraction
import numpy as np

# Reversible A -> B and B -> A: same positive flux through both reactions.
network = Network([[1, 0], [0, 1]], [[0, 1], [1, 0]])
model = EquilibriumModel(network, state=[2, 3], flux=[2, 2])
trajectory = model.simulate([1, 4], np.linspace(0, 4, 41))
print(model.jacobian())

target = [[Fraction(1, 3), 0], [0, Fraction(5, 7)]]
lifted = RationalLift.build(network, [2, 2], target)
assert network.same_skeleton(lifted.network)
assert np.array_equal(lifted.reactivity(), target)
```

`simulate` solves the nonlinear equations in log coordinates with Radau and an
analytic transformed Jacobian; it returns concentrations as species-by-time.
It is tested against a closed-form reversible-pair trajectory and its material
balance. The paper's compact witness is extremely stiff, so long nonlinear runs
are optional. The default figures describe local equilibrium stability and do
not require nonlinear trajectories. Exact determinant and child enumeration
routines are intended for small networks and have combinatorial cost.

## What is checked

* Both padded witnesses and the uniform rational lift retain the original
  stoichiometry, reactant support, and 24 child selections.
* Exact weighted symmetrizations and principal minors cover the child matrices
  and their restrictions. A separate exact scaled-polynomial identity handles
  the singular maximal child, whose eigenvalues include zero for every positive
  diagonal scaling. Thus the appropriate claim is D-noninstability, not strict
  D-stability for every child. Coverage is checked without assuming that each
  child belongs to only one maximal selection.
* For the original network, the isolated negative eigenvalue and the lower
  cubic's coefficient polynomials are computed from the actual Jacobian template.
  The cubic Hurwitz gap expands into seven strictly positive monomials. This
  is a certificate for every positive column scaling, not a sampled claim.
* The principal padded witness has exact eigenvalue 1/500 + (51/500)i. The code
  checks both real and imaginary vector equations with Fraction arithmetic.
* The compact witness at (1,235,585,1/257) reproduces the exact quartic and
  Delta2=2196782093181776/178929 > 0,
  Delta3=-384290295376441702400/327976857 < 0. With positive coefficients these
  Routh signs give two right-half-plane roots. Its unstable pair is numerically
  about 0.000469 +/- 0.2410i.

The sweep changes the two compact-family orders together through even integers
4..240, keeping xbar and v fixed and reconstructing the rate constants at each
point. The first sampled unstable point is 196. This is neither a globally
minimal exponent nor a certified continuous bifurcation threshold. Marker
classifications use exact Hurwitz signs; plotted spectral positions use floating
point. Lines merely connect the sampled integer networks.

The six test groups check the paper's exact quartic values, independent
determinant evaluations, the principal eigenpair, lifting on a second network,
certificate coverage, the seven-term universal stability identity, a complex-step
rate derivative, the analytical reversible-pair solution, and invalid inputs.

## Outputs and evidence boundary

`models.json` exports both sides of every reaction, states, fluxes, reactivities,
Jacobians, log rate constants and reconstruction formulas. `certificates.json`
contains rational eigenpair and child/stability certificates. CSV files give
all child selections, spectra and sweep results. PNG/SVG figures, console output,
a summary, and run metadata complete the package. Metadata records source and
manuscript hashes, dependencies, elapsed time and output hashes. No random
sampling is used; SVG metadata can differ between equivalent reruns.

The exact certificates are numerical algebra plus the stated dissipativity and
Routh--Hurwitz criteria. The general lifting and stability theorems remain those
of the paper. No Lean build is run here; the compact witness is explicitly
outside the manuscript's compiled theorem. Complex-pair instability does not
establish a Hopf bifurcation or a periodic orbit. The examples are not evidence
that bounded-molecularity versions of the conjecture fail.

## License

MIT is proposed for this original example code, pending the owner's license
choice. No license grant is asserted here. This does not change the license of
the manuscript or any third-party dependency.