# Repeated harvesting from an autocatalytic reactor

Withdrawal removes catalyst as well as product. This example follows the paper's
six reversible reaction pairs through conditioning and repeated withdrawal,
food refill, recovery, and collection. It asks whether the actual remaining
chemical mixture recovers, how much free template it exports, and how much of
that production must have been newly synthesized.

The source is *Productive recovery and repeated harvesting in a reversible
autocatalytic reactor*, 15 September 2026: Table 1, the maintained-flow equations,
the intervention and operation definitions, and the main theorem. The manuscript
SHA-256 is recorded in `example.py`. This implementation was written from those
equations; it does not redistribute the author's diagnostic script or Lean code.

## Model and inputs

The state contains foods U and W, free template X, complexes C1 and C2, and duplex
Z. The reaction path X + U -> C1, C1 + W -> C2, C2 -> Z, Z -> 2X is reversible.
Basal formation U + W <-> X and fuel-driven cleavage X <-> U + W complete the
model. Fuel and waste have maintained unit activity. Unit food inflows and unit
washout act continuously on the six internal species.

`Reactor` builds the concentration derivative from `ReactionPair` objects.
`Intervention` removes well-mixed material, applies additional species losses,
and adds food only. `Integrator` advances that state and four integral counters.
`operate` links the actual endpoints and accounts for every pulse. The callable
protocol receives a copy of the actual state, so feedback can be introduced
without changing chemistry or integration. Plotting consumes saved results.

Edit **USER INPUTS** at the top of `example.py`:

| Input | Meaning and admissible model domain | Paper default / conditions |
| --- | --- | --- |
| RELEASE_SPEED | Nonnegative normalized release rate | 19; theorem 19 to 21 |
| CLEAVAGE_SPEED | Nonnegative normalized cleavage rate | .04; theorem .02 to .04 |
| INITIAL_STATE | Six nonnegative normalized concentrations | Paper's explicit operating-region state |
| RETAINED_FRACTIONS | Fraction retained after withdrawal, 0 to 1 | .25, .75, .4, .6; theorem .25 to .75 |
| LOSS_PATTERNS | Fraction surviving additional handling, per species | Paper's two diagnostic patterns; theorem .98 to 1 |
| REFILL_ERRORS | Food refill deviation; actual refill must be nonnegative | Alternating -.005 and .005; theorem absolute error <= .005 |
| CYCLES | Positive integer cycle count | 32, extending the paper's 8-cycle diagnostic |
| Times | Positive normalized durations | 12 conditioning; 3 recovery; 1 collection |
| RTOL / ATOL | Positive integration error targets | 1e-9 / 1e-12 |

All model values are normalized and uncalibrated. The paper illustrates a reading
of one concentration unit as 1 mM, one time unit as 60 seconds, and reactor volume
as 1 mL. That makes one amount unit 1 micromole. It does not identify a laboratory
chemistry. For a reaction of molecularity m, the dimensional rate coefficient is
the model coefficient divided by `time_unit * concentration_unit**(m-1)`; changing
axis labels alone does not perform that conversion.

## Run offline after setup

Python 3.11 or newer. `requirements.txt` pins direct dependencies, not transitive
dependencies. The published run's environment is in `run_metadata.json`.

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
```

`cycles.csv` separates window output, food use, service, and net synthesis.
`trajectories.csv` stores all six species and material/catalytic observables.
`pulses.csv` records both sides of each instantaneous intervention.
`retention_sweep.csv` compares independent cycles from the same preparation at
five retention values. PNG and SVG figures use these same results. `summary.json`
and the actual `console.txt` report inputs and numerical checks. Metadata binds
source, manuscript, environment, runtime, and output hashes.

## Change the operating policy through the public interface

```python
from example import Reactor, Rates, Integrator, Intervention, Schedule, operate, Y

def feedback(cycle, state):
    return Intervention(retained=.75 if Y @ state < .2 else .25)

run = operate(Integrator(Reactor(Rates(release=21, cleavage=.02))),
              schedule=Schedule(cycles=8), protocol=feedback)
print(run.summary)
```

The rates, initial preparation, intervention box, and time schedule determine
whether the paper's bounds apply. Outside them, valid nonnegative inputs still
define an exploratory simulation, but no theorem bound is asserted. For more
complex chemistry, supply a reactor with the same derivative interface and
revisit every observable and balance; the fixed paper's guarantees do not transfer
automatically to new reactions.

## What the calculation establishes

The paper guarantees at least 1/28 template equivalents and 1/540 free X in each
collection window, with at most 951/200 of each food and 9/50 gross fuel/waste
service per routine cycle. Simulated values can exceed these conservative floors.
The program is a numerical illustration, not a uniform proof over the parameter
box, a finite-molecule survival probability, or a test of periodic attraction.

Template inventory counts X, C1, C2 once and Z twice. Its derivative is net
covalent synthesis minus continuous effluent. Consequently total synthesis equals
final minus initial inventory plus all effluent, pulse withdrawals, and handling
losses. Collected product counts only the final unit of each routine cycle;
withdrawals and noncollection effluent are not silently counted as collection.
The paper's synthesis lower bound `m/28 - 1.1` becomes positive at cycle 31.

Tests independently substitute the literal field, check elemental conservation,
pulse and flow accounting, low-stock recovery, invalid inputs, continued states,
and a counterexample to globally positive catalyst growth. Each run compares a
Radau cycle with stricter DOP853 integration. The linear balance residual alone
does not measure integration accuracy. Positivity checks inspect solver nodes and
saved samples; they are not interval enclosures of all intermediate times.

No repository-wide software license was present when this example was prepared.
MIT is proposed for this new original code, pending the owner's licensing choice;
no license is applied to the manuscript or other authors' code by this package.