# Structural autocatalysis and productive reactor operation

Companion to paper 16, *Structural autocatalysis can be common while
productive operation is rare: an exact separation theorem in a capped-Zipf
binary polymer model*. Manuscript SHA-256:
`76c780846acad0324c5ac76cc00fb7d42ca5c5bf7060f9b9a2d13b59f75f5f3f`.

This package connects a literal polymer catalogue, catalytic environments,
RAF detection, and both deterministic and stochastic fed-reactor models.
It also evaluates the capped source law and the paper's conservative finite
probability bounds without allocating astronomical catalogues. The scientific
layers share the same reaction identities. No universal framework is required.

## Run

Python 3.11 or newer:

```sh
python -m venv .venv
# PowerShell: .venv\Scripts\Activate.ps1
# macOS/Linux: source .venv/bin/activate
python -m pip install -r requirements.txt
python example.py --output outputs
python -m unittest -v
```

Defaults finish in roughly tens of seconds on the development machine.
Pinned NumPy, SciPy, mpmath and Matplotlib provide random generation,
integration, high-precision special functions, and figures. Seeded outputs
were checked with Python 3.11 on Windows. Scientific JSON and CSV outputs
replay identically in the checked dependency environment; exact bytes across
other library/platform versions are not promised.

## Inputs and model

Edit the clearly marked block at the top of `example.py`. The default n=4
catalogue contains 30 nonempty binary words and 68 reversible ligation/split
channels. A channel keeps its ordered factors and split position: two splits
of a repetitive word remain distinct. Food is all monomers and dimers.
`PolymerCatalogue.closure` ignores catalyst timing and uses both directions;
`CatalyticEnvironment.max_raf` then removes channels without internally
available reactants/products and catalysts. Its result is the largest RAF,
not a smallest one. `irreducible_raf` returns a deletion-minimal witness,
which need not have globally minimum size.

`CappedZipfSource` samples a degree D=min(K,R)-1 from the full capped Zipf law,
including its tail atom at R-1, then a uniform D-subset of channels for each
molecule. Rows are independent; entries within a row are not. An incidence
catalyses both directions. Degree probabilities are evaluated at high
precision and converted to double precision for categorical sampling.
`statistics(n)` instead evaluates the exact finite-law formulas at 80-digit
working precision, returning 35 significant digits; it does not enumerate
the catalogue. In particular p=E[D]/R, not E[D]/(R-1).

`KineticParameters` draws independent marks S,B,H in {1,1.5,2}, with basal
coefficient epsilon*S*B and catalytic coefficient 4*S*H. The same S is used
across a channel's directions and catalysts. Epsilon=1/500000000 is tiny but
strictly positive. The special figure-mark adapter exactly reproduces the
manuscript's specified host, with one catalyst per selected channel. Use
the default general sampler for multiple catalysts per channel.

`ReactionEvent` retains the full input multiplicity of a catalyst even when
it coincides with a substrate or product. `FedReactor.rates(counts,V)` uses
falling factorials with no factorial divisor: three identical inputs give
N(N-1)(N-2), not N^3 and not a binomial coefficient. Catalysts cancel only
from net stoichiometry. Feed supplies each food molecule at rate V; every
species dilutes at rate N. Food-only initial counts equal V per food species.
Time is measured in dilution-time units and concentrations are counts/V.
These are effective dimensionless manuscript parameters, not calibrated
chemical rates or a physical vessel-volume conversion.

`deterministic` replaces falling-factorial count propensities by concentration
mass-action monomials and integrates the full finite host with DOP853. Its
augmented variables integrate export, signed catalytic input, signed basal
input and positive basal input. `stochastic` is the direct Gillespie method
on the literal feed, outflow, basal and catalytic event list; it neither
uses tau-leaping nor substitutes an ODE. The terminal incomplete waiting
interval produces no event. Every jump checks integer mass conservation
and the complete signed nonfood balance.

## What the default calculations show

The specified five-incidence host is deliberately chosen, not sampled from
the capped-Zipf law. All food rows are empty and the product 0011 catalyses
00+11 <-> 0011, so it belongs to the paper's witness event E and contains a
singleton RAF. Four additional nonfood catalytic assignments remain active.
Deterministic enabled export over (1,100] is about 162.69303; disabling only
catalysis gives 0.0000414748. Continuous feed allows cumulative export to
exceed instantaneous reactor mass. Two integration tolerances and the total
and nonfood mass balances check this numerical illustration.

The same host at finite volume 20 gives zero output in both saved stochastic
runs. With such tiny basal rates, a small count reactor may never create a
seed during the observation window. Those two trajectories are not an
estimate of source-averaged output probability, nor a contradiction of the
large-volume theorem. Enabled and disabled runs reuse the source and kinetic
parameters. Reusing the random seed aids reproducibility but does not pair
individual reaction events once hazards diverge.

The productive observable is the conjunction of maximum total mass <= 11
through time 100 and exported nonfood mass/V > 0.1 during (1,100]. It is
checked from every stochastic jump, even though saved traces are decimated.
Hitting `MAX_EVENTS` returns `completed=False` and `productive=None`; never
count an unfinished trajectory as a failed run. Stochastic volume has an
implementation guard below 10^12. The theorem's enormous sufficient volumes
are handled analytically, not simulated.

The independent source census samples 48 environments at each n=4,...,8;
all 240 saved samples contain RAFs. This small finite census does not imply
probability one or establish the limiting survival probability. Saved sizes
are maximum-RAF sizes. The analytic source plot shows X_n*p_n and E[D]/n
approaching 9/pi^2, and the separate probabilities p_n and q0^6*p_n decaying.
The latter is the exact witness-environment probability, not the full
productive probability. The openness q*=1-exp(-9/pi^2) is also not the RAF
survival probability theta(q*).

`TheoremBounds.structural_budget` evaluates the paper's explicit B_n=A_n^2
in logarithmic form. Its normalized log tends to zero, illustrating the
subexponential cutoff. The theorem that this cutoff captures asymptotically
all RAF-existence mass is an analytic result, not proved by this finite census.
It does not determine typical minimum RAF size or the sizes in most
productive environments.

`TheoremBounds.evaluate` is tied to the manuscript's fixed rates, feed,
window and observable. It does not certify modified reactor parameters.
It rejects V<2n/d, reports the noise term in log space, clips probability
upper bounds at one, and marks a nonpositive lower bound as absent. At the
sampled volume its hypotheses fail. At V=10^60(n+1)^2 the lower bound is
positive and noise negligible, but the catalogue upper constant is so large
that the saved finite-n upper bound remains one. These are sufficient proof
bounds, not a necessary or practical vessel size. No exact productive
probability or asymptotic prefactor is claimed.

## Reuse a sampled environment

```python
import numpy as np
from example import PolymerCatalogue, CappedZipfSource, KineticParameters, FedReactor

catalogue = PolymerCatalogue(4)
source = CappedZipfSource(catalogue)
environment = source.sample(np.random.default_rng(123))
print('Contains a RAF:', bool(environment.max_raf()))
parameters = KineticParameters.sample(environment, seed=456)
reactor = FedReactor(environment, parameters)
run = reactor.stochastic(volume=50, seed=789)
print(run['completed'], run['productive'])
```

To explore witness environments, call `sample(rng, condition_on_witness=True)`.
This conditions food rows to be empty and samples the product row using the
correct size-biased degree law, then forces its required channel. Simply
inserting that channel into an ordinary sample would give the wrong law.
Conditioned trajectories estimate conditional quantities only; weighting
their successes by q0^6*p_n supplies a witness-event contribution, not all
productive environments.

For a volume sweep, reuse one environment and its kinetic parameters and
call `stochastic` at each integer volume with independent seeds. For complete
trials estimating the joint source/marks/trajectory law, resample all three
layers. Repeated trajectories in one environment answer a different question.
Changing epsilon can help explore startup, but then the paper's quantitative
bounds no longer apply automatically. Explicit catalogues are guarded at
1022 species; large-n source-law and bound calculations need no catalogue.

## Saved evidence

- `deterministic.csv` and `disabled_deterministic.csv`: all species and rewards.
- `enabled_stochastic.csv` and `disabled_stochastic.csv`: decimated traces;
  full-jump maxima and exact balance checks are in `results.json`.
- `results.json`: solver diagnostics, finite source census, source-law values,
  subexponential cutoff and finite theorem bounds.
- `reactor.png/.svg` and `source.png/.svg`: generated scientific figures.
- `console.txt` and `run_metadata.json`: summary and source/output hashes.

Six test groups check catalogue identities and literal RAF subsets, the cap
atom against direct sums, conditional source sampling, coincident-reactant
propensities, mass drift, the manuscript's numerical host, exact jump rewards,
event-budget handling, and finite-bound hypotheses. The PDF and its formal
proofs are unchanged; this Python implementation is not extracted from Lean.

License: MIT is proposed, pending the owner's choice. No license grant is
asserted by this package.