""" 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/inheritance/` 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< 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/inheritance/ 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< 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 inheritance.metrics as knowledge_metrics import inheritance.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( "inheritance.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)