giorgio/MachineSex
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Layer 1 core: Wright-Fisher knowledge-transmission model with E1-E2

Browse source 
Scaffold plus the Layer 1 analytical core and the first two experiments.

- knowledge/: truth, metrics, teachers (2.7.1 shared-switch construction),
  step, lineage, experiment, config, seeding (imported as `knowledge`).
- Validation spine green: neutral decay (Pred 1), fixation (Pred 2), exact
  mutation-drift equilibrium (Pred 3), union coverage (Pred 5). 68 tests pass.
- E1 reproduces tail-first collapse. E2 delivers the headline: a grounding
  phase boundary g* << 1, with stationary H tracking the exact H_eq closed
  form (g=0.005 -> 68% of truth diversity; g=0.05 -> 96%).
- Reproducibility: uv venv from a hash-pinned uv.lock is the source of truth;
  every run writes results.parquet + resolved_config.yaml + manifest.json
  (lib versions, git commit, sha256). Figures and manifests tracked; the
  large regenerable parquet is gitignored.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
Giorgio Gilestro 2026-07-04 18:10:18 +02:00
commit a6eb9b7512
33 changed files with 4356 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

491
tests/test_scientific_validation.py Normal file
View file

@ -0,0 +1,491 @@
"""
Scientific-validation test suite for Layer 1 of the Lamarckian Society model.
WHAT THIS FILE IS
-----------------
This is the *spine of trust* for Layer 1. It encodes the four analytic targets of
blueprint section 2.4 as executable assertions. It is simultaneously a scientific
check (the simulator reproduces known population-genetics results) and a code check
(the implementation is correct). If any test here fails, the science is wrong, not
just the code -- do not trust any downstream Layer-1 figure until these pass.
The Layer-1 implementation in `src/knowledge/` is DONE, for validation purposes,
when the conformance tests in Part 5 pass against the real package.
HOW IT IS STRUCTURED
--------------------
Part 1 Analytic theory closed-form ground truth (no simulation)
Part 2 Reference implementation minimal, pinned dynamics used to exercise the theory
Part 3 Tolerances honest Monte-Carlo tolerances, with rationale
Part 4 SPINE tests reference-vs-theory; ALWAYS run; prove the targets
Part 5 CONFORMANCE tests package-vs-theory; SKIP until `knowledge.*` exists
The four analytic targets (blueprint 2.4 / 2.7.1):
1. Neutral heterozygosity decay E[H_t] = H_0 (1 - 1/n)^t
2. Fixation probability P(item i fixes) = p_0^i
3. Mutation-drift equilibrium EXACT stationary H of the immigration model
4. Recombination union coverage U(K_T, rho, q) = T[ rho q + (1-rho)(1-(1-q)^K_T) ]
Every stochastic test is seeded and deterministic. Run with: pytest -q
"""
from __future__ import annotations
import numpy as np
import pytest
# =====================================================================================
# PART 1 -- ANALYTIC THEORY (closed forms; the ground truth these tests defend)
# =====================================================================================
#
# These functions contain NO simulation. They are the right-hand sides of the
# blueprint's analytic predictions. They are what everything else is compared against.
def theory_heterozygosity_decay(H0: float, n: int, t: np.ndarray | int) -> np.ndarray:
"""Blueprint 2.4-1. Expected heterozygosity under neutral Wright-Fisher resampling
of `n` items: E[H_t] = H0 (1 - 1/n)^t.
Derivation (exact): with p_{t+1}^i = c_i/n, c ~ Multinomial(n, p_t),
E[sum_i c_i^2] = n + n(n-1) G_t where G_t = sum_i (p_t^i)^2,
so E[G_{t+1}] = 1/n + (1-1/n) G_t and hence E[H_{t+1}] = (1-1/n) E[H_t].
"""
t = np.asarray(t, dtype=float)
return H0 * (1.0 - 1.0 / n) ** t
def theory_fixation_probability(p0: np.ndarray) -> np.ndarray:
"""Blueprint 2.4-2. Under neutral drift the probability that item i is the one
eventually fixed equals its initial frequency. So the target *is* p0."""
return np.asarray(p0, dtype=float)
def theory_mutation_drift_H_eq(n: int, m: int, H_star: float) -> float:
"""Blueprint 2.4-3, EXACT form for the grounding model actually implemented:
p_{t+1} = ( Multinomial(n, p_t) + Multinomial(m, p*) ) / (n + m).
Derivation. Let G_t = sum(p_t^2), G* = sum(p*^2), M_t = sum(p_t p*). Exact
multinomial moments give the linear mean recursions
E[M_{t+1}] = (n/N) E[M_t] + (m/N) G* => E[M_inf] = G*
E[G_{t+1}] = [ n + n(n-1)E[G_t] + m + m(m-1)G* + 2 n m E[M_t] ] / N^2 , N=n+m
Substituting M_inf = G* and solving the G fixed point, the sum(p*^2) terms cancel
against a factor of (1 - G*) = H*, leaving the clean closed form
H_eq = H* * m (2n + m - 1) / (n + 2 n m + m^2).
Limits: m->0 gives H_eq->0 (collapse to fixation); m->inf gives H_eq->H*
(the truth's own heterozygosity is recovered). In the rare-immigrant / many-types
limit (H*~1, m<<n) this reduces to theta/(1+theta) with theta = 2m, which is the
textbook infinite-alleles approximation the blueprint quotes. We test the EXACT
form, not the approximation, because a test should assert the strongest true thing.
"""
return H_star * m * (2 * n + m - 1) / (n + 2 * n * m + m * m)
def theory_union_coverage_fraction(K_T: int, rho: float, q: float) -> float:
"""Blueprint 2.4-5 / 2.7.1. Expected fraction of tail items retained by at least
one of K_T teachers built by the shared-switch construction (marginal retention q,
exact pairwise retention-correlation rho):
U/T = rho*q + (1 - rho) * (1 - (1 - q)^K_T).
Limits: K_T=1 -> q (single teacher, independent of rho); rho=1 -> q (identical
teachers, union = one); rho=0 -> 1-(1-q)^K_T (independent teachers, maximal union).
"""
return rho * q + (1.0 - rho) * (1.0 - (1.0 - q) ** K_T)
# =====================================================================================
# PART 2 -- REFERENCE IMPLEMENTATION (minimal, pinned dynamics)
# =====================================================================================
#
# This is the smallest correct implementation of the core Layer-1 dynamics. It exists
# so the spine tests can run before src/knowledge/ is written, and so the exact
# semantics the package must reproduce are unambiguous. The package will do far more
# (config, logging, regions, selection, re-minting, per-region metrics); it must agree
# with THIS on the analytic-check subset.
def ref_heterozygosity(p: np.ndarray) -> float:
"""Expected heterozygosity H = 1 - sum_i p_i^2. (Simpson diversity.)"""
p = np.asarray(p, dtype=float)
return float(1.0 - np.sum(p * p))
def ref_neutral_step(p: np.ndarray, n: int, rng: np.random.Generator) -> np.ndarray:
"""One neutral Wright-Fisher / Shumailov resampling step: draw n, refit."""
c = rng.multinomial(n, p)
return c / float(n)
def ref_grounded_step(p: np.ndarray, p_star: np.ndarray, n: int, m: int,
rng: np.random.Generator) -> np.ndarray:
"""One grounded step (immigration): pool n inherited draws with m real draws.
Denominator is exactly n+m since each multinomial's counts sum to its size."""
c = rng.multinomial(n, p) + rng.multinomial(m, p_star)
return c / c.sum()
def ref_make_retention_matrix(T: int, K_T: int, rho: float, q: float,
rng: np.random.Generator) -> np.ndarray:
"""Shared-switch exchangeable-Bernoulli construction (blueprint 2.7.1).
Returns an (K_T, T) 0/1 matrix with marginal retention q and EXACT pairwise
column-correlation rho. For each tail item j: shared switch z_j~Bern(rho),
shared retention s_j~Bern(q), independent u^k_j~Bern(q); r^k_j = s_j if z_j else u^k_j.
"""
z = rng.random(T) < rho # (T,) shared switch per item
s = rng.random(T) < q # (T,) shared retention per item
u = rng.random((K_T, T)) < q # (K_T, T) independent retentions
return np.where(z[None, :], s[None, :], u).astype(np.int8)
# ---- small helpers -----------------------------------------------------------------
def _zipf_p_star(K: int, s: float = 1.1) -> np.ndarray:
"""A true distribution with a genuine heavy tail (Zipf), normalised."""
w = 1.0 / np.arange(1, K + 1, dtype=float) ** s
return w / w.sum()
def _mean_heterozygosity_decay(n: int, K: int, T: int, reps: int,
seed: int) -> tuple[np.ndarray, float]:
"""Run `reps` neutral lineages from the uniform distribution; return the
replicate-mean heterozygosity per generation (length T+1) and H0."""
rng = np.random.default_rng(seed)
p0 = np.full(K, 1.0 / K)
H0 = ref_heterozygosity(p0)
Hbar = np.zeros(T + 1)
for _ in range(reps):
p = p0.copy()
Hbar[0] += ref_heterozygosity(p)
for t in range(1, T + 1):
p = ref_neutral_step(p, n, rng)
Hbar[t] += ref_heterozygosity(p)
return Hbar / reps, H0
def _empirical_fixation(p0: np.ndarray, n: int, reps: int, t_max: int,
seed: int) -> tuple[np.ndarray, int]:
"""Run `reps` neutral lineages to fixation; return empirical fixation frequency
per item and the number that fixed within t_max generations."""
rng = np.random.default_rng(seed)
K = len(p0)
counts = np.zeros(K)
n_fixed = 0
for _ in range(reps):
p = p0.copy()
for _t in range(t_max):
p = ref_neutral_step(p, n, rng)
nz = np.nonzero(p)[0]
if nz.size == 1:
counts[nz[0]] += 1
n_fixed += 1
break
else:
counts[int(np.argmax(p))] += 1 # not fixed in time (should be rare)
return counts / reps, n_fixed
def _stationary_heterozygosity(n: int, m: int, K: int, T: int, t_avg: int,
reps: int, seed: int) -> float:
"""Run `reps` grounded lineages; average H over the final t_avg generations and
over replicates -> an estimate of the stationary heterozygosity."""
rng = np.random.default_rng(seed)
p_star = _zipf_p_star(K)
vals = np.empty(reps)
for r in range(reps):
p = p_star.copy()
acc = 0.0
for t in range(T):
p = ref_grounded_step(p, p_star, n, m, rng)
if t >= T - t_avg:
acc += ref_heterozygosity(p)
vals[r] = acc / t_avg
return float(vals.mean())
# =====================================================================================
# PART 3 -- TOLERANCES (documented; Monte-Carlo error, not fudge factors)
# =====================================================================================
#
# Every tolerance below was calibrated: the empirical error at the given (reps, sizes)
# was measured, and the tolerance set a comfortable multiple above it, so the suite is
# robust to seed changes but still fails on a genuinely wrong implementation.
REL_TOL_DECAY = 0.02 # measured max rel err ~0.005 at reps=3000
ABS_TOL_FIX = 0.03 # measured max abs err ~0.007 at reps=4000
REL_TOL_HEQ = 0.02 # measured rel err <0.001 at the configured reps
ABS_TOL_UNION = 0.02 # measured max abs err ~0.004 at T=5e4
ABS_TOL_MARGINAL = 0.01 # retention marginal q
ABS_TOL_CORR = 0.03 # retention pairwise correlation rho
# =====================================================================================
# PART 4 -- SPINE TESTS (reference vs theory; ALWAYS run)
# =====================================================================================
class TestTheorySelfConsistency:
"""The closed forms must satisfy their own limits. Pure algebra; no simulation.
If these fail, the theory functions are miswritten and every other test is moot."""
def test_decay_at_t0_equals_H0(self):
assert theory_heterozygosity_decay(0.9, 100, 0) == pytest.approx(0.9)
def test_decay_is_monotone_nonincreasing(self):
H = theory_heterozygosity_decay(0.9, 50, np.arange(0, 100))
assert np.all(np.diff(H) <= 1e-15)
def test_Heq_zero_grounding_is_zero(self):
assert theory_mutation_drift_H_eq(100, 0, 0.9) == pytest.approx(0.0)
def test_Heq_infinite_grounding_recovers_truth(self):
# As m -> inf, H_eq -> H_star.
H_star = 0.9
big = theory_mutation_drift_H_eq(100, 10_000_000, H_star)
assert big == pytest.approx(H_star, abs=1e-3)
def test_Heq_rare_immigrant_limit_matches_textbook(self):
# H*~1, m<<n => H_eq ~ theta/(1+theta), theta = 2m.
n, m = 100_000, 3
H_star = 1.0
exact = theory_mutation_drift_H_eq(n, m, H_star)
theta = 2 * m
assert exact == pytest.approx(theta / (1 + theta), rel=1e-3)
@pytest.mark.parametrize("K_T", [1, 2, 3, 5])
def test_union_single_teacher_is_q(self, K_T):
# K_T=1 must give q for every rho.
for rho in (0.0, 0.5, 1.0):
assert theory_union_coverage_fraction(1, rho, 0.3) == pytest.approx(0.3)
def test_union_identical_teachers_no_benefit(self):
# rho=1 gives q regardless of K_T.
for K_T in (1, 2, 3, 5):
assert theory_union_coverage_fraction(K_T, 1.0, 0.3) == pytest.approx(0.3)
def test_union_independent_teachers_maximal(self):
# rho=0 gives 1-(1-q)^K_T.
q, K_T = 0.3, 5
assert theory_union_coverage_fraction(K_T, 0.0, q) == pytest.approx(
1 - (1 - q) ** K_T)
class TestHeterozygosityDecay:
"""Target 1: neutral drift decays heterozygosity geometrically at rate 1/n.
This is the quantitative form of 'collapse is tail-first and its rate is set by
the distillation sample size n'. Falsifier of the HARNESS (not the theory): if the
simulator's decay does not match, the simulator is wrong -- fix before proceeding."""
def test_decay_matches_geometric(self):
n, K, T, reps = 100, 50, 40, 3000
Hbar, H0 = _mean_heterozygosity_decay(n, K, T, reps, seed=101)
theory = theory_heterozygosity_decay(H0, n, np.arange(T + 1))
rel_err = np.abs(Hbar - theory) / theory
assert rel_err.max() < REL_TOL_DECAY, (
f"max rel err {rel_err.max():.4f} exceeds {REL_TOL_DECAY}")
def test_rate_scales_with_n(self):
# Larger n -> slower decay. Compare one-step drop for two n values.
K, reps = 50, 3000
H_small, H0 = _mean_heterozygosity_decay(50, K, 1, reps, seed=102)
H_large, _ = _mean_heterozygosity_decay(500, K, 1, reps, seed=103)
drop_small = H0 - H_small[1]
drop_large = H0 - H_large[1]
assert drop_small > drop_large > 0
class TestFixationProbability:
"""Target 2: under neutral drift, P(item i fixes) = p_0^i. A direct check that the
resampling has no hidden bias toward any item (which would silently distort collapse)."""
def test_fixation_equals_initial_frequency(self):
p0 = np.array([0.2, 0.3, 0.5])
freq, n_fixed = _empirical_fixation(p0, n=40, reps=4000, t_max=2000, seed=201)
assert n_fixed >= 0.99 * 4000, "lineages did not reach fixation within t_max"
assert np.max(np.abs(freq - theory_fixation_probability(p0))) < ABS_TOL_FIX
class TestMutationDriftEquilibrium:
"""Target 3: with grounding, heterozygosity reaches a positive stationary value
given EXACTLY by the immigration-model closed form. This is the phase boundary in
closed form -- the quantitative heart of 'a little grounding protects a lot of
inheritance'. Falsifier of the claim: if stationary H is flat in m, grounding buys
nothing (that scientific falsifier is exercised by experiment E2; here we validate
that the simulator hits the analytic curve)."""
@pytest.mark.parametrize("n,m,reps,T,t_avg", [
(50, 10, 150, 1200, 400),
(100, 20, 120, 2000, 600),
])
def test_stationary_H_matches_exact_formula(self, n, m, reps, T, t_avg):
K = 100
H_star = ref_heterozygosity(_zipf_p_star(K))
H_sim = _stationary_heterozygosity(n, m, K, T, t_avg, reps, seed=300 + n + m)
H_eq = theory_mutation_drift_H_eq(n, m, H_star)
assert H_sim == pytest.approx(H_eq, rel=REL_TOL_HEQ), (
f"n={n} m={m}: sim {H_sim:.4f} vs theory {H_eq:.4f}")
def test_grounding_raises_stationary_H_monotonically(self):
# More grounding -> higher stationary heterozygosity (closer to the truth).
K = 100
H_star = ref_heterozygosity(_zipf_p_star(K))
prev = -1.0
for m in (2, 10, 40):
H_eq = theory_mutation_drift_H_eq(100, m, H_star)
assert H_eq > prev
prev = H_eq
class TestRecombinationUnionCoverage:
"""Target 4: multi-teacher recombination. The shared-switch construction must hit
its target marginal q and pairwise correlation rho, and the resulting union tail
coverage must match the closed form. This is what makes 'collapse suppression is
proportional to decorrelation' an exact, checkable statement rather than a slogan."""
@pytest.mark.parametrize("rho", [0.0, 0.5, 1.0])
def test_marginal_retention_equals_q(self, rho):
rng = np.random.default_rng(400)
q, T = 0.2, 40_000
R = ref_make_retention_matrix(T, K_T=4, rho=rho, q=q, rng=rng)
assert R.mean() == pytest.approx(q, abs=ABS_TOL_MARGINAL)
@pytest.mark.parametrize("rho", [0.0, 0.25, 0.5, 0.75, 1.0])
def test_pairwise_correlation_equals_rho(self, rho):
rng = np.random.default_rng(401)
q, T, K_T = 0.3, 40_000, 4
R = ref_make_retention_matrix(T, K_T, rho, q, rng)
corrs = [np.corrcoef(R[a], R[b])[0, 1]
for a in range(K_T) for b in range(a + 1, K_T)]
assert np.mean(corrs) == pytest.approx(rho, abs=ABS_TOL_CORR)
@pytest.mark.parametrize("K_T", [1, 2, 3, 5])
@pytest.mark.parametrize("rho", [0.0, 0.5, 1.0])
def test_union_coverage_matches_closed_form(self, K_T, rho):
rng = np.random.default_rng(402)
q, T = 0.3, 50_000
R = ref_make_retention_matrix(T, K_T, rho, q, rng)
union_emp = (R.sum(axis=0) > 0).mean()
union_theory = theory_union_coverage_fraction(K_T, rho, q)
assert union_emp == pytest.approx(union_theory, abs=ABS_TOL_UNION)
class TestMetricSanity:
"""Cheap guards on the diversity metric itself; a broken metric invalidates every
curve above."""
def test_heterozygosity_bounds_and_extremes(self):
assert ref_heterozygosity(np.array([1.0, 0.0, 0.0])) == pytest.approx(0.0)
p = np.full(10, 0.1)
assert ref_heterozygosity(p) == pytest.approx(1 - 1 / 10)
assert 0.0 <= ref_heterozygosity(_zipf_p_star(100)) <= 1.0
# =====================================================================================
# PART 5 -- CONFORMANCE TESTS (package vs theory; SKIP until knowledge.* is built)
# =====================================================================================
#
# These are the acceptance gate for the real Layer-1 implementation. They import the
# normative interfaces of blueprint 2.7 and assert the package reproduces the same
# analytic targets as the reference above. They SKIP cleanly until the package exists,
# then must PASS. Do not weaken the assertions; if the package's config object differs
# from the minimal one built here, adapt the *construction* of cfg, never the tolerance.
try:
import knowledge.metrics as knowledge_metrics
import knowledge.teachers as knowledge_teachers
HAVE_KNOWLEDGE = True
except Exception: # package not built yet -> conformance layer skips, spine still runs
knowledge_metrics = knowledge_teachers = None
HAVE_KNOWLEDGE = False
requires_package = pytest.mark.skipif(
not HAVE_KNOWLEDGE, reason="Layer-1 package (knowledge.*) not implemented yet")
@requires_package
class TestPackageMetricsConform:
def test_package_heterozygosity_matches_reference(self):
for p in (np.array([1.0, 0.0, 0.0]), np.full(10, 0.1), _zipf_p_star(100)):
assert knowledge_metrics.heterozygosity(p) == pytest.approx(
ref_heterozygosity(p), abs=1e-12)
@requires_package
class TestPackageRetentionConform:
@pytest.mark.parametrize("rho", [0.0, 0.5, 1.0])
def test_package_retention_marginal_and_correlation(self, rho):
rng = np.random.default_rng(500)
q, T, K_T = 0.3, 40_000, 4
R = np.asarray(knowledge_teachers.make_retention_matrix(T, K_T, rho, q, rng))
assert R.shape == (K_T, T)
assert R.mean() == pytest.approx(q, abs=ABS_TOL_MARGINAL)
corrs = [np.corrcoef(R[a], R[b])[0, 1]
for a in range(K_T) for b in range(a + 1, K_T)]
assert np.mean(corrs) == pytest.approx(rho, abs=ABS_TOL_CORR)
@pytest.mark.parametrize("K_T", [1, 2, 3, 5])
def test_package_union_matches_closed_form(self, K_T):
rng = np.random.default_rng(501)
q, T, rho = 0.3, 50_000, 0.0
R = np.asarray(knowledge_teachers.make_retention_matrix(T, K_T, rho, q, rng))
union_emp = (R.sum(axis=0) > 0).mean()
assert union_emp == pytest.approx(
theory_union_coverage_fraction(K_T, rho, q), abs=ABS_TOL_UNION)
@requires_package
class TestPackageDynamicsConform:
"""Validate the package's generational dynamics via run_lineage. The cfg below is
the minimal contract run_lineage must honour (a mapping matching blueprint 2.7's
schema). It must return a tidy per-generation frame with a 'heterozygosity' column.
If the package uses a typed config object instead of a mapping, wrap the dict here;
do not change what is asserted."""
def _run_lineage(self, cfg_overrides, seed):
lineage = pytest.importorskip(
"knowledge.lineage", reason="Layer-1 package not implemented yet")
base = {
"truth": {"K": 50, "R": 1, "tail": "zipf", "zipf_s": 1.1,
"tail_frac": 0.5, "tail_threshold": 1e-3},
"dynamics": {
"n": 100,
"teachers": {"K_T": 1, "rho": 0.0, "q": 1.0},
"grounding": {"m": 0, "policy": "uniform"},
"selection": {"mode": "none", "novelty_alpha": 0.0},
"remint": {"enabled": False, "period": None, "H_gate": None},
},
"generations": 40,
"metrics": {"kl_floor": 1e-9},
}
# shallow-merge overrides
for k, v in cfg_overrides.items():
if isinstance(v, dict):
base[k] = {**base.get(k, {}), **v}
else:
base[k] = v
return lineage.run_lineage(base, seed)
def test_neutral_decay_via_package(self):
# Average H_t over replicate seeds and compare to the geometric law.
n, K, T, reps = 100, 50, 40, 200
H0 = 1 - 1 / K
Hsum = np.zeros(T + 1)
for s in range(reps):
df = self._run_lineage(
{"truth": {"K": K}, "dynamics": {"n": n}, "generations": T}, seed=s)
Hsum += df["heterozygosity"].to_numpy()[: T + 1]
Hbar = Hsum / reps
theory = theory_heterozygosity_decay(H0, n, np.arange(T + 1))
rel_err = np.abs(Hbar - theory) / theory
assert rel_err.max() < 0.05 # looser: fewer reps than the reference spine test
def test_grounded_equilibrium_via_package(self):
n, m, K, T, t_avg, reps = 50, 10, 100, 1200, 400, 60
H_star = ref_heterozygosity(_zipf_p_star(K))
vals = []
for s in range(reps):
df = self._run_lineage(
{"truth": {"K": K}, "dynamics": {"n": n, "grounding": {"m": m}},
"generations": T}, seed=s)
H = df["heterozygosity"].to_numpy()
vals.append(H[-t_avg:].mean())
H_sim = float(np.mean(vals))
H_eq = theory_mutation_drift_H_eq(n, m, H_star)
assert H_sim == pytest.approx(H_eq, rel=0.05)