# Reliable repeated harvesting in a finite-molecule autocatalytic reactor

This package implements the literal molecule-count reactor and molecular withdrawal law, continues actual returned states, and evaluates the paper's repeated-operation certificate. It follows manuscript entry 43, SHA-256 `235f5a2a613e48b9059f90bf38afa0b502225d32956e541096cd3141b2ea4db2`.

## Run and edit

Use Python 3.11 or newer:

```sh
python -m venv .venv
# Activate the environment for your shell.
python -m pip install -r requirements.txt
python example.py --output outputs
python -m unittest -v test_example
```

Inputs at the top of `example.py` control the certified copy scale and mission length, target failure, paired kinetic speeds, initial density, withdrawal schedule, species losses, refill errors, density cycles, a separate small-copy stochastic diagnostic and physical units. Keep rational inputs as strings for exact coefficient and rounding calculations. The default stochastic diagnostic has V=1000 and three cycles, whereas the certificate uses V=200,000,000,000 and 100 cycles. They are intentionally different calculations, labeled separately throughout.

The main run ordinarily takes several seconds. Direct simulation at the certificate scale would require far too many events. There is no tau-leaping substitution: `StochasticCycle` uses the chronological exponential-waiting-time/reaction-label algorithm with floating-point propensities. A budget limit returns an unfinished record with no success verdict. The limited multinomial sampler rejects individual counts at or above 2^53; exact rational pulse enumeration is available for tiny populations.

## Chemical source and reusable components

`Reactor` builds twenty directed channels for six internal species: foods U,W, template X, complexes C1,C2 and duplex Z. Twelve channels form these reversible pairs:

```text
U+W <-> X           epsilon, epsilon/10
X+U <-> C1          20, 20
C1+W <-> C2         20, 20
C2 <-> Z           20, 2
Z <-> 2X           r, r
X+F <-> U+W+P      d, d*eta, with maintained F=P=1
```

Here epsilon=1/500,000,000 and eta=1/8,000,000,000. The theorem uses fixed r in [19,21] and d in [1/50,1/25] throughout a mission. The remaining channels feed U and W at rate V each and wash each internal species out at its molecule count. Each `Channel` stores reactant labels, product labels, coefficient, stoichiometric jump and event marks. The count propensity for 2X association is r*nX*(nX-1)/V, with no extra factorial divisor; its concentration counterpart is r*x^2.

The source freshly verifies exact symbolic identities for the two material totals and the weighted catalyst:

```text
A = U+X+2C1+2C2+2Z          L A = V-A
B = W+X+C1+2C2+2Z          L B = V-B
Y = X+(9/8)C1+(7/5)C2+(9/5)Z
L Y = V*Y(f(N/V)) + r*N_X/(5V).
I = X+C1+C2+2Z             template-equivalent inventory
```

The last correction is the favorable finite-count falling-factorial term. The drift identities do not imply that individual stochastic trajectories grow monotonically or stay in a corridor.

## Restart, pulse and actual history

The restart set requires A and B in [159V/160,161V/160] and Y>=V/20. The integer comparison uses 40Y=40X+45C1+56C2+72Z. All-free template (0,0,V,0,0,0) is admitted. At multiples of 1000, the richer state proportional to (934,935,60,1,1,1) is admitted and has I=8V/125. The main certificate uses this richer state. This is established-stock preparation, not food-only startup; no twelve-unit conditioning period is prepended.

An `Intervention` specifies q in [1/4,3/4], six survival factors ell_i in [49/50,1], and two refill errors of magnitude at most 1/200. Each molecule of species i independently enters one of three categories with probabilities

```text
retained: q*ell_i       withdrawn: 1-q       extra loss: q*(1-ell_i).
```

Only retained molecules remain. The pulse then adds floor(V*(1-q+eU)) food-U molecules and floor(V*(1-q+eW)) food-W molecules. It never adds template or complexes. The pulse is instantaneous, and V is maintained externally. `MolecularPulse.sample()` uses a multinomial draw for each species; `enumerate()` returns the complete rational product law for at most twelve molecules. Outcomes losing all catalyst keep their mass. Tests verify normalization and that extinction probability explicitly.

`StochasticCycle` runs this pulse followed by the four-unit source. `HistoryMission` calls a user-supplied controller with the complete recorded history, then starts the next cycle from the actual endpoint. It continues after a failed completed cycle, without rejection, resampling or resetting to a representative restart state. Controller observations are copied so accidental mutation cannot rewrite stored history. The example `ReturnedStockController` switches to high retention after a zero-free-output cycle and otherwise follows the configured schedule. Its decisions remain inside the admitted intervention box. Reaction speeds do not change within or between cycles.

## Same-event output and supply ledgers

During [3,4], washouts of X,C1,C2,Z add 1,1,1,2 to collected template equivalents QI. Only a washout of X adds to free-template output QX. Thus QX is included in QI; adding them would double count free template. Before time 3 these output marks are off, although washout and chemistry continue naturally. An event exactly at 3 is collected; the continuous-time law assigns such an exact event zero probability.

Food counters CU,CW include pulse doses and all feed arrivals over [0,4]. Gross driven service G counts both pair-5 directions, not net fuel consumption. `CycleLedger` also records all-time template washout and signed internal synthesis. Thresholds do not stop or cap the process. A completed record succeeds exactly when:

```text
endpoint belongs to the restart set
QI >= ceil(V/56), QX >= ceil(V/1080)
CU,CW <= 5V, G <= floor(V/5).
```

No corridor or favorable-pulse conditioning is imposed on the simulated law. An unfinished record cannot be classified as successful merely because its partial counters look good.

## Repeated-operation certificate and sizing

`MissionCertificate` requires integer V>=1,000,000 and fixed paired speeds in the theorem rectangle. It evaluates the full Appendix A error budget e(V), preserving repeated material terms. It does not substitute the convenient envelope for the full sum or silently drop exponentials outside ordinary floating-point range. The calculation uses high-precision logarithms and outward interval arithmetic.

From any admitted restart, one-cycle success is at least max(0,1-e(V)). Conditioning successively on actual successful returns gives the joint lower bound [max(0,1-e(V))]^m; it does not require independent cycles. The code clamps uninformative bounds to zero, including when e(V)>1. An empty mission has probability one. The nine-decimal printed claim is rounded downward and checked against the interval endpoint. The probability theorem and its chronological-process identification are imported from the manuscript, not proved by evaluating the formula.

For V>=200,000,000,000 the paper gives e(V)<=101*exp(-V/10^10). The default full budget certifies 100-cycle success >=0.999979388 and 48,000-cycle success above 0.99. The logarithmic rule

```text
V >= max(200,000,000,000, ceil(10^10*log(101*m/failure)))
```

certifies a mission at the requested failure tolerance. `sizing()` takes the ceiling of an outward interval upper endpoint and then checks a rational Taylor lower sum for exp(V/10^10) against 101*m/failure. For one million cycles and failure 1/100, this gives 230,358,012,608. The paper's displayed 230,359,000,000 is a coarser outward-rounded sufficient example, not a disagreement.

Optional demand inputs add sufficient scales ceil(56*template_demand/m) and ceil(1080*free_demand/m). These are sufficient demand rules; exploiting every possible integer-ceiling plateau is not an optimization objective. The earlier linear sizing rule is also exported for comparison. A lower bound e(V)>=100*exp(-V/10^10) identifies the bottleneck of this certificate, not a physical minimum reactor size.

The coefficient 100 in the dominant term comes from a low-free-template occupation estimate. `PhaseTransport` checks the exact adjoint recurrence and evaluates its closed-form finite-clock weights. With qc=3000V, a=1-70/qc, the weights after n steps are

```text
a^n * (1, 20n/(qc*a), 580n(n-1)/(qc^2*a^2), 38n/(qc*a)).
```

For 89V<=n<=91V, each coordinate is at least 1/35 of the corresponding Y weight. Thus catalyst initially bound in any phase can fund later free-template observations. Rational window margins and an explicit Taylor upper bound for exp(11/5) are checked. The plotted V=100 weights illustrate this algebra; they are not a probability claim for a V=100 reactor.

`stopping_boundary_example()` supplies a separate exact three-state chain illustrating why a closed restart event cannot simply be transferred from a stopped process. After two steps the stopped closed-event probability is one, while the physical value is one half; the strictly active probability is one quarter under both laws. This is not a reactor counterexample. The paper's strict-activity and renewal/exhaustion arguments supply the actual process-identification theorem; Python does not reprove them.

## Internal synthesis and physical reading

Every sampled cycle freshly audits the pathwise balance

```text
net internal synthesis = final I - initial I
                       + all template washout
                       + withdrawn I + extra-lost I.
```

Chemical pairs 0 and 3 contribute +1 forward/-1 reverse, pair 5 contributes -1 forward/+1 reverse, and the other chemical pairs contribute zero. Food feeds add no template inventory. Collection is only part of total washout, and pulse withdrawals are not falsely counted as collected output.

On m successful cycles, final restart gives I_final>=ceil(V/28). Therefore net synthesis is at least m*ceil(V/56)+ceil(V/28)-I_initial. Uniformly over admitted starts, replace I_initial by floor(161V/160). This becomes positive after 55 successful cycles; the richer default start already gives a positive lower bound after two. At V=200 billion and m=100:

| Quantity | Conditional bound |
|---|---:|
| Collected template equivalents | >=357,142,857,200 |
| Included free template molecules | >=18,518,518,600 |
| Each food expenditure | <=100,000,000,000,000 |
| Gross driven events | <=4,000,000,000,000 |
| Final template inventory | >=7,142,857,143 |
| Net synthesis, arbitrary admitted start | >=163,035,714,343 |
| Net synthesis, richer configured start | >=351,485,714,343 |

This proves net synthesis on the successful event; it does not identify which exported molecules came from which pathway or establish causal knockout necessity. Costs are allowances on successful missions, not bounds on spending after failure.

At 1 mM and a one-minute time unit, V=200 billion corresponds to about 0.332108 nL. One hundred cycles require 400 minutes of source-running time, excluding preparation and handling. Guaranteed collected output is about 0.593050 pmol of template equivalents, including 0.0307507 pmol of free template. The chemistry and units are schematic, not a calibrated biochemical protocol. Maintained reservoir activities, instantaneous handling and the thinning law require physical justification in any implementation. Gross driving counts are not a complete energy bill; separation, purification, pumping and maintenance are outside the model.

## Reuse, diagnostics and verification

```python
from example import Reactor, MissionCertificate, HistoryMission, StochasticCycle, Intervention

reactor = Reactor(release='19', cleavage='1/25')
certificate = MissionCertificate(200_000_000_000, reactor)
print(certificate.evaluate(48_000))
print(certificate.sizing(1_000_000, '1/100', free_demand=10**12))

# A small-count exploration, not a useful paper-scale probability certificate.
def controller(history):
    return Intervention(q='3/4' if history and not history[-1]['success'] else '1/4')

path = HistoryMission(StochasticCycle(reactor)).run(
    (0,0,1000,0,0,0), 1000, 3, controller, seed=7, event_budget=250000)
```

`DensityExperiment` is a distinct mean-pulse mass-action illustration. It uses q*ell times the actual previous density, continuous food doses without integer rounding, and full six-dimensional ODE evolution, with separate integration intervals for recovery and collection. Its configured schedule is independent of the stochastic diagnostic controller. Radau and BDF agree to about 8e-10 in the saved example; the template-inventory residual is below 2e-15. These are numerical consistency checks, not probability estimates.

The saved default stochastic path completes all three small-count cycles and fails the joint success test. Failed outcomes are retained. Results include its entire cycle history and inventory ledger, one sampled paper-scale pulse, all channel definitions, density trajectories and cycle CSVs, phase weights, the complete error-budget sweep, certificate/sizing/accounting JSON, figures and reproducibility hashes.

Seven scientific test groups check symbolic source/phase identities, the exact normalized molecular pulse including catalyst-loss mass, event-level inventory and collection boundaries, failed-history continuation and budget semantics, the full sum against an independently nested error formula, integer sizing and synthesis, and independent density solvers plus a completed small-count trace. The mathematical certificate remains conditional on the paper's source, intervention and initial-state assumptions. Changing or adding a reaction requires new theorem justification. No Lean build is rerun; finite tests do not replace the formal probability, nonexplosion or process-identification proofs.

MIT is proposed for the example code, pending owner approval; this README grants no license and does not change manuscript rights.