# Prediction coverage with inherited cell states

This package implements the actual two-type branching process and its
mean-composition scalar closure. Both match the expected count and expected birth
and death fluxes; their count distributions can still give different prediction
coverage. The code provides general constant-rate models, founder/plating laws,
low-count probabilities, moment diagnostics, numerical corrected thresholds and
the paper's exact instance-specific certificate.

## Run and modify

Python 3.11 was used. In a virtual environment:

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

`--simulate` adds one optional molecule-free, integer cell-count Gillespie
realization. An event-budget interruption preserves the current population and
reports unfinished, never a terminal probability. No simulation is used to prove
the default coverage claims.

The top input block sets the two birth rates, two death rates, S->R and R->S
switching, resistant-founder probability, horizon, target coverage and founder
counts. All rates are per day. Defaults are the synthetic paper witness:
bS=0, bR=.1, dS=.3, dR=0, alpha=.001, beta=.00001, pR=5/24,
T=log(2)/.10001, about 6.93078 days. The isolated switching times of 1,000 and
100,000 days are deliberately extreme persistence, not fitted cell-state rates.
The positive-rate extension and fast-switching control are computed separately.

Changing inputs affects numerical predictions. It does not transfer the fixed
published witness's exact constants to another preparation, horizon, target level
or interval rule. In particular the certified endpoint 3 is not a universal
add-one correction.

## Model architecture

`Rates` specifies constant nonnegative per-cell rates. Division adds one cell of
the parent's current type; switching changes type while preserving total count.
`FounderPreparation` specifies a number of independent potential founders, their
independent plating probability, and one independent type draw per present founder.
Descendants inherit the current parental type; the founder type is not redrawn at
every division. Fixed numbers of each initial type are a different preparation.

`MomentClosure` integrates exact closed type means and their covariance matrix,
the accumulated-birth integral, scalar count variance, scalar mean independently,
and the integrated covariance causing the variance defect. For one randomly typed
founder the initial type covariance is p(1-p)*[[1,-1],[-1,1]], not zero; the initial
total-count variance is zero because that matrix sums to zero.

`ScalarTerminalLaw` constructs Kendall's one-founder law from total mean m and
accumulated-birth parameter A:

```
P(Z=0) = 1 - m/(1+A)
P(Z=n) = m/(1+A)^2 * (A/(1+A))^(n-1), n>=1.
```

`BackwardCountLaw` integrates the two low-degree generating-function systems.
Their degree-n derivative involves only degrees 0..n, so a cap on polynomial
degree is **not a cap on population size**. Families that grow beyond that count
and later lose descendants still contribute correctly. The only numerical error
in those low coefficients is the ODE calculation; unrepresented higher terminal
counts are not silently removed by a stopped population model.

`FounderPreparation.compose` forms the low coefficients of g(z)^F for fixed
independent founders, or (1-f+f*g(z))^M for binomially plated founders. Binary
polynomial powering truncates high degrees without renormalizing low ones.
Reported mass deficits are numerical diagnostics of the resulting terminal law.

`ForwardOverflowLaw` independently solves the joint forward equation. It sends
every exit beyond its population cap, and any excess initial count, to an
absorbing overflow state. In exact arithmetic, retained event mass is a lower
bound and adding overflow is an upper bound. The floating matrix exponential is
not validated interval arithmetic, so the package labels these checks as
diagnostics rather than rigorous numerical enclosures.

`PredictionAssessment` evaluates the scalar equal-tailed integer interval against
the actual law, then computes the actual one-sided upper endpoint separately.
If a coefficient cap cannot resolve a requested quantile it raises rather than
returning a made-up endpoint. Numerical endpoints for new parameters are
recommendations under the specified model, not certified coverage guarantees.

`EventSimulator` is a separate exact-event count model, useful for extensions and
trajectory illustrations. Event rates depend on the actual numbers of each type.
It shares the physical preparation and rates, not the scalar mean closure.

For example:

```python
from example import Rates, scenario
rates = Rates(birth=(.01,.1), death=(.3,.01), switching=(.1,.1))
moments, one_founder, scalar, assessments = scenario(
    rates, p=5/24, duration=6.9307787277, founders=(1,100), cap=192)
```

For time-dependent rates, the backward equation must be formulated in starting
time with a fixed terminal time. The elapsed-time equation implemented here is
specific to constant rates; passing varying rate functions is not supported.

## Exact witness and the repair

`exact_witness()` freshly expands the positive polynomial certificate bounding
the birth integral, checks its twelve positive **unnormalized** Bernstein
coefficients, integrates the cubic majorant exactly, verifies killed-history
ODE identities, and computes the rational tail, repair and variance bounds.
These finite checks use the paper's conventional process arguments to connect
their algebra to probability; they do not rerun Lean.

The analytic mean/birth box is m<=11/20, A<=2/5. It gives scalar coverage of
{0,1,2} at least 332/343. The resistant killed process keeps paths with no switch,
with q_n(T)=1/2*(c/2)^(n-1), c=10000/10001. Its seven live count states suffice to
extract the tail contributions at counts 3..7:

```
L = 151415045117539070312500 / 3001800450060004500180003 > .05.
P(actual count <= 2) <= 1-L < .95.
```

Selected sensitive no-event/death histories and resistant no-switch counts 1..3
give coverage at endpoint 3 at least
467880212635/481696324816 > .9713. Since coverage increases with the endpoint,
3 is the smallest successful endpoint **among {0,...,k}** for this preparation.
The full no-switch geometric tail gives the stronger conventional upper bound
about .9479323 at endpoint 2. A mean-only Markov bound would require endpoint 10.

The numerical witness gives actual/scalar coverage .9475977178/.9767606151 at
endpoint 2 and actual coverage .9738271218 at endpoint 3. The scalar central
95% interval is numerically [0,2]. The analytic bound 332/343 alone does not prove
that equal-tailed quantile choice, because its complementary bound exceeds .025.
The code and outputs keep that numerical identification separate.

## Why more founders do not fix the error

The exact generator identities LN=G and LN^2=2NG+B+D imply

```
Var(actual)-Var(scalar)
 = 2*m(t)^2 * integral Cov(N, G-gbar*N)/m(t)^2 dt
```

for equal initial count variance. Switching preserves N but changes composition;
inheritance matters when larger families covary with greater composition-weighted
net growth. The moment solver integrates this identity independently and compares
it with the directly evolved variances. Equal net growth in all types makes this
defect zero, even when their turnover rates differ; equal variances need not mean
equal full laws. Inheritance does not universally inflate variance.

For the witness, actual/scalar variances are about 1.087594/.583256. The scalar
central interval's actual coverage approaches
2*Phi(z_.975*sqrt(v_scalar/v_actual))-1, numerically .848799756.
The paper's exact variance inequalities imply an upper bound below .9308 on that
limit. They do not bound every finite-founder point: F=2 and F=3 exceed .95,
while F=100 gives [38,67] with actual coverage about .851146. Discrete endpoints
produce jumps, so coverage need not be monotone in founder number.

The default sweep evaluates both switching rates and recomputes quantiles at each
point. The numerical extension with positive birth/death rates in both types and
both switching rates .1/day gives about .906285 coverage at 100 founders; the
exact witness constants do not apply to it.

## A control where the closure performs well

The code also implements the paper's source-linked Figure 8 two-type preparation
from Browning's pinned repository commit
`401c054dcf588262edfa1263767ffa7ce366533b`: alpha=1, beta=.01, bR=.1, dS=.3,
horizon 7 days, Binomial(1000,f) plated founders with
f=817*614/(pi*4500^2), and resistant fraction the smaller root of
.05*theta^2-.57*theta+.02=0. This control has a different founder law and is not a
point on the witness founder curve. Its interval [3,24] has actual coverage about
.9532165, independently cross-checked with a forward overflow model.

Both ratio conventions are explicit in the output. For this control,
v_actual/v_scalar is about 1.034 (the source script's convention), while
v_scalar/v_actual is its reciprocal (the convention in the founder-limit formula).
It would be incorrect to insert 1.034 directly into that formula. This reproduces
the paper's two-type preparation, not the source publication's full continuous-
phenotype inference or experimental fit.

## Outputs, evidence and scope

`results.json` contains exact witness bounds, numerical laws and moments, founder
assessments, switching sweep, positive-rate examples and source-linked control.
The four CSV files expose count laws, moment/covariance histories, founders and
switching results. PNG/SVG figures, console text and output hashes are included.
Seven test groups check exact certificates, witness values, independent forward
agreement, grow-then-shrink contributions, type-independent special cases, equal
net-growth variance, founder/plating composition, the control, budget interruptions
and invalid inputs. Saved probabilities are never renormalized to hide missing mass.

The paper has conventional proofs with machine-checked components. The unrestricted
process-to-mean and scalar-process-to-explicit-law bridges remain conventional;
the variance-defect and general low-count results are conventional too. This code
performs exact algebra and numerical cross-checks, not a new formal verification.
Prediction regions concern repeated future counts under known fixed parameters;
they are not confidence or posterior credible intervals. No density dependence,
shared environment, related founders, measurement error or parameter uncertainty
is included. No cell line, molecular mechanism or treatment is inferred.

Manuscript SHA-256:
`a30348c7b611dc86e0435ca9aa69eeaafaeb26bd6992488937114d9f71f78f1d`.
The full PDF, its check script and the local `browning_count_boundary.py` source
were read. Equations and constants are implemented locally; no saved result or
success flag is loaded. MIT is proposed for the new code pending the owner's
license choice; this README grants no license and changes no manuscript terms.