society: multi-locus recombination frame — the vertical claim (E7/E8)

Enter the Lamarckian society with a robust theoretical frame. The single-
locus, fixed-p* model can only express recovery toward a ceiling; the
society's load-bearing claim is vertical -- capability that EXCEEDS any
component. Generalize knowledge to a distribution over genotypes (L
biallelic loci, K=2^L, additive fitness = # correct loci), reusing all the
K-mode machinery. The one new operator is recombination: free recombination
sends p -> product of per-locus marginals (linkage equilibrium).

E8 (star, kind: society) -- the vertical claim / Fisher-Muller: decorrelated
PARENTS (specialists, expert on their loci, agnostic elsewhere) are
recombined; sexual merge assembles a genotype fitter than any parent,
climbing to the optimum (12/12, a genotype no parent had) as parent count
grows and rho->0, while the best single parent (~8.7) and the mean-mixture
"model soup" (~11.6) plateau below. Reuses make_retention_matrix (locus
mastery replaces tail-item retention).

E7 (kind: genotype_lineage) -- the advantage of sex: a single population
adapts toward the optimum; the sexual lineage adapts faster than asexual
(clonal interference) by keeping loci in linkage equilibrium (LD->0 vs LD
spike). Honest scope: a speed advantage, not a permanent Muller's-ratchet
gap (subtle to force); E8 carries the headline.

Metaphor shift (per GG): the society is sexual reproduction with UNBOUNDED
parents, not teacher->pupil. Teacher->pupil caps at the ceiling; n-parent
recombination is combinatorial and generative, and unlike biology there is
no two-parent limit. Collapse = asexual degradation; the cure = sex. This
unifies E4 (merge != average) + E6 (irreversibility) under evolution-of-sex
theory and reaches ground Riis's single-locus n-grams cannot.

New: knowledge/{genotype,genotype_lineage,society}.py, configs/layer1/{E7,
E8}.yaml, figures/plot_{E7,E8}.py, READMEs, tests/test_genotype.py (+7).
experiment.py dispatch (kind in {genotype_lineage, society}); make layer1
wired. 112 tests green.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

src/knowledge/experiment.py

@@ -127,6 +127,46 @@ def run_experiment(cfg: dict) -> pd.DataFrame:
     return out
 
 
+_GENOTYPE_KEYS = ("genotype", "generations")
+
+
+def run_genotype_experiment(cfg: dict) -> pd.DataFrame:
+    """Run a genotype lineage across a sweep x replicates (E7, advantage of sex).
+
+    Mirrors ``run_experiment`` (paired replicate seeds) but assembles the base from the
+    ``genotype``/``generations`` blocks and calls ``run_genotype_lineage``. Sweeps use the same
+    dotted-path ``_apply_param`` (e.g. ``genotype.recomb_rate`` for asexual vs sexual).
+    """
+    from .genotype_lineage import run_genotype_lineage
+
+    base = {k: copy.deepcopy(cfg[k]) for k in _GENOTYPE_KEYS if k in cfg}
+    sweeps = cfg.get("sweep", [])
+    if isinstance(sweeps, dict):
+        sweeps = [sweeps]
+    params = [s["param"] for s in sweeps]
+    value_lists = [list(s["values"]) for s in sweeps]
+    combos = [({}, base)] if not sweeps else []
+    for values in itertools.product(*value_lists):
+        lin = copy.deepcopy(base)
+        label: dict = {}
+        for param, val in zip(params, values):
+            label.update(_apply_param(lin, param, val))
+        combos.append((label, lin))
+
+    seeds = spawn_seeds(int(cfg["seed"]), int(cfg["n_replicates"]))
+    frames: list[pd.DataFrame] = []
+    for label, lin in combos:
+        for rep, ss in enumerate(seeds):
+            df = run_genotype_lineage(lin, int(ss.generate_state(1)[0]))
+            for col, val in label.items():
+                df[col] = val
+            df["replicate"] = rep
+            frames.append(df)
+    out = pd.concat(frames, ignore_index=True)
+    out.insert(0, "experiment", cfg["experiment"])
+    return out
+
+
 def run_coverage(cfg: dict) -> pd.DataFrame:
     """E4 runner: multi-teacher recombination coverage (blueprint 2.5-E4 / 2.7.1).
 
@@ -278,7 +318,16 @@ def run_and_save(config_path: str | Path) -> Path:
     config_path = Path(config_path)
     cfg = yaml.safe_load(config_path.read_text())
     out_dir = Path(cfg.get("output", {}).get("dir", f"results/{cfg['experiment']}"))
-    df = run_coverage(cfg) if cfg.get("kind") == "coverage" else run_experiment(cfg)
+    kind = cfg.get("kind", "lineage")
+    if kind == "coverage":
+        df = run_coverage(cfg)
+    elif kind == "genotype_lineage":
+        df = run_genotype_experiment(cfg)                     # E7: advantage of sex
+    elif kind == "society":
+        from .society import run_society                      # E8: multi-parent recombination
+        df = run_society(cfg)
+    else:
+        df = run_experiment(cfg)
     save_artifacts(cfg, df, out_dir)
     return out_dir
 

src/knowledge/genotype.py

@@ -0,0 +1,202 @@
+"""Multi-locus genotype space + recombination — the theoretical frame for the society.
+
+The single-locus, fixed-`p*` model can only express *recovery toward a ceiling*. The society's
+vertical claim — capability that *exceeds* any component — needs combinatorial structure. A
+**genotype** is `L` biallelic loci (allele 1 = "correct", 0 = "wrong"); a model's knowledge is a
+distribution over the `2^L` genotypes (a K-vector with `K = 2^L`, so all of Layer 1's K-mode
+machinery — drift, grounding, selection, metrics — applies unchanged). Fitness is additive (the
+number of correct loci); the optimum is the all-correct genotype.
+
+The one genuinely new operator is **recombination**. Free recombination replaces the joint genotype
+distribution by the product of its per-locus marginals (linkage equilibrium / Robbins proportions) —
+the population-level image of meiotic reassortment. This is what halts Muller's ratchet (it
+reconstitutes low-load genotypes from complementary high-load parents) and drives the Fisher–Muller
+effect (it assembles beneficial alleles that arose in different lineages into a genotype fitter than
+any parent). `rate=0` is asexual/clonal (Layer 1's regime); `rate=1` is free recombination.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .metrics import heterozygosity
+
+
+def genotype_bits(L: int) -> np.ndarray:
+    """Return the ``(2^L, L)`` matrix of genotype bit-vectors (row g = g in binary, MSB first).
+
+    Args:
+        L (int): Number of biallelic loci.
+
+    Returns:
+        np.ndarray: ``(2^L, L)`` int8 array; ``bits[g, l]`` is allele of locus ``l`` in genotype ``g``.
+    """
+    g = np.arange(1 << L, dtype=np.int64)
+    shifts = np.arange(L - 1, -1, -1, dtype=np.int64)
+    return ((g[:, None] >> shifts[None, :]) & 1).astype(np.int8)
+
+
+def additive_fitness(L: int) -> np.ndarray:
+    """Additive fitness vector: ``f(g)`` = number of correct (allele-1) loci in genotype ``g``.
+
+    Args:
+        L (int): Number of loci.
+
+    Returns:
+        np.ndarray: Length-``2^L`` fitness vector in ``{0, 1, ..., L}``; optimum = all-ones genotype.
+    """
+    return genotype_bits(L).sum(axis=1).astype(float)
+
+
+def locus_marginals(p: np.ndarray, L: int) -> np.ndarray:
+    """Per-locus frequency of the correct (allele-1) variant under distribution ``p``.
+
+    Args:
+        p (np.ndarray): Genotype distribution (length ``2^L``, sums to 1).
+        L (int): Number of loci.
+
+    Returns:
+        np.ndarray: Length-``L`` array; entry ``l`` is ``P(locus l is correct)``.
+    """
+    return genotype_bits(L).T.astype(float) @ np.asarray(p, dtype=float)   # (L, 2^L) @ (2^L,)
+
+
+def linkage_equilibrium(marginals: np.ndarray) -> np.ndarray:
+    """Build the product (linkage-equilibrium) genotype distribution from per-locus marginals.
+
+    ``p_LE(g) = ∏_l [q_l if g_l==1 else (1−q_l)]`` — the joint under free recombination, where the
+    loci are statistically independent (Robbins proportions).
+
+    Args:
+        marginals (np.ndarray): Length-``L`` correct-allele frequencies ``q_l``.
+
+    Returns:
+        np.ndarray: Length-``2^L`` product distribution (sums to 1).
+    """
+    q = np.asarray(marginals, dtype=float)
+    L = q.size
+    bits = genotype_bits(L)                                   # (2^L, L)
+    per = np.where(bits == 1, q[None, :], 1.0 - q[None, :])   # (2^L, L)
+    p = per.prod(axis=1)
+    s = p.sum()
+    return p / s if s > 0 else p
+
+
+def recombine(p: np.ndarray, L: int, rate: float) -> np.ndarray:
+    """Apply recombination at ``rate`` to a genotype distribution.
+
+    ``rate=0`` returns ``p`` unchanged (asexual / clonal); ``rate=1`` returns the full
+    linkage-equilibrium product of ``p``'s marginals (free recombination); intermediate rates
+    interpolate ``(1−rate)·p + rate·p_LE``.
+
+    Args:
+        p (np.ndarray): Genotype distribution (length ``2^L``).
+        L (int): Number of loci.
+        rate (float): Recombination rate in ``[0, 1]``.
+
+    Returns:
+        np.ndarray: The recombined distribution (sums to 1).
+    """
+    p = np.asarray(p, dtype=float)
+    if rate <= 0.0:
+        return p
+    p_le = linkage_equilibrium(locus_marginals(p, L))
+    out = p_le if rate >= 1.0 else (1.0 - rate) * p + rate * p_le
+    s = out.sum()
+    return out / s if s > 0 else out
+
+
+def recombine_teachers(teachers: list[np.ndarray], L: int, rate: float) -> np.ndarray:
+    """Merge ``K_T`` teacher genotype distributions (the multi-teacher / Fisher–Muller operator).
+
+    ``rate=0`` = clonal mean-mixture (``mean(teachers)`` — keeps genotypes intact, cannot create a
+    genotype no teacher had). ``rate=1`` = sexual merge: the linkage-equilibrium product of the
+    *pooled* per-locus marginals, which assembles the best allele of each locus across teachers into
+    genotypes none of them held — the Fisher–Muller effect and "merge, don't average" at the genotype
+    level.
+
+    Args:
+        teachers (list[np.ndarray]): ``K_T`` genotype distributions (each length ``2^L``).
+        L (int): Number of loci.
+        rate (float): Recombination rate in ``[0, 1]``.
+
+    Returns:
+        np.ndarray: The merged genotype distribution (sums to 1).
+    """
+    stack = np.asarray(teachers, dtype=float)
+    mean_mix = stack.mean(axis=0)
+    if rate <= 0.0:
+        return mean_mix
+    pooled_marginals = locus_marginals(mean_mix, L)          # pooled per-locus correct-allele freq
+    p_le = linkage_equilibrium(pooled_marginals)
+    out = p_le if rate >= 1.0 else (1.0 - rate) * mean_mix + rate * p_le
+    s = out.sum()
+    return out / s if s > 0 else out
+
+
+def mutate(p: np.ndarray, L: int, mu: float) -> np.ndarray:
+    """Apply per-locus symmetric mutation at rate ``mu`` to a genotype distribution.
+
+    Each locus independently flips with probability ``mu``. At the distribution level this convolves
+    ``p`` with the per-locus flip kernel; implemented as ``L`` independent locus mixings. ``mu=0`` is
+    a no-op. Provides the mutational pressure that Muller's ratchet grinds against.
+
+    Args:
+        p (np.ndarray): Genotype distribution (length ``2^L``).
+        L (int): Number of loci.
+        mu (float): Per-locus flip probability in ``[0, 0.5]``.
+
+    Returns:
+        np.ndarray: The mutated distribution (sums to 1).
+    """
+    p = np.asarray(p, dtype=float)
+    if mu <= 0.0:
+        return p
+    q = p.reshape([2] * L) if L > 0 else p
+    for axis in range(L):
+        flipped = np.flip(q, axis=axis)
+        q = (1.0 - mu) * q + mu * flipped
+    out = q.reshape(-1)
+    s = out.sum()
+    return out / s if s > 0 else out
+
+
+def locus_metrics(p: np.ndarray, L: int, fitness: np.ndarray, alive_eps: float = 1e-9) -> dict:
+    """Genotype-aware metrics for one generation (the multi-locus analogue of the K-mode row).
+
+    Args:
+        p (np.ndarray): Genotype distribution (length ``2^L``).
+        L (int): Number of loci.
+        fitness (np.ndarray): Length-``2^L`` fitness (typically :func:`additive_fitness`).
+        alive_eps (float): A genotype/allele is "present" if its probability exceeds this.
+
+    Returns:
+        dict: ``mean_fitness``, ``best_fitness`` (max fitness of any present genotype),
+        ``opt_freq`` (mass on the all-correct optimum), ``min_load`` (fewest wrong loci among present
+        genotypes = the Muller-ratchet state), ``mean_locus_correct`` (mean per-locus correct-allele
+        freq), ``locus_H`` (mean per-locus heterozygosity), and ``ld`` (mean pairwise linkage
+        disequilibrium |D|).
+    """
+    p = np.asarray(p, dtype=float)
+    present = p > alive_eps
+    q = locus_marginals(p, L)                                 # per-locus correct freq
+    bits = genotype_bits(L)
+    row = {
+        "mean_fitness": float(fitness @ p),
+        "best_fitness": float(fitness[present].max()) if present.any() else 0.0,
+        "opt_freq": float(p[-1]),                             # genotype 2^L-1 = all-ones optimum
+        "min_load": int(L - fitness[present].max()) if present.any() else L,
+        "mean_locus_correct": float(q.mean()),
+        "locus_H": float(np.mean([heterozygosity(np.array([1 - qi, qi])) for qi in q])),
+    }
+    # Mean pairwise linkage disequilibrium |D_ij| = |P(1_i,1_j) - q_i q_j| over locus pairs.
+    if L >= 2:
+        ds = []
+        for i in range(L):
+            for j in range(i + 1, L):
+                p_ij = float(p[(bits[:, i] == 1) & (bits[:, j] == 1)].sum())
+                ds.append(abs(p_ij - q[i] * q[j]))
+        row["ld"] = float(np.mean(ds))
+    else:
+        row["ld"] = 0.0
+    return row

src/knowledge/genotype_lineage.py

@@ -0,0 +1,82 @@
+"""Single-population genotype evolution — the advantage of sex (E7).
+
+A population (distribution over the ``2^L`` genotypes) adapts toward a multi-locus optimum under
+the composed generational step: **selection** (fitness-proportional, favouring correct alleles) +
+**drift** (finite resample of ``n``) + **mutation** (per-locus flips) + **recombination** (asexual
+``rate=0`` vs sexual ``rate>0``). Recombination reassorts beneficial alleles that arise in different
+sub-lineages into one genotype; without it (asexual) those alleles suffer *clonal interference* and
+adaptation is slower. So a sexual lineage climbs toward the optimum faster than an asexual one — the
+classical advantage of sex, and the dynamic counterpart of E8's one-shot multi-parent assembly.
+
+Reuses ``step.apply_selection`` (fitness = number of correct loci) and the ``genotype`` operators;
+emits the same tidy per-generation DataFrame contract as ``run_lineage`` (with genotype-aware
+columns from ``genotype.locus_metrics``).
+"""
+
+from __future__ import annotations
+
+from typing import Any, Mapping
+
+import numpy as np
+import pandas as pd
+
+from .genotype import additive_fitness, locus_metrics, mutate, recombine
+from .step import apply_selection
+
+
+def run_genotype_lineage(cfg: Mapping[str, Any], seed: int) -> pd.DataFrame:
+    """Run one genotype lineage and return per-generation genotype metrics.
+
+    Args:
+        cfg (Mapping): Config with a ``genotype`` block (``L``; drift ``n``; mutation ``mu``;
+            selection ``base`` for multiplicative fitness ``base^#correct``; recombination
+            ``recomb_rate``; ``init`` in {``wrong``, ``uniform``, ``optimum``}) and
+            ``generations``.
+        seed (int): Replicate seed; the run is a pure function of (cfg, seed).
+
+    Returns:
+        pd.DataFrame: One row per generation 0..T with ``generation`` plus the
+        ``genotype.locus_metrics`` columns (``mean_fitness``, ``best_fitness``, ``opt_freq``,
+        ``min_load``, ``mean_locus_correct``, ``locus_H``, ``ld``).
+    """
+    g = cfg["genotype"]
+    L = int(g["L"])
+    n = int(g["n"])
+    mu = float(g.get("mu", 0.0))
+    base = float(g.get("base", 1.0))
+    rate = float(g.get("recomb_rate", 0.0))
+    generations = int(cfg.get("generations", 100))
+    K = 1 << L
+
+    report_fitness = additive_fitness(L)             # # correct loci (0..L), for metrics
+    sel_fitness = base ** report_fitness             # multiplicative selection weight
+
+    init = g.get("init", "wrong")
+    p = np.zeros(K)
+    if init == "wrong":
+        p[0] = 1.0                                   # all-wrong genotype (load L); adapt upward
+    elif init == "optimum":
+        p[-1] = 1.0                                  # all-correct (for degradation studies)
+    elif init == "uniform":
+        p[:] = 1.0 / K
+    else:
+        raise ValueError(f"unknown genotype init {init!r} (expected wrong|optimum|uniform)")
+
+    rng = np.random.default_rng(seed)
+    rows: list[dict] = []
+
+    def record(t: int) -> None:
+        row = {"generation": t}
+        row.update(locus_metrics(p, L, report_fitness, alive_eps=1.0 / n))
+        rows.append(row)
+
+    record(0)
+    for t in range(1, generations + 1):
+        p = apply_selection(p, sel_fitness, "greedy", 0.0)   # fitness-proportional selection
+        counts = rng.multinomial(n, p)                       # drift
+        p = counts / counts.sum()
+        p = mutate(p, L, mu)                                 # per-locus mutation
+        p = recombine(p, L, rate)                            # asexual (0) vs sexual (>0)
+        record(t)
+
+    return pd.DataFrame(rows)

src/knowledge/society.py

@@ -0,0 +1,95 @@
+"""Multi-parent recombination — the Fisher–Muller vertical claim (E8).
+
+The society's headline: **an offspring recombined from many decorrelated parents can be fitter
+than any parent** — capability that *exceeds* every component, not just recovers a ceiling. This
+is the Fisher–Muller effect, and unlike biological sex it has **no two-parent limit**: an offspring
+here can have arbitrarily many parents (``recombine_teachers`` pools per-locus marginals over all of
+them).
+
+Each parent is a *specialist*: confident-correct on the loci it has mastered, agnostic (~0.5) on the
+rest — the realistic picture of an expert. Which loci each parent masters comes from the exact
+shared-switch construction (``teachers.make_retention_matrix``), so parent count ``K_T`` and
+decorrelation ``rho`` are clean, independently-swept knobs (mastery of a *locus* replaces retention
+of a *tail item*). Three deployable capabilities are compared, all as the fitness of the *mode*
+(most-probable) genotype — what you would actually ship:
+
+* **best_parent** — the single fittest specialist (no combination).
+* **average** — the mean-mixture "model soup" (combine, but don't recombine loci).
+* **sexual** — the union-preserving recombination (assemble the best allele of each locus).
+
+Sexual climbs to the optimum as parents accumulate and decorrelate; the other two plateau.
+"""
+
+from __future__ import annotations
+
+import itertools
+
+import numpy as np
+import pandas as pd
+
+from .genotype import additive_fitness, linkage_equilibrium, recombine_teachers
+from .seeding import spawn_seeds
+from .teachers import make_retention_matrix
+
+
+def make_specialist(mastery: np.ndarray, hi: float, lo: float) -> np.ndarray:
+    """Build a specialist's genotype distribution from its per-locus mastery mask.
+
+    Args:
+        mastery (np.ndarray): Length-``L`` bool; True where this parent is expert.
+        hi (float): Correct-allele probability on mastered loci (confident).
+        lo (float): Correct-allele probability on unmastered loci (agnostic, ~0.5).
+
+    Returns:
+        np.ndarray: Length-``2^L`` genotype distribution (product / linkage-equilibrium form).
+    """
+    q = np.where(np.asarray(mastery, dtype=bool), hi, lo)
+    return linkage_equilibrium(q)
+
+
+def _mode_fitness(p: np.ndarray, fitness: np.ndarray) -> float:
+    """Fitness of the most-probable (deployed consensus) genotype."""
+    return float(fitness[int(np.argmax(p))])
+
+
+def run_society(cfg: dict) -> pd.DataFrame:
+    """Sweep parent count x decorrelation; compare best-parent vs average vs sexual capability.
+
+    Args:
+        cfg (dict): Parsed experiment config with a ``society`` block (``L``, ``q`` mastery
+            fraction, ``hi``, ``lo``), a ``sweep`` (``K_T`` x ``rho``), ``seed``, ``n_replicates``.
+
+    Returns:
+        pd.DataFrame: One row per (K_T, rho, replicate) with ``best_parent``, ``average``,
+        ``sexual`` (mode-genotype fitness, 0..L) and ``L``.
+    """
+    soc = cfg["society"]
+    L = int(soc["L"])
+    q = float(soc.get("q", 0.5))
+    hi, lo = float(soc.get("hi", 0.9)), float(soc.get("lo", 0.45))
+    fitness = additive_fitness(L)
+
+    sweeps = cfg["sweep"]
+    if isinstance(sweeps, dict):
+        sweeps = [sweeps]
+    params = [s["param"] for s in sweeps]
+    value_lists = [list(s["values"]) for s in sweeps]
+    seeds = spawn_seeds(int(cfg["seed"]), int(cfg["n_replicates"]))
+
+    rows: list[dict] = []
+    for combo in itertools.product(*value_lists):
+        d = dict(zip(params, combo))
+        K_T, rho = int(d["K_T"]), float(d["rho"])
+        for rep, ss in enumerate(seeds):
+            child = int(ss.generate_state(1)[0])
+            rng = np.random.default_rng(child)
+            mastery = make_retention_matrix(L, K_T, rho, q, rng).astype(bool)  # (K_T, L)
+            parents = [make_specialist(mastery[k], hi, lo) for k in range(K_T)]
+            rows.append({
+                "experiment": cfg["experiment"], "K_T": K_T, "rho": rho, "q": q, "L": L,
+                "replicate": rep,
+                "best_parent": max(_mode_fitness(p, fitness) for p in parents),
+                "average": _mode_fitness(recombine_teachers(parents, L, 0.0), fitness),
+                "sexual": _mode_fitness(recombine_teachers(parents, L, 1.0), fitness),
+            })
+    return pd.DataFrame(rows)