# POPCORN — certified population-genetics likelihoods

Tool page: https://bootloops.ai/tools/popcorn.html

KIND: package (`import popcorn`; CLI `certsfs.py`; battery `selftest.py`;
flat layout — the engine modules ship beside the package files and are
aliased in place by identity).

PURPOSE: The Poisson Random Field selection-SFS/DFE stack of the
polyDFE/fitdadi/fastDFE class made exact and certified:

- **Selection SFS, exact**: M_i = 1F1(n-i; n; -S) for the whole vector via ONE
  exact-rational tridiagonal solve — M_i = alpha_i + beta_i e^{-S} with alpha,
  beta exact Fractions (n <= 2000 practical); exact gradient stack, one extra
  solve per order, same matrix. Certified-precision assembly: the Fraction
  heights BOUND the cancellation a priori, so the assembly dps is computed,
  not guessed.
- **Selection SFS, biobank-n**: stable-direction contiguity recursion in arb
  ball arithmetic (python-flint) — ball-certified vectors + gradients to
  n = 10^5; radii track the directional contamination exactly, and precision
  doubles until the target is met.
- **Certified DFE mixing** (`popcorn.dfe.CertifiedKernel`): polyDFE model-C
  class mixtures with a certified kernel cache — likelihood evaluation is
  matrix-vector fast while every number carries a proof budget (GL degree-pair
  self-check, analytic S->0 and S->infinity tail bounds, exact analytic
  parameter gradients).
- **Dominance h != 1/2** (`popcorn.dominance`): E(S,h) is an entire
  q-deformation of the Kummer line. A certified 2-fold nested-quadrature
  oracle (degree-pair self-checks at both levels, h=1/2 collapse gate) plus a
  stable q-series evaluator, cross-checked against pinned reference values.
- **Exact certificate family** (`popcorn.certificates`): positivity of
  rational polynomials on [0,1] (two independent routes), pure-rational LP
  membership with Farkas witnesses, float-propose/exact-decide LP, exact
  integer-cone dual simplex, certify-after-float hull membership with
  certificate reuse (CertHull), and exact-rational REGION certificates
  (Bernstein box positivity, certified sup over a chart, Cauchy-Schwarz
  chi^2 lower bounds, rational reconstruction + Sturm root counting).
- **Exact Lambda-coalescent machinery** (`popcorn.lambda_coalescent`):
  merger rates lambda_{b,k} for Kingman, Beta(2-alpha, alpha) and Dirac(psi)
  coalescents (and the Kingman+Dirac two-atom mixture) as exact rationals, and
  the exact expected branch-length spectrum E[L_i], i = 1..n-1 (hence the
  expected unfolded SFS) by the standard first-transition recursion, all in
  Fractions with built-in exact identity checks (Kingman 2/i closed form,
  alpha=2 -> Kingman, alpha=1 -> Bolthausen-Sznitman, psi=0 -> Kingman,
  psi=1 -> star); ships an n = 20 exact reference grid (Kingman, 7 Beta,
  6 Dirac) and an exact linear functional separating that grid from the
  variable-population-size Kingman class.
- **Two-locus branch-length moments** (`popcorn.twolocus`): the
  second moments E[T_i^A T_j^B] (i, j = 1..n-1) of the branch lengths
  subtending i samples at locus A and j at locus B, for a sample of n under
  the coalescent with recombination (scaled rate rho = 4 N_ref R, each doubly
  ancestral lineage splitting at rho/2) and piecewise-constant N(t); i.e. the
  expected joint two-locus frequency spectrum up to scale, plus the Kingman
  first moments. Deterministic: the labeled two-locus configuration chain is
  enumerated once per n, the terminal epoch is one sparse LU solve and each
  finite epoch one expm_multiply — no simulation, no time stepping.
- **Transient selected SFS, large samples** (`popcorn.transient`): the
  expected unfolded SFS E_i (i = 1..n-1, theta = 1) under genic selection S
  after a piecewise-constant size history [(nu, T), ...] following an
  ancestral Wright equilibrium, by direct integration of the Kimura forward
  diffusion in u = x(1-x)f — Scharfetter-Gummel exponentially fitted fluxes
  on a fixed log/lin/log grid, Crank-Nicolson with Rannacher startup, banded
  LAPACK solves with one factorization per epoch, binomial projection to
  sample size n (n ~ 1000-2500 design range) — plus hypergeometric down-sampling
  n -> m (formula-exact, float) and folding. Not a moment closure: n enters only
  through the projection. FLOAT-VALIDATED register (measured trust radius,
  see FOOTGUNS), never certified.
- **Ancestral-state join** (`popcorn.ancestral`, engine `epo_join.py`):
  stream a sites table (chrom pos ref alt, or pos ref alt) against an
  Ensembl-EPO-style ancestral-allele FASTA, classify each site by EPO
  confidence (uppercase = high, lowercase = low, `N`/`-`/`.` = no call),
  orient REF/ALT to ancestral/derived over biallelic SNVs, and write the
  polarized `pos ref alt anc conf` table with counters (confidence
  histogram, anc==ref / anc==alt / third-allele mismatch by confidence,
  unpolarizable) and the summary rates. Pure Python, standard library;
  measured about 0.6 million sites per second on a 5 Mb synthetic ancestor.

USE-WHEN:

- Certified selection-SFS entries/vectors/exact gradients for PRF DFE
  likelihoods; auditing a float-pipeline popgen fit; biobank-scale n where
  float pipelines are structurally unavailable.
- Certified DFE-mixing kernel; dominance h != 1/2.
- Exact hull/positivity/LP certificates, and region certificates — "the data
  point is provably >= this far from the WHOLE model region", not just from a
  grid of samples.
- Exact neutral expected SFS under any Lambda-coalescent for n up to a couple
  of hundred (constant population size); an exact reference against which to
  validate a float or Monte Carlo implementation (e.g. msprime's Beta and
  Dirac models); exact sign decisions on linear functionals of the spectrum.
- Two-locus: you need the recombining two-site spectrum (or its rho = 0 /
  rho -> infinity limits) at small n as a smooth function of rho and of a
  piecewise-constant size history — as the model side of a two-site
  composite likelihood, to calibrate LD-decay summaries, or as the
  deterministic comparand for a coalescent simulator. n <= 7 is sub-second;
  n = 8 at rho > 0 costs ~13 s and ~0.8 GB per evaluation (states grow ~2.9x
  per added sample).
- Non-equilibrium selected spectra at sample sizes where moment closures
  degrade (n ~ 10^3), e.g. a DFE likelihood under a bottleneck/growth
  history; stationary questions go to the certified engine (`popcorn.sfs`,
  `certsfs.py`) instead.
- Polarizing a VCF-derived sites list to ancestral/derived before building
  an unfolded SFS or any derived-allele statistic; measuring the
  ancestral-mismatch and low-confidence rates of a call set against the EPO
  ancestor; producing the `pos ref alt anc conf` table other tools consume.

NOT-FOR: certified demography or demographic inference. The certified
SFS/DFE/dominance stack is equilibrium, constant-N only; piecewise-constant
N(t) is available only at FLOAT-VALIDATED register (`popcorn.transient`:
single-locus selected SFS at large n; `popcorn.twolocus`: neutral two-locus
branch-length moments), and `popcorn.lambda_coalescent` is neutral and
constant-N. No migration or multi-population spectra, no continuous size
change, no fitting driver. The 1F1 closed form for
the fixed-S SFS is folklore (Zivkovic et al. 2015), not claimed here. The
package's built-in self-checks certify internal consistency; a production
claim should additionally be gated against an oracle that shares no code with
the leg being checked.

INVOKE:

```sh
python3 tools/popcorn/certsfs.py selftest --dps 40      # pins + reference values
python3 tools/popcorn/certsfs.py entry  20 3 -100 --dps 40
python3 tools/popcorn/certsfs.py check  100 50 -1000 --dps 40   # two routes + agreement
python3 tools/popcorn/certsfs.py vector 100 -1000 --dps 40 --grad
python3 tools/popcorn/lambda_exact.py                       # exact identity checks + n=20 reference grid (prints only)
python3 tools/popcorn/twolocus_engine.py 6 1.0               # two-locus demo: states, E[T_i], E[T_i^A T_j^B] at n=6, rho=1
python3 tools/popcorn/transient_sfs_engine.py                # transient engine smoke (~1 s)
python3 tools/popcorn/epo_join.py --fasta ancestor_21.fa --sites sites.tsv --out anc_chr21.tsv.gz --chrom chr21 --summary counts.json
```

```python
import popcorn                       # with tools/ on PYTHONPATH
popcorn.verify()                     # sha pins; [] iff clean — treat any hit as fatal
alpha, beta = popcorn.sfs.solve_M_exact(50, -10)         # exact Fractions
K = popcorn.dfe.CertifiedKernel(20, dps=40); K.build(-1); K.build(1)
Ev, selfcons_d, tail_rel = K.mix({'b': 0.4, 'Sd': -1000, 'pb': 0.02, 'Sb': 10})
E, sc = popcorn.dominance.E_dominant(20, 3, -100.0, 0.3, dps=50)
C = popcorn.certificates             # lazy: scipy only if cert_lp touched
status, cert = C.exact_lp.membership(points, q)          # exact Farkas witness
R = C.region                         # region certificates (see its docstring)
from fractions import Fraction
L = popcorn.lambda_coalescent        # exact, standard library
h = L.expected_lengths(50, L.beta_rate(Fraction(3, 2)))[50]   # {i: E[L_i]} exact Fractions
x = L.xi_hat(50, L.dirac_rate(Fraction(1, 10)))               # normalized expected SFS, sums to 1 exactly
W2 = popcorn.twolocus                # numpy + scipy; float-validated
M, a, b = W2.moments(6, 0.9)         # E[T_i^A T_j^B], E[T_i] at n=6, rho=0.9; W2.pool_pairs(M)
T = popcorn.transient                # numpy + scipy; float-validated
eng = T.TransientSFSEngine(n=1000)   # build once per n
E, meta = eng.expected_sfs(-5.0, [(0.33, 0.47), (3.38, 0.017)]); assert meta['negative_entries'] == 0
from popcorn.ancestral import join_sites, polarization_rates
counts = join_sites('ancestor_21.fa', 'sites.tsv', 'anc_chr21.tsv.gz', chrom='chr21'); rates = polarization_rates(counts)
```

INPUTS: S = 4*Ne*s (S>0 advantageous), unfolded derived-allele SFS, theta=1
normalization. polyDFE/fastDFE S == this S; dadi/fitdadi gamma == S/2.
Dominance conventions: fitnesses 1 : 1+2sh : 1+2s. entry/check/selftest need
mpmath only; vector needs python-flint above n=2000. Certificate polynomials:
dict {exponent: rational} on [0,1] (positivity) or {(es, eu): Fraction} on
[0,1]^2 (region). Lambda-coalescent parameters as Fractions or 'p/q' strings;
two-locus time in 2 N_ref generations, eta = N_ref/N, rho = 4 N_ref R;
transient epochs [(nu, T), ...] past -> present, genic S only; EPO FASTA one
chromosome per file, sites as `chrom pos ref alt` or `pos ref alt`.

OUTPUTS: CLI prints values at the requested dps (check mode: two routes +
agreement digits; below dps-5 agreement is a named CertSFSGateError, rc=3 —
never a silent fallback). Library returns exact Fractions / certified balls /
kernel objects / exact certificates that are independently re-substituted
before being returned.

REQUIREMENTS: python3 + mpmath (core). Optional, each unlocking legs that
otherwise SKIP BY NAME: python-flint (ball route, kernel, cone/positivity
route R), gmpy2 (LP certificate family), numpy + scipy (float proposers;
`popcorn.twolocus` and `popcorn.transient` — legs twolocus_* and transient_*
SKIP BY NAME without them), sympy (Sturm root counting, positivity route S).
`popcorn.lambda_coalescent` and `popcorn.ancestral` need the standard library
only.

GATES: every pinned file byte-identical (`popcorn.verify()`, and the CLI
re-verifies before any engine import); certsfs check agreement >= dps-5 or
rc=3; battery default tier = the fast legs below; `--full` adds the held-out
two-route gate (gate_phase1.py: >= 40 matched digits at held-out points to
n = 10007, gradient gate via the parameter-shift identity), the dominance
oracle's full selftest, and the n=20 kernel selftest.

Battery: `python3 tools/popcorn/selftest.py` → one line per leg + OVERALL
PASS, rc=0. Measured ~19 s with all optional engines present; legs missing an
engine SKIP BY NAME stating what to install (skips do not fail the battery).
Includes planted negative controls (exterior LP points, negative polynomials,
and the region selftest's MUST-FAIL plants). `--full` is minutes-scale. The
battery writes nothing into the package tree (scratch goes to a temp dir).

FOOTGUNS:

- The shipped engine bytes are sha-pinned; `popcorn.verify()` and the CLI
  fail closed on drift. Regenerate pins only when deliberately changing
  shipped bytes (`python3 -c "import popcorn._pins as p; p.regen()"`), and
  update `certsfs.PINNED_SHA256` in the same change.
- Exact assembly: alpha/beta are exact for a SPECIFIC rational S — assembling
  with any other nearby S (e.g. the binary float of a decimal) is amplified
  by the full cancellation factor (measured: thousands of digits lost).
  Thread ONE exact rational S through solve and assembly.
- Never substitute mp.quad for certquad: mp.quad's error estimate is absolute
  and it silently underconverges on exponential-ramp panels — the defect
  certquad exists to cure. A wrong lnf estimator costs panels, not
  correctness, BUT an lnf that reports the floor over live support silently
  drops mass — sanity-check lnf against f.
- Convert parameters to mpf INSIDE the workdps block: converting outside
  rounds them at ambient precision (measured: a uniform 16-digit ceiling).
- Degree-pair self-consistency between correlated legs is not proven digits;
  the h=1/2 collapse gate and the two-route CLI check are the independent
  anchors shipped here.
- `bern_nonneg` (region) is subdivision-only: a polynomial that TOUCHES zero
  on the box needs touch-root deflation — positivity route B.
- `cone_lp_b`'s one-RREF bit-sorted basis completion is load-bearing: the
  greedy per-column variant of the same simplex measured ~40x slower. Do not
  "simplify" it back.
- q-series dominance points at extreme h and large |S| cost the full computed
  precision pad (deliberate; can run minutes) — prefer the oracle there.
- Consumers fitting likelihoods built on these outputs: optimize on the
  deviance scale C*(d - log1p(d)), not raw logL — when the objective's
  dynamic range swamps the curvature, L-BFGS-B/FD-Hessian breaks on raw logL
  and the deviance form recovered ~1e-11 relative accuracy in testing.
- Register EXACT but pure-Python Fraction arithmetic: cost is O(n^3) rational
  operations whose operand size grows with n and with the family (Beta(3/2):
  ~2 s at n=100, ~17 s at n=200; Dirac(1/2): ~3 min at n=200; Kingman is
  O(n^2)). Pass parameters as Fractions or strings ('3/2'), never as binary
  floats, or the rates are exact for the wrong rational. Constant population
  size and expected values only. The shipped n=20 Kingman-class functional's
  nonnegativity over the variable-size Kingman class is carried by the
  reference file as given, not re-derived by this module; the module only
  evaluates it exactly.
- Two-locus register is FLOAT-VALIDATED, not certified — double-precision
  sparse linear algebra with no enclosure; validated at rho = 0 against the
  exact-rational Kingman E[L_i L_j] (n = 8: ~2e-16 on the normalized pair
  vector, ~1e-13 relative on raw entries), at n = 2 against the closed-form
  two-locus covariance (rho + 18)/(rho^2 + 13 rho + 18) for all rho, and by
  self-consistency across rho (Kingman marginals 2/i at any rho, A/B
  symmetry, epoch-splitting invariance, time-rescaling covariance, ~1e-15).
  Never mix its output into an exact/certified table without the label.
  Units: time in 2 N_ref generations, eta = N_ref/N, rho = 4 N_ref R;
  `epochs_from_Nt` divides generations by 2 N_ref (default N_ref = 1e4).
  `pool_pairs` doubles the diagonal (it reads M + M^T on i <= j): it is a
  pooling convention, not the class law of an unordered pair of sites —
  compare only against vectors pooled the same way.
- `popcorn.transient` is FLOAT-VALIDATED, not certified: double precision
  with a measured trust radius and nothing more. Measured (data shipped
  under `reference/transient/`): at n = 1000 on the default grid, the
  projected equilibrium and the equilibrium held T = 0.5 through the
  integrator agree with the certified two-route stationary references
  (shipped to 20 significant digits) to <= 9.4e-5 max-relative in every
  frequency band for S in {0, -1, -5, -20, -100}, quadrature-dominated and
  concentrated in the log-spaced bands; transient self-convergence on a
  two-epoch bottleneck-then-growth history, entries E_i > 1e-5: dt
  refinement 4e-4 -> 1e-4 <= 3.9e-4 for |S| <= 20, 4.0e-3 at S = -100,
  1.1e-2 at S = -200 (grid refinement <= 5.2e-3 throughout); the
  exponentially suppressed tail (E_i < 1e-5) moves by up to 9e-2 at
  S = -200. Other n, |S| > 200, custom grids that break seam spacing
  continuity (an abrupt jump costs about 2.5e-4), larger dt0 and sharper
  histories are unmeasured. `expected_sfs` never raises on negative tail
  entries — check `meta["negative_entries"]`. Selection is genic (h = 1/2)
  only; for dominance use `popcorn.dominance` (stationary). Never quote a
  number from this module beside certified ones without labeling each.
- Register is DATA-UTILITY: integer counters and string columns, no
  estimator and nothing certified; the summary rates are plain ratios of the
  counters. One chromosome per FASTA: only the first record is read (further
  records are ignored and flagged as `fasta_extra_sequences`), so always
  pass `--chrom` with a multi-chromosome sites table or rows of other
  chromosomes are joined against the wrong sequence. Positions are 1-based;
  a position past the sequence end yields anc `.` and counts in
  `out_of_fasta_range`, it is not an error. Orientation case-folds the
  ancestral base (a low-confidence `t` orients like `T`); filter on `conf`
  downstream if only high-confidence calls are wanted. Compare `.gz`
  outputs by decompressed payload, never by file hash (gzip stores an
  mtime).

CREDIT: Implements-and-certifies the Poisson Random Field selection-SFS/DFE
methodology of the Sawyer & Hartl (1992) line as used by polyDFE (Tataru et
al. 2017), dadi (Gutenkunst et al. 2009) / fitdadi (Kim et al. 2017), and
fastDFE; the fixed-S 1F1 closed form is due to the literature (Zivkovic et
al. 2015). The certificate instruments are classical exact methods (Bernstein
/ de Casteljau subdivision, Farkas certificates, Sturm's theorem) arranged
under a float-proposes/exact-decides discipline.

Lambda-coalescents: Pitman (1999), Sagitov (1999); Beta(2-alpha, alpha)
family: Schweinsberg (2003); psi-coalescent: Eldon & Wakeley (2006); expected
SFS recursions under Lambda-coalescents: Birkner, Blath & Eldon (2013);
Kingman closed form E[L_i] = 2/i: Fu (1995); Bolthausen & Sznitman (1998).

Two-locus: the two-locus ancestral process with recombination is due
to Griffiths (1981) and Hudson (1983); the n = 2 covariance formula is
theirs (see McVean 2002, Genetics 162:987); the rho = 0 comparand is the
Kingman branch-length covariance of Fu (1995, Theor. Popul. Biol. 48:172);
two-locus sampling theory under constant and variable N(t) is Hudson (2001,
Genetics 159:1805) and Kamm, Spence, Chan & Song (2016, Genetics 203:1381,
arXiv:1510.06017), of which this chain is the branch-length-moment
counterpart.

Transient: integrates the Kimura (1964) forward diffusion in the
Poisson-random-field setting (Sawyer & Hartl 1992) for the non-equilibrium
frequency spectrum (Evans, Shvets & Slatkin 2007, Theor. Popul. Biol.
71:109), using Scharfetter & Gummel (1969, IEEE Trans. Electron Devices
16:64) exponential fitting and Rannacher (1984, Numer. Math. 43:309)
startup. dadi (Gutenkunst et al. 2009) is the standard finite-difference
diffusion route and moments (Jouganous et al. 2017, Genetics 206:1549) the
moment-closure route for the same problem; this module takes the PDE route
so that the sample size enters only through the final projection.

Ancestral join: The ancestral calls and their upper/lower-case confidence convention
are those of the Ensembl EPO (Enredo-Pecan-Ortheus) pipeline: Paten et al.
(2008) Genome Res. 18:1814 (doi:10.1101/gr.076554.108) and 18:1829
(doi:10.1101/gr.076521.108); Herrero et al. (2016) Database bav096
(doi:10.1093/database/bav096); the same convention annotates the ancestral
allele in the 1000 Genomes Project releases (Nature 526:68, 2015,
doi:10.1038/nature15393). This module only joins and counts.

## License

MIT License. Copyright (c) 2026 Anthropic, PBC. Created by Matthew D. Schwartz; code written by Claude (Anthropic) under his supervision. Documentation is licensed CC BY 4.0. See LICENSE, LICENSE-CONTENT and NOTICE at the repository root; third-party components keep their own licenses (THIRD_PARTY.md).