society: make the sexual-transmission model rigorous (E9 epistasis, E10 directed sex)

Deepen the sexual-reproduction frame before entering the full society, on
the two facets GG chose: landscape robustness and directed recombination.
Adds a Kauffman NK landscape (genotype.nk_fitness, tunable ruggedness),
finite n-parent crossover (genotype.crossover, per-gap recombination rate),
and hill-climb (parents = local optima = trained models).

E9 (recomb_landscape) -- the "why sex?" test: E8's dramatic super-parent
result used an ADDITIVE landscape. On rugged/epistatic landscapes, blindly
recombining local optima causes OUTBREEDING DEPRESSION -- offspring fall
below the parents, worse with both ruggedness and recombination rate (K=8,
free recomb: ~ -0.23), and the optimal recombination rate shrinks as
ruggedness grows. Design rule: merge freely when skills are complementary/
additive; sparingly (and with selection) when entangled.

E10 (directed_sex) -- directed sex beats biological sex: biology is stuck
with 2 random-mating parents and no offspring preview; an AI can choose
complementary mates, evaluate many recombinant offspring, keep the fittest,
and use unbounded parents (iterated recombine-then-select). Random
("biological") sex craters with ruggedness (0.66->0.51); directed sex
tracks/exceeds the best parent at every ruggedness -- converting the
outbreeding-depression catastrophe into a win. No biological analog.

Complete sexual-transmission picture: dramatic super-parent offspring when
skills are complementary (E8); outbreeding-depression risk when entangled
(E9); directed sex resolves the risk (E10). configs/layer1/{E9,E10}.yaml,
figures/plot_{E9,E10}.py, READMEs, +5 tests (117 green).

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

src/knowledge/experiment.py

@@ -326,6 +326,12 @@ def run_and_save(config_path: str | Path) -> Path:
     elif kind == "society":
         from .society import run_society                      # E8: multi-parent recombination
         df = run_society(cfg)
+    elif kind == "recomb_landscape":
+        from .society import run_recomb_landscape             # E9: landscape robustness / epistasis
+        df = run_recomb_landscape(cfg)
+    elif kind == "directed_sex":
+        from .society import run_directed_sex                 # E10: directed sex beats biology
+        df = run_directed_sex(cfg)
     else:
         df = run_experiment(cfg)
     save_artifacts(cfg, df, out_dir)

src/knowledge/genotype.py

@@ -48,6 +48,90 @@ def additive_fitness(L: int) -> np.ndarray:
     return genotype_bits(L).sum(axis=1).astype(float)
 
 
+def nk_fitness(L: int, K: int, seed: int) -> np.ndarray:
+    """Kauffman NK fitness landscape over the ``2^L`` genotypes (tunable ruggedness).
+
+    Each locus contributes a random value in [0,1] that depends on its own allele plus the ``K``
+    following loci (adjacent ring neighbourhood); total fitness = mean of the ``L`` contributions.
+    ``K=0`` is additive/smooth (a single peak, recombination unambiguously helps); larger ``K`` is
+    epistatic/rugged (many local optima, co-adapted allele blocks that recombination can break —
+    the regime where sex can hurt).
+
+    Args:
+        L (int): Number of loci.
+        K (int): Epistatic interactions per locus (``0..L-1``); ruggedness knob.
+        seed (int): Landscape seed (the landscape is a deterministic function of it).
+
+    Returns:
+        np.ndarray: Length-``2^L`` fitness vector in [0,1].
+    """
+    rng = np.random.default_rng(seed)
+    bits = genotype_bits(L).astype(np.int64)
+    F = np.zeros(1 << L)
+    for i in range(L):
+        idx = [(i + j) % L for j in range(K + 1)]              # locus i + its K ring-neighbours
+        table = rng.random(1 << (K + 1))                       # random contribution per pattern
+        pat = np.zeros(1 << L, dtype=np.int64)
+        for b, locus in enumerate(idx):
+            pat |= bits[:, locus] << (K - b)
+        F += table[pat]
+    return F / L
+
+
+def hill_climb(fitness: np.ndarray, L: int, start: int) -> int:
+    """Greedy single-locus-flip hill-climb to a local optimum (a "trained specialist" parent).
+
+    Args:
+        fitness (np.ndarray): Length-``2^L`` fitness vector.
+        L (int): Number of loci.
+        start (int): Starting genotype index.
+
+    Returns:
+        int: A local-optimum genotype index (no single-locus flip improves fitness).
+    """
+    g = int(start)
+    while True:
+        neighbours = [g ^ (1 << l) for l in range(L)]
+        best = max(neighbours, key=lambda x: fitness[x])
+        if fitness[best] <= fitness[g]:
+            return g
+        g = best
+
+
+def crossover(parents_bits: np.ndarray, rate: float, rng: np.random.Generator) -> np.ndarray:
+    """One recombinant offspring from ``n`` parents by per-gap switching (finite, stochastic).
+
+    Walks the loci left to right inheriting from a current parent; at each gap the current parent is
+    re-drawn uniformly among the ``n`` parents with probability ``rate``. ``rate=0`` clones one
+    parent (asexual); ``rate=0.5`` gives near-independent per-locus inheritance (free recombination);
+    small ``rate`` preserves linked blocks of co-adapted alleles (the knob that matters on rugged
+    landscapes). Generalises 2-parent crossover to arbitrarily many parents (no two-parent limit).
+
+    Args:
+        parents_bits (np.ndarray): ``(n, L)`` bit-matrix of the parent genotypes.
+        rate (float): Per-gap recombination probability in ``[0, 0.5]``.
+        rng (np.random.Generator): Explicit RNG.
+
+    Returns:
+        np.ndarray: Length-``L`` offspring bit-vector.
+    """
+    n, L = parents_bits.shape
+    cur = int(rng.integers(n))
+    child = np.empty(L, dtype=parents_bits.dtype)
+    for l in range(L):
+        if l > 0 and rng.random() < rate:
+            cur = int(rng.integers(n))
+        child[l] = parents_bits[cur, l]
+    return child
+
+
+def bits_to_index(bits: np.ndarray) -> int:
+    """Convert a genotype bit-vector (MSB first) to its integer index."""
+    L = bits.size
+    weights = 1 << np.arange(L - 1, -1, -1)
+    return int(np.asarray(bits) @ weights)
+
+
 def locus_marginals(p: np.ndarray, L: int) -> np.ndarray:
     """Per-locus frequency of the correct (allele-1) variant under distribution ``p``.
 

src/knowledge/society.py

@@ -27,7 +27,10 @@ import itertools
 import numpy as np
 import pandas as pd
 
-from .genotype import additive_fitness, linkage_equilibrium, recombine_teachers
+from .genotype import (
+    additive_fitness, bits_to_index, crossover, genotype_bits, hill_climb, linkage_equilibrium,
+    nk_fitness, recombine_teachers,
+)
 from .seeding import spawn_seeds
 from .teachers import make_retention_matrix
 
@@ -52,6 +55,120 @@ def _mode_fitness(p: np.ndarray, fitness: np.ndarray) -> float:
     return float(fitness[int(np.argmax(p))])
 
 
+def _local_optima(fitness: np.ndarray, L: int, n_parents: int,
+                  rng: np.random.Generator) -> np.ndarray:
+    """``n_parents`` local-optimum genotypes (trained specialists) as an ``(n_parents, L)`` matrix."""
+    bits = genotype_bits(L)
+    opts = [hill_climb(fitness, L, int(rng.integers(1 << L))) for _ in range(n_parents)]
+    return np.stack([bits[g] for g in opts])
+
+
+def _random_sex(fitness, pbits, rate, n_off, rng) -> float:
+    """Mean fitness of ``n_off`` blind recombinant offspring (biology: random mating, no selection)."""
+    return float(np.mean([fitness[bits_to_index(crossover(pbits, rate, rng))]
+                          for _ in range(n_off)]))
+
+
+def _directed_sex(fitness, pbits, rate, pop, keep, rounds, rng) -> float:
+    """Best fitness reachable by *directed* sex: iterated recombine-then-select-offspring.
+
+    The AI superpower biology lacks — evaluate many recombinants and keep only the fittest, over
+    several rounds (mate choice + offspring selection + unbounded parents). Returns the best fitness
+    found.
+    """
+    pool = pbits.copy()
+    best = max(float(fitness[bits_to_index(b)]) for b in pool)
+    for _ in range(rounds):
+        offs = np.stack([crossover(pool, rate, rng) for _ in range(pop)])
+        fits = np.array([fitness[bits_to_index(o)] for o in offs])
+        pool = offs[np.argsort(fits)[-keep:]]
+        best = max(best, float(fits.max()))
+    return best
+
+
+def _sweep_grid(cfg):
+    """Return (params, value_lists, seeds) for the config's sweep."""
+    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"]))
+    return params, value_lists, seeds
+
+
+def run_recomb_landscape(cfg: dict) -> pd.DataFrame:
+    """E9 — landscape robustness: when does recombination help, and at what rate?
+
+    Parents are local optima ("trained models") of a Kauffman NK landscape whose ruggedness ``K``
+    (epistasis) is swept together with the recombination rate. On smooth/mildly-rugged landscapes
+    recombination helps; on rugged ones free recombination breaks co-adapted blocks and offspring
+    fall *below* the parents (outbreeding depression), with the optimal rate shrinking as ruggedness
+    grows. Sweep ``K`` x ``rate``.
+
+    Returns:
+        pd.DataFrame: rows with ``K``, ``rate``, ``best_parent``, ``mean_offspring``,
+        ``best_offspring``, ``global_opt`` (all NK fitness in [0,1]).
+    """
+    soc = cfg["society"]
+    L, nP, pop = int(soc["L"]), int(soc["n_parents"]), int(soc.get("pop", 200))
+    params, value_lists, seeds = _sweep_grid(cfg)
+    rows: list[dict] = []
+    for combo in itertools.product(*value_lists):
+        d = dict(zip(params, combo))
+        K, rate = int(d["K"]), float(d["rate"])
+        for rep, ss in enumerate(seeds):
+            seed = int(ss.generate_state(1)[0])
+            fitness = nk_fitness(L, K, seed)
+            rng = np.random.default_rng(seed)
+            pbits = _local_optima(fitness, L, nP, rng)
+            offs = np.array([fitness[bits_to_index(crossover(pbits, rate, rng))]
+                             for _ in range(pop)])
+            best_parent = max(float(fitness[bits_to_index(b)]) for b in pbits)
+            rows.append({
+                "experiment": cfg["experiment"], "K": K, "rate": rate, "replicate": rep,
+                "best_parent": best_parent, "mean_offspring": float(offs.mean()),
+                "best_offspring": float(offs.max()), "global_opt": float(fitness.max())})
+    return pd.DataFrame(rows)
+
+
+def run_directed_sex(cfg: dict) -> pd.DataFrame:
+    """E10 — directed sex beats biological sex: mate choice + offspring selection rescue ruggedness.
+
+    Across landscape ruggedness ``K``, compare the deployed capability of: the best single parent;
+    **random sex** (blind mating, no selection — biology's default, which suffers outbreeding
+    depression on rugged landscapes); and **directed sex** (iterated recombine-then-select, the
+    AI-only move). Directed sex avoids the random-sex catastrophe and matches or exceeds the best
+    parent even when skills are entangled.
+
+    Returns:
+        pd.DataFrame: rows with ``K``, ``best_parent``, ``random_sex``, ``directed_sex``,
+        ``global_opt``.
+    """
+    soc = cfg["society"]
+    L, nP = int(soc["L"]), int(soc["n_parents"])
+    rate = float(soc.get("rate", 0.2))
+    n_off = int(soc.get("pop", 200))
+    keep, rounds = int(soc.get("keep", 8)), int(soc.get("rounds", 5))
+    params, value_lists, seeds = _sweep_grid(cfg)
+    rows: list[dict] = []
+    for combo in itertools.product(*value_lists):
+        d = dict(zip(params, combo))
+        K = int(d["K"])
+        for rep, ss in enumerate(seeds):
+            seed = int(ss.generate_state(1)[0])
+            fitness = nk_fitness(L, K, seed)
+            rng = np.random.default_rng(seed)
+            pbits = _local_optima(fitness, L, nP, rng)
+            rows.append({
+                "experiment": cfg["experiment"], "K": K, "replicate": rep,
+                "best_parent": max(float(fitness[bits_to_index(b)]) for b in pbits),
+                "random_sex": _random_sex(fitness, pbits, 0.5, n_off, rng),
+                "directed_sex": _directed_sex(fitness, pbits, rate, n_off, keep, rounds, rng),
+                "global_opt": float(fitness.max())})
+    return pd.DataFrame(rows)
+
+
 def run_society(cfg: dict) -> pd.DataFrame:
     """Sweep parent count x decorrelation; compare best-parent vs average vs sexual capability.