worldhistory: initial public history

19 files changed, 2478 insertions, 0 deletions

diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/tests/__init__.py
diff --git a/tests/helpers.py b/tests/helpers.py
new file mode 100644
index 0000000..28b6266
--- /dev/null
+++ b/tests/helpers.py
@@ -0,0 +1,68 @@
+import functools
+
+import numpy as np
+
+from worldhistory.world import World, neighbours_from_ids
+
+DEFAULTS = dict(ocean=False, elevation_m=100.0, T_mean=15.0, P_ann=1000.0, holdridge=20, landform=1, river=False,
+                strahler=0, lake=False, ice=0, coal_potential=False, iron_potential=False, vent_potential=0.0,
+                seabed_type=0, lithology=1, ground=0, gravity_g=1.0, o2_fraction=0.21, pressure_bar=1.0,
+                bottom_temp_c=2.0, plate=0)
+
+
+def fields(n, **over):
+    out = {}
+    for k, v in {**DEFAULTS, **over}.items():
+        a = np.asarray(v)
+        out[k] = np.broadcast_to(a, (n,)).copy() if a.ndim == 0 else a.copy()
+    return out
+
+
+def line_world(n, spacing_km=100.0, radius_km=6371.0, **over):
+    """Cells on the equator in a row, each linked to the next: 1-D worlds with hand-checkable answers."""
+    lon = np.arange(n) * np.degrees(spacing_km / radius_km)
+    lat = np.zeros(n)
+    xyz = np.stack([np.cos(np.radians(lon)), np.sin(np.radians(lon)), np.zeros(n)], 1)
+    nb = np.full((n, 6), -1, np.int64)
+    nb[1:, 0] = np.arange(n - 1)
+    nb[:-1, 1] = np.arange(1, n)
+    return World(nb=nb, lat=lat, lon=lon, xyz=xyz, area=np.full(n, spacing_km ** 2), radius_km=radius_km,
+                 base=fields(n, **over))
+
+
+@functools.lru_cache(maxsize=4)
+def _h3_grid(res):
+    import h3.api.basic_int as h3
+    ids = sorted(c for r0 in h3.get_res0_cells() for c in h3.cell_to_children(r0, res))
+    ll = np.array([h3.cell_to_latlng(c) for c in ids])
+    area = np.array([h3.cell_area(c, unit="rads^2") for c in ids])
+    return np.array(ids, np.uint64), ll, area
+
+
+def globe_world(res=1, radius_km=6371.0):
+    """A whole H3 globe (res 1 = 842 cells) with default fields; tests edit world.base afterwards."""
+    ids, ll, area = _h3_grid(res)
+    lat, lon = ll[:, 0].copy(), ll[:, 1].copy()
+    la, lo = np.radians(lat), np.radians(lon)
+    xyz = np.stack([np.cos(la) * np.cos(lo), np.cos(la) * np.sin(lo), np.sin(la)], 1)
+    return World(nb=neighbours_from_ids(ids), lat=lat, lon=lon, xyz=xyz, area=area * radius_km ** 2,
+                 radius_km=radius_km, base=fields(len(ids)))
+
+
+from worldhistory.config import parse_history, parse_race
+
+
+def make_race(id="a", **over):
+    """A land race with density 1/km² at q=1, growth 0.1/step, no Allee, founder 1, no mobility; sections override."""
+    d = {"id": id,
+         "demography": {"density": 1.0, "growth": 0.1, "allee": 0.0, "founder": 1.0, "mobility": 0.0},
+         "habitat": {"realm": "land", "terms": [{"p": "land", "w": 1.0}]}}
+    for k, v in over.items():
+        d[k] = {**d.get(k, {}), **v} if isinstance(v, dict) and isinstance(d.get(k), dict) else v
+    return parse_race(d, f"test:{id}")
+
+
+def make_history(**over):
+    d = {"run": {"years": 100, "step": 10, "snapshot_every": 50, "seed": 1}}
+    d.update(over)
+    return parse_history(d, "test:history")
diff --git a/tests/test_adaptation.py b/tests/test_adaptation.py
new file mode 100644
index 0000000..117a39e
--- /dev/null
+++ b/tests/test_adaptation.py
@@ -0,0 +1,119 @@
+import unittest
+
+import numpy as np
+
+from tests.helpers import line_world, make_race
+from worldhistory.adaptation import adapt
+from worldhistory.config import CONDITIONS
+from worldhistory.habitat import environment, fitness, native_optima
+from worldhistory.state import new_state
+
+D = CONDITIONS.index("depth")
+TOL = {"depth": {"optimum": 0, "width": 50, "lo": 0, "hi": 600, "rate": 2.0}}
+
+
+def world_and_state(race, pops):
+    w = line_world(4, ocean=[False, True, True, True], elevation_m=[10, -200, -600, -1000])
+    st = new_state(1, 4, native_optima([race]), seed=0)
+    st.P[0] = pops
+    return w, st
+
+
+class AdaptationTest(unittest.TestCase):
+    def test_moves_at_most_rate_times_years(self):
+        r = make_race(tolerance=TOL)
+        w, st = world_and_state(r, [0, 10, 0, 0])
+        adapt(st, environment(w), [r], np.ones((1, 4)), 10)
+        self.assertAlmostEqual(float(st.O[0, D, 1]), 20.0)         # 2 m/yr × 10 yr toward the deep end (400)
+        self.assertEqual(float(st.O[0, D, 0]), 0.0)                # empty cells unchanged
+
+    def test_multiplier_speeds_adaptation(self):
+        r = make_race(tolerance=TOL)
+        w, st = world_and_state(r, [0, 10, 0, 0])
+        adapt(st, environment(w), [r], np.ones((1, 4)), 10, mult=np.full((1, 4), 2.0))
+        self.assertAlmostEqual(float(st.O[0, D, 1]), 40.0)         # pressed: 2 × (2 m/yr × 10 yr)
+
+    def test_turnover_slows_poor_habitat(self):
+        r = make_race(tolerance=TOL)
+        w, st = world_and_state(r, [0, 10, 0, 0])
+        adapt(st, environment(w), [r], np.zeros((1, 4)), 10)
+        self.assertAlmostEqual(float(st.O[0, D, 1]), 6.0)          # × 0.3
+
+    def test_soft_limit_brakes_but_does_not_stop(self):
+        r = make_race(tolerance=TOL)                               # hi 600, width 50
+        w, st = world_and_state(r, [0, 0, 0, 10])
+        st.O[0, D, 3] = 650.0                                      # 1 width past hi
+        adapt(st, environment(w), [r], np.ones((1, 4)), 10)
+        self.assertAlmostEqual(float(st.O[0, D, 3]), 650.0 + 20.0 * np.exp(-0.5), places=4)
+        for _ in range(200):
+            adapt(st, environment(w), [r], np.ones((1, 4)), 10)
+        self.assertGreater(float(st.O[0, D, 3]), 650.0)
+        self.assertLess(float(st.O[0, D, 3]), 800.0)               # far slower than the free rate (2 m/yr)
+
+    def test_inward_moves_are_not_braked(self):
+        r = make_race(tolerance=TOL)
+        w = line_world(3)                                          # all land: target depth 0
+        st = new_state(1, 3, native_optima([r]), seed=0)
+        st.P[0] = 10
+        st.O[0, D] = 700.0                                         # past hi, moving back in
+        adapt(st, environment(w), [r], np.ones((1, 3)), 10)
+        np.testing.assert_allclose(st.O[0, D], 680.0)
+
+    def test_rain_does_not_drift_at_sea(self):
+        r = make_race(habitat={"realm": "both", "terms": [{"p": "sea", "w": 1}]},
+                      tolerance={"rain": {"optimum": 3.0, "width": 0.3, "lo": 1.5, "hi": 3.5, "rate": 0.01}})
+        w = line_world(3, ocean=True, elevation_m=-100)
+        st = new_state(1, 3, native_optima([r]), seed=0)
+        st.P[0] = 10
+        adapt(st, environment(w), [r], np.ones((1, 3)), 10)
+        np.testing.assert_allclose(st.O[0, CONDITIONS.index("rain")], 3.0)
+
+    def test_inland_returns_to_native(self):
+        r = make_race(tolerance=TOL)
+        w = line_world(3)                                          # all land: depth span [0, 0]
+        st = new_state(1, 3, native_optima([r]), seed=0)
+        st.P[0] = 10
+        st.O[0, D] = 100.0
+        adapt(st, environment(w), [r], np.ones((1, 3)), 10)
+        np.testing.assert_allclose(st.O[0, D], 80.0)
+
+    def test_adapting_costs_the_old(self):
+        r = make_race(tolerance=TOL)
+        w, st = world_and_state(r, [0, 0, 0, 10])
+        env = environment(w)
+        O_before = st.O[0].copy()
+        for _ in range(60):
+            adapt(st, env, [r], np.ones((1, 4)), 10)
+        at_surface_before = fitness(env, O_before, r)[0]
+        at_surface_after = fitness(env, np.repeat(st.O[0][:, 3:4], 4, 1), r)[0]
+        self.assertAlmostEqual(float(at_surface_before), 1.0)
+        self.assertLess(float(at_surface_after), 0.1)              # deep-adapted people cannot live at the surface
+
+
+class ProspectiveTest(unittest.TestCase):
+    def test_empty_cells_take_neighbour_optimum(self):
+        from worldhistory.adaptation import prospective
+        r = make_race(tolerance=TOL)
+        w, st = world_and_state(r, [0, 30, 0, 10])
+        st.O[0, D] = [0.0, 100.0, 0.0, 400.0]
+        prospective(st, w)
+        # cell 0: only neighbour 1 occupied -> 100; cell 2: (30*100 + 10*400) / 40 = 175; occupied cells unchanged
+        np.testing.assert_allclose(st.O[0, D], [100.0, 100.0, 175.0, 400.0])
+
+    def test_isolated_empty_cells_keep_their_value(self):
+        from worldhistory.adaptation import prospective
+        r = make_race(tolerance=TOL)
+        w, st = world_and_state(r, [0, 0, 0, 0])
+        st.O[0, D] = [5.0, 6.0, 7.0, 8.0]
+        prospective(st, w)
+        np.testing.assert_allclose(st.O[0, D], [5.0, 6.0, 7.0, 8.0])
+
+    def test_fractional_groups_keep_their_optimum(self):
+        # all groups count (spec §1.3): a thin frontier (P < 1) is not "empty"
+        from worldhistory.adaptation import prospective
+        r = make_race(tolerance=TOL)
+        w, st = world_and_state(r, [100, 0.5, 0, 0])
+        st.O[0, D] = [0.0, 900.0, 0.0, 0.0]
+        prospective(st, w)
+        self.assertEqual(float(st.O[0, D, 1]), 900.0)
+        self.assertEqual(float(st.O[0, D, 2]), 900.0)          # the empty cell beyond takes the frontier's optimum
diff --git a/tests/test_cli.py b/tests/test_cli.py
new file mode 100644
index 0000000..e5a464c
--- /dev/null
+++ b/tests/test_cli.py
@@ -0,0 +1,44 @@
+import json
+import tempfile
+import unittest
+from pathlib import Path
+
+import numpy as np
+
+from tests.helpers import globe_world
+from tests.test_engine import small_setup
+from worldhistory.cli import main
+from worldhistory.engine import Engine
+from worldhistory.preview import compare, preview
+
+
+class PreviewTest(unittest.TestCase):
+    def test_preview_and_compare(self):
+        with tempfile.TemporaryDirectory() as d:
+            a, b = Path(d) / "a", Path(d) / "b"
+            w, h, races = small_setup(1)
+            Engine(w, h, races).run(a)
+            w2, h2, races2 = small_setup(2)
+            Engine(w2, h2, races2).run(b)
+            files = preview(a, world=w)
+            names = {p.name for p in files}
+            self.assertIn("race_y0060.png", names)
+            self.assertIn("chart.png", names)
+            self.assertIn("contact.png", names)
+            same = compare(a, a, world=w)
+            self.assertEqual(same["agreement"], 1.0)
+            diff = compare(a, b, world=w, out=Path(d) / "cmp")
+            self.assertLessEqual(diff["agreement"], 1.0)
+            self.assertTrue((Path(d) / "cmp" / "compare.png").exists())
+
+
+class CliTest(unittest.TestCase):
+    def test_check_config_errors(self):
+        with tempfile.TemporaryDirectory() as d:
+            (Path(d) / "races").mkdir()
+            (Path(d) / "history.toml").write_text("[run]\nyers = 5\n")
+            (Path(d) / "races" / "a.toml").write_text('id = "a"\n[demography]\ndensity = 1\ngrowth = 0.1\n'
+                                                      '[habitat]\nterms = [{ p = "land", w = 1 }]\n')
+            self.assertEqual(main(["check-config", d]), 1)
+            (Path(d) / "history.toml").write_text("[run]\nyears = 5\n")
+            self.assertEqual(main(["check-config", d]), 0)
diff --git a/tests/test_config.py b/tests/test_config.py
new file mode 100644
index 0000000..2a76d8d
--- /dev/null
+++ b/tests/test_config.py
@@ -0,0 +1,194 @@
+import tempfile
+import tomllib
+import unittest
+from pathlib import Path
+
+from tests.helpers import make_history, make_race
+from worldhistory.config import ConfigError, dumps_toml, link, load_config, parse_race, resolved
+
+RACE = """
+id = "stone"
+name = "Stone folk"
+colour = [140, 140, 150]
+[demography]
+density = 1.2
+growth = 0.02
+allee = 75
+[seed]
+clusters = 25
+heads = 200
+[habitat]
+realm = "land"
+terms = [{ p = "mountain", w = 0.7 }, { p = ["land", "coal"], w = 0.5 }]
+gates = [{ penalty = ["above:o2_fraction:0.3", "lithology:3,4"], factor = 0.6 }]
+[tolerance.temperature]
+optimum = 10
+width = 20
+lo = -30
+hi = 40
+rate = 0.001
+[special]
+veins = 0.002
+"""
+
+SUB = """
+id = "deep"
+[demography]
+density = 0.5
+growth = 0.08
+[habitat]
+realm = "sea"
+terms = [{ p = "vent", w = 1.0 }]
+qmax = 0.7
+[emerge]
+parent = "stone"
+birth_min = 0.5
+[tolerance.temperature]
+optimum = 30
+width = 10
+lo = 0
+hi = 50
+"""
+
+HISTORY = """
+[run]
+years = 7500
+seed = 3
+[world]
+eras = ["late"]
+[regions.zone]
+kind = "changed"
+era = "late"
+fields = ["gravity_g"]
+[[event]]
+year = 0
+kind = "seed"
+[[event]]
+year = 20
+kind = "cull"
+region = "zone"
+share = 0.6
+by_race = { stone = 0.8 }
+[[event]]
+year = 20
+kind = "era_switch"
+era = "late"
+"""
+
+
+def write(d, name, text):
+    p = Path(d) / name
+    p.parent.mkdir(parents=True, exist_ok=True)
+    p.write_text(text)
+
+
+class ConfigTest(unittest.TestCase):
+    def test_load_defaults_and_inheritance(self):
+        with tempfile.TemporaryDirectory() as d:
+            write(d, "history.toml", HISTORY)
+            write(d, "races/stone.toml", RACE)
+            write(d, "races/deep.toml", SUB)
+            h, races = load_config(d)
+        self.assertEqual([r.id for r in races], ["stone", "deep"])
+        stone, deep = races
+        self.assertEqual(h.step, 10)                       # default
+        self.assertEqual(h.years, 7500)
+        self.assertEqual(stone.founder, 75)                # founder defaults to allee
+        self.assertEqual(stone.family, "stone")
+        self.assertEqual(deep.family, "stone")             # sub-species joins the parent's family
+        self.assertEqual(deep.tolerance["temperature"].width, 10)   # own curve
+        self.assertEqual(deep.clusters, 0)                 # sub-species are never seeded
+        self.assertEqual(stone.travel["river"], -1.0)      # default
+        self.assertEqual(stone.gates[0]["factor"], 0.6)
+        self.assertEqual(h.events[1]["by_race"], {"stone": 0.8})
+
+    def test_unknown_key(self):
+        with self.assertRaisesRegex(ConfigError, r"test:x.*demography.*desnity"):
+            parse_race({"id": "x", "demography": {"desnity": 1}}, "test:x")
+
+    def test_unknown_predicate(self):
+        with self.assertRaisesRegex(ConfigError, "nosuch"):
+            make_race(habitat={"terms": [{"p": "nosuch", "w": 1}]})
+
+    def test_bad_tolerance_condition(self):
+        with self.assertRaisesRegex(ConfigError, "wind"):
+            make_race(tolerance={"wind": {"optimum": 0, "width": 1, "lo": 0, "hi": 1}})
+
+    def test_tolerance_needs_widths(self):
+        with self.assertRaisesRegex(ConfigError, "width"):
+            make_race(tolerance={"depth": {"optimum": 0, "lo": 0, "hi": 10, "width_lo": 5}})
+
+    def test_link_checks_references(self):
+        h = make_history(event=[{"year": 20, "kind": "cull", "region": "nowhere", "share": 0.5}])
+        with self.assertRaisesRegex(ConfigError, "nowhere"):
+            link(h, [make_race()])
+        h = make_history(event=[{"year": 20, "kind": "cull", "share": 0.5, "by_race": {"ghost": 0.1}}])
+        with self.assertRaisesRegex(ConfigError, "ghost"):
+            link(h, [make_race()])
+        h = make_history(event=[{"year": 20, "kind": "era_switch", "era": "late"}])
+        with self.assertRaisesRegex(ConfigError, "late"):
+            link(h, [make_race()])
+        sub = make_race("s", emerge={"parent": "p"})
+        with self.assertRaisesRegex(ConfigError, "p"):
+            link(make_history(), [sub])
+        same = {"gravity": {"optimum": 1.0, "width": 0.5, "lo": 0.3, "hi": 1.3}}
+        with self.assertRaisesRegex(ConfigError, "never be born"):            # progress stays 0 < birth_min
+            link(make_history(), [make_race("p", tolerance=same), make_race("s", tolerance=same, emerge={"parent": "p"})])
+        link(make_history(), [make_race("p", tolerance=same),                 # birth_min 0 or ritual: fine
+                              make_race("s", emerge={"parent": "p", "birth_min": 0.0, "birth_sure": 0.0})])
+        link(make_history(), [make_race("p", tolerance=same), make_race("s", emerge={"parent": "p", "mode": "ritual"})])
+
+    def test_emerge_defaults_and_merge(self):
+        with tempfile.TemporaryDirectory() as d:
+            write(d, "history.toml", HISTORY)
+            write(d, "races/stone.toml", RACE.replace("[special]", "[tolerance.depth]\noptimum = 0\nwidth = 60\n"
+                                                                    "lo = 0\nhi = 6000\n[special]"))
+            write(d, "races/deep.toml", SUB)
+            h, (stone, deep) = load_config(d)
+        em = deep.emerge
+        self.assertEqual((em["mode"], em["birth_min"], em["birth_sure"], em["birth_rate"]), ("adapt", 0.5, 0.8, 0.2))
+        self.assertEqual((em["return_min"], em["return_sure"]), (0.5, 0.8))      # default to the birth values
+        self.assertEqual((em["mix_min"], em["displace"], em["spread"]), (0.3, 0.9, "contact"))
+        self.assertEqual(deep.tolerance["temperature"].optimum, 30)              # own curve
+        self.assertEqual(deep.tolerance["depth"].width_hi, 60)                   # inherited condition
+
+    def test_removed_keys_name_the_replacement(self):
+        with self.assertRaisesRegex(ConfigError, "removed.*birth_min"):
+            make_race("s", emerge={"parent": "p", "class": "deep_sea"})
+        with self.assertRaisesRegex(ConfigError, "removed.*birth_min"):
+            make_race("s", emerge={"parent": "p", "years": [1, 2]})
+        with self.assertRaisesRegex(ConfigError, "exposure.*removed"):
+            make_history(exposure={"hot": "above:T_mean:10"})
+        with self.assertRaisesRegex(ConfigError, "trigger"):
+            make_race("s", emerge={"parent": "p", "trigger": 0.1})               # ritual only
+        with self.assertRaisesRegex(ConfigError, "crisis_drop"):
+            make_race("s", emerge={"parent": "p", "crisis_drop": 0.9})           # ritual only
+        with self.assertRaisesRegex(ConfigError, "victims"):
+            make_race("s", emerge={"parent": "p", "mode": "ritual", "victims": "everyone"})
+        with self.assertRaisesRegex(ConfigError, "crisis_drop"):
+            make_race("s", emerge={"parent": "p", "mode": "ritual", "crisis_drop": 1.0})
+        self.assertEqual(make_race("s", emerge={"parent": "p", "birth_rate": 0.4}).emerge["return_rate"], 0.4)
+        for key, bad in [("displace", 1.5), ("dominance", -0.1), ("birth_shift", 1.2), ("convert", 0.0),
+                         ("mix_min", 1.1), ("birth_rate", 2.0), ("return_rate", -1.0), ("after", -5)]:
+            with self.assertRaisesRegex(ConfigError, key):
+                make_race("s", emerge={"parent": "p", key: bad})
+        for key, bad in [("ritual_rate", -0.1), ("ritual_cost", -1), ("ritual_converts", -1), ("trigger", 1.5),
+                         ("raid_size", 0), ("raid_km", -1), ("raid_rate", -1), ("ritual_power", 0),
+                         ("backlash_years", -1), ("crisis_years", 0)]:
+            with self.assertRaisesRegex(ConfigError, key):
+                make_race("s", emerge={"parent": "p", "mode": "ritual", key: bad})
+        with self.assertRaisesRegex(ConfigError, "mode"):
+            make_race("s", emerge={"parent": "p", "mode": "magic"})
+        with self.assertRaisesRegex(ConfigError, "birth_min"):
+            make_race("s", emerge={"parent": "p", "birth_min": 0.9, "birth_sure": 0.5})
+
+    def test_dumps_roundtrip(self):
+        h = make_history(regions={"z": {"kind": "box", "lat": [0, 1], "lon": [2, 3]}},
+                         event=[{"year": 0, "kind": "seed"}, {"year": 20, "kind": "cull", "share": 0.3,
+                                                                "by_race": {"a": 0.5}}])
+        doc = resolved(h, [make_race()], seed=5, engine="abc")
+        back = tomllib.loads(dumps_toml(doc))
+        self.assertEqual(back["seed"], 5)
+        self.assertEqual(back["history"]["events"][1]["by_race"], {"a": 0.5})
+        self.assertEqual(back["race"][0]["id"], "a")
+        self.assertEqual(back["history"]["regions"]["z"]["lat"], [0, 1])
diff --git a/tests/test_conflict.py b/tests/test_conflict.py
new file mode 100644
index 0000000..7fd7f6e
--- /dev/null
+++ b/tests/test_conflict.py
@@ -0,0 +1,77 @@
+import unittest
+
+import numpy as np
+
+from tests.helpers import make_race
+from worldhistory.conflict import apply_conflict, apply_curse, conflict
+from worldhistory.state import new_state
+
+
+def two(ca, cb):
+    return [make_race("a", conflict=ca), make_race("b", conflict=cb)]
+
+
+class ConflictTest(unittest.TestCase):
+    def test_internal_scales_with_crowding(self):
+        r = [make_race("a", conflict={"internal": 0.01})]
+        loss, _ = conflict(np.array([[100.0, 100.0]]), np.ones((1, 2)), np.array([[0.5, 1.0]]), r)
+        np.testing.assert_allclose(loss[0], [0.5, 1.0])
+
+    def test_no_border_losses_alone_or_within_family(self):
+        r = two({"border": 0.1}, {"border": 0.1})
+        r[1].family = "a"
+        loss, press = conflict(np.array([[50.0], [50.0]]), np.ones((2, 1)), np.zeros((2, 1)), r)
+        np.testing.assert_array_equal(loss, 0)
+        np.testing.assert_array_equal(press, 0)
+
+    def test_warlike_neighbour_hurts_more(self):
+        P = np.array([[50.0], [50.0]])
+        r = two({"border": 0.1, "aggression": 0.1}, {"border": 0.1, "aggression": 0.9})
+        loss, press = conflict(P, np.zeros((2, 1)), np.zeros((2, 1)), r)
+        self.assertGreater(loss[0, 0], loss[1, 0])
+        self.assertGreater(press[0, 0], press[1, 0])
+        # a: attacked = 0.9*1/1*0.5 = 0.45; attacking b outside b's core = 0 -> loss 0.1*0.45*50
+        self.assertAlmostEqual(float(loss[0, 0]), 0.1 * 0.45 * 50)
+
+    def test_defend_in_core_raises_attacker_losses(self):
+        P = np.array([[50.0], [50.0]])
+        base = two({"border": 0.1, "aggression": 0.8}, {"border": 0.1, "defend": 1.0})
+        hard = two({"border": 0.1, "aggression": 0.8}, {"border": 0.1, "defend": 3.0})
+        core = np.array([[0.0], [1.0]])                            # the cell is b's core habitat
+        l1, _ = conflict(P, core, np.zeros((2, 1)), base)
+        l3, _ = conflict(P, core, np.zeros((2, 1)), hard)
+        self.assertGreater(l3[0, 0], l1[0, 0])
+
+    def test_yield_turns_losses_into_moves(self):
+        st = new_state(1, 1, np.zeros((1, 6)), seed=0)
+        st.P[0, 0] = 100.0
+        apply_conflict(st, np.array([[10.0]]), [make_race(conflict={"yield": 0.8})])
+        self.assertAlmostEqual(float(st.P[0, 0]), 98.0)
+        self.assertAlmostEqual(float(st.push[0, 0]), 8.0)
+
+
+class DreadCurseTest(unittest.TestCase):
+    def test_dread_deters_attacks_and_repels(self):
+        P = np.array([[50.0], [50.0]])
+        plain = two({"border": 0.1}, {"border": 0.1, "aggression": 0.9})
+        feared = two({"border": 0.1, "dread": 0.8}, {"border": 0.1, "aggression": 0.9})
+        l0, p0 = conflict(P, np.zeros((2, 1)), np.zeros((2, 1)), plain)
+        l1, p1 = conflict(P, np.zeros((2, 1)), np.zeros((2, 1)), feared)
+        self.assertAlmostEqual(float(l1[0, 0]), 0.2 * float(l0[0, 0]))   # b's attacks on a × (1 − 0.8)
+        self.assertAlmostEqual(float(p1[1, 0]), float(p0[1, 0]) + 0.8 * 0.5)  # b is repelled by a's dread × a's share
+
+    def test_blame_falls_on_attackers_of_cursing_race(self):
+        P = np.array([[50.0], [50.0]])
+        r = two({"border": 0.1, "curse": 2.0}, {"border": 0.1, "aggression": 0.9})
+        blame = np.zeros((2, 1))
+        conflict(P, np.zeros((2, 1)), np.zeros((2, 1)), r, blame=blame)
+        # share of a's clan killed by b = border 0.1 × attacked 0.9·1·0.5 = 0.045; × curse 2
+        np.testing.assert_allclose(blame[:, 0], [0.0, 0.09])
+
+    def test_curse_kills_and_fades(self):
+        st = new_state(1, 2, np.zeros((1, 6)), seed=0)
+        st.P[0] = [100.0, 100.0]
+        apply_curse(st, np.array([[0.1, 0.0]]), decay=0.5)
+        np.testing.assert_allclose(st.P[0], [90.0, 100.0])
+        apply_curse(st, np.zeros((1, 2)), decay=0.5)
+        np.testing.assert_allclose(st.P[0], [85.5, 100.0])                 # curse halves: 5 % this step
diff --git a/tests/test_currents.py b/tests/test_currents.py
new file mode 100644
index 0000000..2b01a43
--- /dev/null
+++ b/tests/test_currents.py
@@ -0,0 +1,109 @@
+import unittest
+
+import numpy as np
+
+from tests.helpers import line_world, make_history, make_race
+from worldhistory.habitat import environment, suitability
+from worldhistory.migration import drift_and_bud, long_jumps
+from worldhistory.predicates import evaluate
+from worldhistory.state import new_state
+
+MP = make_history().migration
+
+
+def east(w, speed=1.0):
+    """Uniform current along the line (the line runs along the equator, eastward = +lon)."""
+    lon = np.radians(w.lon)
+    return speed * np.stack([-np.sin(lon), np.cos(lon), np.zeros(w.n)], 1)
+
+
+def sea_world(n, spacing_km=100.0, current=True, **over):
+    w = line_world(n, spacing_km=spacing_km, ocean=True, elevation_m=-3000.0, **over)
+    if current:
+        w.base["current"] = east(w)
+    return w
+
+
+def swimmer(bias):
+    return make_race(demography={"founder": 2.0, "allee": 1.0, "jump_path": "sea", "current_bias": bias,
+                                 "mobility": 0.2},
+                     habitat={"realm": "sea", "terms": [{"p": "sea", "w": 1.0}]})
+
+
+class SailCostTest(unittest.TestCase):
+    def test_downstream_cheaper_than_still_cheaper_than_upstream(self):
+        w = sea_world(3)
+        c = w.sail_cost(5.0)
+        k_east = list(w.nb[1]).index(2)
+        k_west = list(w.nb[1]).index(0)
+        edge_h = 100.0 * 1000.0 / 3600.0                                # 100 km in hours per (m/s)
+        self.assertAlmostEqual(c[1, k_east], edge_h / 6.0, delta=0.01)  # 5 + 1 m/s
+        self.assertAlmostEqual(c[1, k_west], edge_h / 4.0, delta=0.01)  # 5 − 1 m/s
+        still = sea_world(3, current=False).sail_cost(5.0)
+        self.assertAlmostEqual(still[1, k_east], edge_h / 5.0, delta=0.01)
+        self.assertTrue(np.isinf(c[0, list(w.nb[0]).index(-1)]))        # no neighbour: no edge
+
+    def test_land_edges_are_impassable(self):
+        w = line_world(3, ocean=[True, True, False])
+        c = w.sail_cost(5.0)
+        self.assertTrue(np.isinf(c[1, list(w.nb[1]).index(2)]))
+
+
+class FieldsTest(unittest.TestCase):
+    def test_missing_current_is_still_water(self):
+        w = sea_world(3, current=False)
+        self.assertEqual(w.current.shape, (3, 3))
+        self.assertTrue(np.all(w.current == 0))
+
+    def test_productivity_predicate(self):
+        w = sea_world(3, productivity=np.array([0.0, 0.5, 1.0]))
+        np.testing.assert_allclose(evaluate("productivity", w), [0.0, 0.5, 1.0])
+        np.testing.assert_allclose(evaluate("productivity", sea_world(3)), [0.0, 0.0, 0.0])   # absent → 0
+
+
+class JumpBiasTest(unittest.TestCase):
+    def landings(self, bias, current=True):
+        w = sea_world(400, spacing_km=10.0, current=current)
+        r = swimmer(bias)
+        env, suit = environment(w), suitability(w, r)[None]
+        east_n = west_n = 0
+        for seed in range(150):
+            st = new_state(1, w.n, np.zeros((1, 6)), seed=0)
+            st.P[0, 200] = 1e6
+            long_jumps(st, w, [r], np.full((1, w.n), 1e6), np.ones((1, w.n), bool), np.random.default_rng(seed),
+                       {**MP, "jump_rate": 1.0}, np.ones((1, w.n)), env=env, suit=suit, f_min=0.0)
+            hit = np.flatnonzero(st.P[0] > 0)
+            east_n += int(np.sum(hit > 200))
+            west_n += int(np.sum(hit < 200))
+        return east_n, west_n
+
+    def test_downstream_jumps_dominate_with_bias(self):
+        e, w = self.landings(3.0)
+        self.assertGreater(e, 3 * max(w, 1))
+
+    def test_no_bias_is_symmetric(self):
+        e, w = self.landings(0.0)
+        self.assertLess(abs(e - w), 0.35 * (e + w) + 5)
+
+    def test_no_current_field_ignores_bias(self):
+        e, w = self.landings(3.0, current=False)
+        self.assertLess(abs(e - w), 0.35 * (e + w) + 5)
+
+
+class DriftBiasTest(unittest.TestCase):
+    def drift(self, bias):
+        w = sea_world(3)
+        r = swimmer(bias)
+        st = new_state(1, 3, np.zeros((1, 6)), seed=0)
+        st.P[0] = [100.0, 1000.0, 100.0]                          # neighbours settled: drift, not founding
+        drift_and_bud(st, w, [r], np.full((1, 3), 1e5), np.ones((1, 3), bool), np.zeros((1, 3)),
+                      np.zeros((1, 3)), np.ones((1, 3)), np.random.default_rng(0), MP)
+        return st.P[0]
+
+    def test_drift_follows_the_current(self):
+        P = self.drift(2.0)
+        self.assertGreater(P[2] - 100.0, 3.0 * max(P[0] - 100.0, 1e-9))   # the downstream neighbour gains most
+
+    def test_no_bias_drifts_evenly(self):
+        P = self.drift(0.0)
+        self.assertAlmostEqual(P[0], P[2], delta=1e-6 * P.sum())
diff --git a/tests/test_demography.py b/tests/test_demography.py
new file mode 100644
index 0000000..26a5f04
--- /dev/null
+++ b/tests/test_demography.py
@@ -0,0 +1,123 @@
+import unittest
+
+import numpy as np
+
+from tests.helpers import line_world, make_race
+from worldhistory.demography import (allee, capacity, crowding, grow, harvest_shocks, overlap_matrix,
+                                     stochastic_round, veins)
+
+
+class DemographyTest(unittest.TestCase):
+    def test_capacity_formula(self):
+        w = line_world(2)
+        r = make_race(demography={"density": 2.0})
+        K = capacity(w, r, np.array([1.0, 0.5]), np.ones(2), np.ones(2), np.array([11.0, 11.0]), np.ones(2),
+                     np.ones(2))
+        # density * area * q * (1 + gain q²): 2*10000*1*12, 2*10000*0.5*(1+11*0.25)
+        np.testing.assert_allclose(K, [240000, 37500])
+
+    def test_logistic_reaches_K(self):
+        P = np.array([[10.0]])
+        K = np.array([[1000.0]])
+        r = [make_race(demography={"growth": 0.5})]
+        for _ in range(200):
+            grow(P, K, crowding(P, K, np.zeros((1, 1)), r, 2.0), np.ones((1, 1)), r, np.ones((1, 1)))
+        self.assertAlmostEqual(float(P[0, 0]), 1000.0, places=3)
+
+    def test_growth_rate_low_gear(self):
+        P = np.array([[10.0, 10.0]])
+        K = np.array([[1e9, 1e9]])
+        r = [make_race(demography={"growth": 0.1})]
+        grow(P, K, crowding(P, K, np.zeros((1, 1)), r, 2.0), np.array([[1.0, 0.0]]), r, np.ones((1, 2)))
+        np.testing.assert_allclose(P[0], [11.0, 10.3], rtol=1e-6)      # r*(0.3+0.7q)
+
+    def test_grow_returns_births(self):
+        r = [make_race()]
+        P = np.array([[100.0, 100.0]])
+        b = grow(P, np.array([[1e6, 50.0]]), np.array([[0.0, 2.0]]), np.ones((1, 2)), r, np.ones((1, 2)))
+        # gross births (generational turnover), not net growth: a crowded, shrinking group still bears children
+        np.testing.assert_allclose(b[0], [0.1 * 100, 0.1 * 100])
+
+    def test_allee(self):
+        self.assertLess(allee(np.array(10.0), 75.0), 0)
+        self.assertGreater(allee(np.array(100.0), 75.0), 0)
+        self.assertAlmostEqual(float(allee(np.array(1.0), 1000.0)), -0.3)
+        P = np.array([[10.0]])
+        r = [make_race(demography={"growth": 0.1, "allee": 75.0})]
+        grow(P, np.array([[1e6]]), np.zeros((1, 1)), np.ones((1, 1)), r, np.ones((1, 1)))
+        self.assertLess(float(P[0, 0]), 10.0)
+
+    def test_no_capacity_halves(self):
+        P = np.array([[10.0]])
+        grow(P, np.array([[0.0]]), np.zeros((1, 1)), np.ones((1, 1)), [make_race()], np.ones((1, 1)))
+        self.assertAlmostEqual(float(P[0, 0]), 5.0)
+
+    def test_competition(self):
+        P = np.array([[100.0], [100.0]])
+        K = np.array([[1000.0], [1000.0]])
+        rs = [make_race("a"), make_race("b")]
+        c0 = crowding(P, K, np.zeros((2, 2)), rs, 1.0)
+        c1 = crowding(P, K, np.array([[0, 1.0], [1.0, 0]]), rs, 1.0)
+        np.testing.assert_allclose(c0[:, 0], [0.1, 0.1])
+        np.testing.assert_allclose(c1[:, 0], [0.2, 0.2])
+
+    def test_tolerated_share_sharpens(self):
+        P = np.array([[100.0], [100.0]])
+        K = np.array([[1000.0], [1000.0]])
+        rs = [make_race("a", temperament={"tolerated_share": 0.1}), make_race("b", temperament={"tolerated_share": 0.9})]
+        c = crowding(P, K, np.array([[0, 1.0], [1.0, 0]]), rs, 2.0)
+        np.testing.assert_allclose(c[:, 0], [0.3, 0.2])               # a: others' share 0.5 > 0.1 → ×2
+
+    def test_overlap_matrix(self):
+        S = np.array([[1.0, 0, 0], [1.0, 0, 0], [0, 0, 1.0]])
+        calm = {"xenophobia": 0.0, "pace": 0.5, "structure": 0.5, "warlike": 0.0}
+        rs = [make_race("a", temperament=calm), make_race("b", temperament=calm), make_race("c", temperament=calm)]
+        a = overlap_matrix(S, rs)
+        self.assertAlmostEqual(a[0, 1], 0.5)          # sim 1 × friction 0.5
+        self.assertAlmostEqual(a[0, 2], 0.0)          # disjoint habitats
+        rs[0].temperament = {**calm, "xenophobia": 1.0}
+        a = overlap_matrix(S, rs)
+        self.assertAlmostEqual(a[0, 1], 1.0)          # a minds b more than b minds a
+        self.assertAlmostEqual(a[1, 0], 0.5)
+        rs[1].overlap = {"a": 0.05}
+        self.assertAlmostEqual(overlap_matrix(S, rs)[1, 0], 0.05)
+
+    def test_shocks_mean_one_and_scale(self):
+        w = line_world(20000)
+        e = np.random.default_rng(1).normal(size=w.n)
+        m = harvest_shocks(w, np.ones(w.n), np.ones(w.n), e)
+        self.assertAlmostEqual(float(m.mean()), 1.0, places=2)
+        np.testing.assert_array_equal(harvest_shocks(w, np.zeros(w.n), np.ones(w.n), e), 1.0)
+
+    def test_veins_need_viable_community(self):
+        np.testing.assert_allclose(veins(np.array([50.0, 100.0]), np.array([1000.0, 1000.0]), 0.01, 75.0),
+                                   [50.0, 109.0])
+
+    def test_stochastic_round_expectation(self):
+        P = np.full((1, 100000), 2.3)
+        stochastic_round(P, np.random.default_rng(0))
+        self.assertTrue(set(np.unique(P)) <= {2.0, 3.0})
+        self.assertAlmostEqual(float(P.mean()), 2.3, places=2)
+        Q = np.array([[7.5, 0.0]])
+        stochastic_round(Q, np.random.default_rng(0))
+        np.testing.assert_array_equal(Q, [[7.5, 0.0]])
+
+
+class RangeTest(unittest.TestCase):
+    def test_wide_range_crowding(self):
+        # a clan drawing on a wide range: crowding = people over capacity summed across the neighbourhood
+        from tests.helpers import line_world
+        w = line_world(3)
+        P, K = np.array([[300.0, 0.0, 0.0]]), np.array([[100.0, 100.0, 100.0]])
+        narrow = crowding(P, K, np.zeros((1, 1)), [make_race()], 2.0)
+        wide = crowding(P, K, np.zeros((1, 1)), [make_race(demography={"range": 1})], 2.0, world=w)
+        self.assertAlmostEqual(float(narrow[0, 0]), 3.0)
+        self.assertAlmostEqual(float(wide[0, 0]), 1.5)                 # 300 over cells 0+1 (200)
+
+    def test_world_range_crowding(self):
+        # range -1: crowding is the whole people's numbers against its whole capacity (few, widely scattered clans)
+        from tests.helpers import line_world
+        w = line_world(3)
+        P, K = np.array([[300.0, 0.0, 0.0]]), np.array([[100.0, 100.0, 0.0]])
+        c = crowding(P, K, np.zeros((1, 1)), [make_race(demography={"range": -1})], 2.0, world=w)
+        np.testing.assert_allclose(c[0], [1.5, 1.5, 0.0])                # 300 / 200 where it can live
diff --git a/tests/test_engine.py b/tests/test_engine.py
new file mode 100644
index 0000000..b0178e2
--- /dev/null
+++ b/tests/test_engine.py
@@ -0,0 +1,128 @@
+import os
+import tempfile
+import tomllib
+import unittest
+from pathlib import Path
+
+import numpy as np
+
+from tests.helpers import fields, globe_world, line_world, make_history, make_race
+from worldhistory.config import CONDITIONS, link, load_config
+from worldhistory.engine import Engine
+from worldhistory.nudges import nudge_multipliers
+from worldhistory.regions import Regions
+from worldhistory.snapshot import load_snapshot
+
+
+def small_setup(seed=1):
+    w = globe_world(1)
+    w.base["ocean"] = w.lat < -40
+    w.eras["late"] = fields(w.n, ocean=w.lat < -30, gravity_g=np.where(w.lon > 90, 0.35, 1.0))
+    h = make_history(run={"years": 60, "step": 10, "snapshot_every": 20, "seed": seed},
+                     world={"eras": ["late"]},
+                     regions={"zone": {"kind": "changed", "era": "late", "fields": ["gravity_g"], "land": "base"}},
+                     event=[{"year": 0, "kind": "seed"},
+                            {"year": 20, "kind": "cull", "region": "zone", "share": 0.6},
+                            {"year": 20, "kind": "era_switch", "era": "late"}],
+                     hazard=[{"kind": "roaming", "start": 20, "home_region": "zone", "local": 3, "wanderers": 1}],
+                     stats={"regions": ["zone"]})
+    races = link(h, [make_race("a", seed={"clusters": 10, "heads": 500},
+                               demography={"density": 2.0, "growth": 0.1, "mobility": 0.1}),
+                     make_race("b", seed={"clusters": 10, "heads": 500}, conflict={"border": 0.01}),
+                     make_race("c", emerge={"parent": "a", "birth_min": 0.0, "birth_sure": 0.0, "birth_rate": 1.0,
+                                            "return_min": 1.0, "return_sure": 1.0})])
+    return w, h, races
+
+
+class EngineTest(unittest.TestCase):
+    def test_run_writes_outputs(self):
+        w, h, races = small_setup()
+        with tempfile.TemporaryDirectory() as d:
+            stats = Engine(w, h, races).run(d)
+            files = sorted(p.name for p in (Path(d) / "snap").iterdir())
+            self.assertEqual(files, ["y0000.npz", "y0020.npz", "y0040.npz", "y0060.npz"])
+            run = tomllib.loads((Path(d) / "run.toml").read_text())
+            self.assertEqual(run["seed"], 1)
+            snap = load_snapshot(Path(d) / "snap" / "y0060.npz", w.n)
+            self.assertEqual(snap["races"], ["a", "b", "c"])
+            self.assertTrue((Path(d) / "stats.json").exists())
+        self.assertEqual([s["year"] for s in stats["steps"]], [0, 10, 20, 30, 40, 50, 60])
+        self.assertEqual(stats["events"][0]["after"]["a"], 5000.0)     # 10 clusters × 500 at the gifting
+        kinds = [e["kind"] for e in stats["events"]]
+        self.assertEqual(kinds, ["seed", "cull", "era_switch"])
+        self.assertIn("zone", stats["steps"][-1]["regions"])
+        self.assertEqual(len(stats["emergence"]), 1)                  # exposure 10 yrs → c emerges once
+        self.assertGreater(stats["steps"][-1]["pop"]["c"], 0)
+        self.assertEqual(stats["steps"][-1]["family"]["a"],
+                         stats["steps"][-1]["pop"]["a"] + stats["steps"][-1]["pop"]["c"])
+
+    def test_no_birth_where_native_curves_do_not_fit(self):
+        w, h, _ = small_setup()
+        races = link(h, [make_race("a", seed={"clusters": 10, "heads": 500},
+                                   demography={"density": 2.0, "growth": 0.1, "mobility": 0.1}),
+                         make_race("c", emerge={"parent": "a", "birth_min": 0.0, "birth_sure": 0.0,
+                                                "birth_rate": 1.0},
+                                   tolerance={"temperature": {"optimum": -80, "width": 1, "lo": -90, "hi": 60}})])
+        with tempfile.TemporaryDirectory() as d:
+            stats = Engine(w, h, races).run(d)
+        self.assertEqual(stats["emergence"], [])
+
+    def test_progress_and_returns_in_stats(self):
+        w, h, races = small_setup()
+        with tempfile.TemporaryDirectory() as d:
+            stats = Engine(w, h, races).run(d)
+        self.assertIn("a>c", stats["steps"][-1]["progress"])
+        self.assertIsInstance(stats["returns"], list)
+        self.assertIsInstance(stats["raids"], list)
+        self.assertIn("a", stats["emergence"][0])
+
+    def test_deterministic(self):
+        a = Engine(*small_setup(1))
+        b = Engine(*small_setup(1))
+        c = Engine(*small_setup(2))
+        for e in (a, b, c):
+            for _ in range(5):
+                e.step()
+        np.testing.assert_array_equal(a.state.P, b.state.P)
+        self.assertFalse(np.array_equal(a.state.P, c.state.P))
+
+    def test_depth_creep_rate(self):
+        n = 16
+        w = line_world(n, ocean=[False] + [True] * (n - 1), elevation_m=[10] + [-100.0 * i for i in range(1, n)])
+        tol = {"depth": {"optimum": 0, "width": 30, "lo": 0, "hi": 5000, "rate": 2.0}}
+        h = make_history(run={"years": 600, "step": 10, "seed": 1}, event=[])
+        r = make_race("d", habitat={"realm": "both", "terms": [{"p": "land", "w": 1}, {"p": "sea", "w": 1}]},
+                      demography={"density": 5.0, "growth": 0.3, "founder": 2.0, "mobility": 0.2}, tolerance=tol)
+        e = Engine(w, h, link(h, [r]))
+        e.state.P[0, 0] = 1000.0
+        front = []
+        for _ in range(60):
+            e.step()
+            occ = np.flatnonzero(e.state.P[0] >= 1)
+            front.append(100.0 * occ.max())
+        self.assertEqual(front, sorted(front))                        # only deeper with time
+        self.assertGreater(front[-1], 0.4 * 2.0 * 600)                # ≈ rate × years (1,200 m)
+        self.assertLess(front[-1], 1.3 * 2.0 * 600)
+
+    def test_nudges(self):
+        w = line_world(3)
+        g, c, x, a = nudge_multipliers([{"race": "a", "region": "all", "years": [0, 100], "growth": 1.1,
+                                         "capacity": 1.2, "expansion": 1.3, "adapt": 2.0}],
+                                       [make_race("a")], Regions(w, {}), 50, 3)
+        np.testing.assert_allclose(g, 1.1)
+        np.testing.assert_allclose(c, 1.2)
+        np.testing.assert_allclose(a, 2.0)
+        g, c, x, a = nudge_multipliers([{"race": "a", "region": "all", "years": [0, 10], "growth": 1.1, "capacity": 1,
+                                         "expansion": 1, "adapt": 1}], [make_race("a")], Regions(w, {}), 50, 3)
+        np.testing.assert_allclose(g, 1.0)
+
+
+@unittest.skipUnless(os.environ.get("WORLDHISTORY_SMOKE_WORLD"), "set WORLDHISTORY_SMOKE_WORLD and _CONFIG")
+class SmokeTest(unittest.TestCase):
+    def test_100_years_on_a_real_world(self):
+        from worldhistory.world import load_world
+        h, races = load_config(os.environ["WORLDHISTORY_SMOKE_CONFIG"])
+        w = load_world(os.environ["WORLDHISTORY_SMOKE_WORLD"], eras=h.eras, extra_fields=h.world_fields)
+        with tempfile.TemporaryDirectory() as d:
+            stats = Engine(w, h, races).run(d, years=100)
+        self.assertGreater(sum(stats["steps"][-1]["pop"].values()), 0)
diff --git a/tests/test_events.py b/tests/test_events.py
new file mode 100644
index 0000000..89ef9e5
--- /dev/null
+++ b/tests/test_events.py
@@ -0,0 +1,85 @@
+import unittest
+
+import numpy as np
+
+from tests.helpers import fields, line_world, make_history, make_race
+from worldhistory.events import apply_event, cull, era_switch, seed
+from worldhistory.regions import Regions
+from worldhistory.state import new_state
+
+
+class EventTest(unittest.TestCase):
+    def test_seed_totals_and_fringe(self):
+        w = line_world(200)
+        r = make_race(seed={"clusters": 20, "heads": 100})
+        q = np.linspace(0, 1, 200)[None]
+        st = new_state(1, 200, np.zeros((1, 6)), seed=0)
+        placed = seed(st, w, [r], q, np.ones((1, 200), bool), 0.25, np.random.default_rng(0))
+        self.assertAlmostEqual(float(st.P.sum()), 2000.0)
+        cells = placed["a"]
+        fringe_cells = [c for c in cells if 0.05 < q[0, c] < 0.3]
+        self.assertGreaterEqual(len(fringe_cells), 5)                     # 25 % of 20 on the fringe
+        self.assertEqual(float(st.P[0, 0]), 0.0)                          # q = 0: never
+
+    def test_seed_realm_land_skips_sea(self):
+        ocean = np.arange(200) >= 100                                     # right half is sea, and the best q
+        w = line_world(200, ocean=ocean)
+        r = make_race(seed={"clusters": 20, "heads": 100, "realm": "land"}, habitat={"realm": "both"})
+        q = np.linspace(0.1, 1, 200)[None]
+        st = new_state(1, 200, np.zeros((1, 6)), seed=0)
+        placed = seed(st, w, [r], q, np.ones((1, 200), bool), 0.25, np.random.default_rng(0))
+        self.assertTrue(placed["a"])
+        self.assertTrue(all(c < 100 for c in placed["a"]))
+
+    def test_cull_exact_and_by_race(self):
+        rs = [make_race("a"), make_race("b")]
+        st = new_state(2, 3, np.zeros((2, 6)), seed=0)
+        st.P[:] = 100.0
+        cull(st, rs, np.array([True, True, False]), 0.6, {"b": 0.8})
+        np.testing.assert_allclose(st.P, [[40, 40, 100], [20, 20, 100]])
+
+    def test_cull_noise(self):
+        w = line_world(500)
+        st = new_state(1, 500, np.zeros((1, 6)), seed=0)
+        st.P[:] = 100.0
+        cull(st, [make_race()], np.ones(500, bool), 0.3, {}, noise=w.noise(1), amp=0.15)
+        self.assertAlmostEqual(float(st.P.mean()), 70.0, delta=2.0)
+        self.assertGreater(float(st.P.std()), 5.0)
+
+    def test_era_switch_drowns_land_races(self):
+        w = line_world(3)
+        w.eras["late"] = fields(3, ocean=[False, True, True])
+        land, sea = make_race("l"), make_race("s", habitat={"realm": "sea", "terms": [{"p": "sea", "w": 1}]})
+        st = new_state(2, 3, np.zeros((2, 6)), seed=0)
+        st.P[0] = 10.0
+        st.P[1, 2] = 5.0
+        era_switch(st, w, [land, sea], "late")
+        np.testing.assert_array_equal(st.P[0], [10, 0, 0])
+        np.testing.assert_array_equal(st.P[1], [0, 0, 5])
+        self.assertEqual(w.era, "late")
+
+    def test_targeted_seed_places_one_race_in_region(self):
+        w = line_world(200)                                               # lon 0 .. ~178°, 0.9° apart
+        h = make_history(regions={"west": {"kind": "box", "lat": [-1, 1], "lon": [-1, 45]}},
+                         event=[{"year": 0, "kind": "seed", "race": "b", "region": "west", "clusters": 3,
+                                 "heads": 50}])
+        rs = [make_race("a", seed={"clusters": 5, "heads": 100}), make_race("b", seed={"clusters": 5, "heads": 100})]
+        q = np.ones((2, 200)) * 0.5
+        st = new_state(2, 200, np.zeros((2, 6)), seed=0)
+        log = apply_event(st, w, rs, h.events[0], Regions(w, h.regions), q, np.ones((2, 200), bool), h, 0)
+        self.assertEqual(float(st.P[0].sum()), 0.0)                       # race a: untouched
+        self.assertAlmostEqual(float(st.P[1].sum()), 150.0)               # 3 clusters × 50
+        west = np.flatnonzero(w.lon <= 45)
+        self.assertTrue(set(np.flatnonzero(st.P[1])) <= set(west))
+        self.assertEqual(len(log["placed"]["b"]), 3)
+
+    def test_apply_event_log(self):
+        w = line_world(2)
+        h = make_history(regions={"west": {"kind": "box", "lat": [-1, 1], "lon": [-1, 0.5]}},
+                         event=[{"year": 20, "kind": "die_off", "region": "west", "share": 0.9}])
+        st = new_state(1, 2, np.zeros((1, 6)), seed=0)
+        st.P[:] = 100.0
+        log = apply_event(st, w, [make_race()], h.events[0], Regions(w, h.regions), None, None, h, 0)
+        np.testing.assert_allclose(st.P[0], [10, 100])
+        self.assertEqual(log["before"]["a"], 200.0)
+        self.assertEqual(log["after"]["a"], 110.0)
diff --git a/tests/test_fastpaths.py b/tests/test_fastpaths.py
new file mode 100644
index 0000000..893be9b
--- /dev/null
+++ b/tests/test_fastpaths.py
@@ -0,0 +1,314 @@
+"""The speed-ups that compute only where values can change give the same floats, bit for bit, as the full-array
+versions they replaced (kept here as oracles)."""
+import unittest
+
+import numpy as np
+
+from worldhistory import state as S
+from worldhistory.state import ATTRS, new_state
+
+
+def move_full(st, r, src, dst, amt):
+    """The original move: every cell updated."""
+    src, dst, amt = np.asarray(src, np.int64), np.asarray(dst, np.int64), np.asarray(amt, float)
+    keep = amt > 0
+    src, dst, amt = src[keep], dst[keep], amt[keep]
+    if not len(amt):
+        return
+    n = st.P.shape[1]
+    P = st.P[r]
+    out = np.bincount(src, amt, n)
+    scale = np.where(out > P, P / np.maximum(out, 1e-12), 1.0)
+    amt = amt * scale[src]
+    stay = np.maximum(P - np.bincount(src, amt, n), 0.0)
+    inflow = np.bincount(dst, amt, n)
+    tot = stay + inflow
+    for name in ATTRS:
+        A = getattr(st, name)[r]
+        if A.shape[0] == 0:
+            continue
+        Sm = np.stack([np.bincount(dst, amt * A[x, src], n) for x in range(A.shape[0])])
+        A[:] = np.where(tot > 0, (stay * A + Sm) / np.maximum(tot, 1e-12), A)
+    st.P[r] = tot
+
+
+def random_state(n=4000, seed=0):
+    rng = np.random.default_rng(seed)
+    st = new_state(2, n, np.full((2, 6), 10.0, np.float32), seed)
+    st.P[:] = np.where(rng.random((2, n)) < 0.4, rng.random((2, n)) * 1e4, 0.0)
+    st.P[0, :20] = rng.random(20) * 1e-13               # thin groups (below the 1e-12 floor)
+    st.O[:] = (rng.normal(10, 5, st.O.shape)).astype(np.float32)
+    st.T[:] = rng.random(st.T.shape).astype(np.float32)
+    st.C[:, 0] = np.where(rng.random((2, n)) < 0.1, rng.random((2, n)) * 0.3, 0.0)   # curses: float64, nonzero
+    return st, rng
+
+
+class MoveTest(unittest.TestCase):
+    def test_move_matches_full_update(self):
+        for seed in range(4):
+            a, rng = random_state(seed=seed)
+            b, _ = random_state(seed=seed)
+            n = a.P.shape[1]
+            for _ in range(5):
+                k = int(rng.integers(1, 3000))
+                src, dst = rng.integers(0, n, k), rng.integers(0, n, k)
+                amt = rng.random(k) * 5e3 * (rng.random(k) < 0.9)
+                S.move(a, 0, src, dst, amt)
+                move_full(b, 0, src, dst, amt)
+            for name in ("P", *ATTRS):
+                with self.subTest(seed=seed, attr=name):
+                    self.assertTrue(np.array_equal(getattr(a, name), getattr(b, name)))
+
+
+class MoveSumOrderTest(unittest.TestCase):
+    def test_many_moves_into_few_cells_sum_in_move_order(self):
+        """Thousands of moves into a few cells, float64 curses: any other summation order shows in the last bits."""
+        for seed in range(3):
+            a, rng = random_state(seed=seed)
+            b, _ = random_state(seed=seed)
+            n = a.P.shape[1]
+            c = rng.random(n) * 0.4
+            for st in (a, b):
+                st.C[0, 0] = c
+            k = 20000
+            src, dst = rng.integers(0, n, k), rng.integers(0, 50, k)
+            amt = np.exp(rng.normal(0, 3, k))
+            S.move(a, 0, src, dst, amt)
+            move_full(b, 0, src, dst, amt)
+            for name in ("P", *ATTRS):
+                with self.subTest(seed=seed, attr=name):
+                    self.assertTrue(np.array_equal(getattr(a, name), getattr(b, name)))
+
+
+class ConflictTest(unittest.TestCase):
+    def test_conflict_matches_full_arrays_many_races(self):
+        """≥ 8 races: numpy sums a column-major block pairwise — the subset must still add race by race."""
+        from tests.helpers import make_race
+        from worldhistory.conflict import conflict
+        rng = np.random.default_rng(9)
+        races = [make_race(id=f"r{k}", conflict={"aggression": rng.random(), "power": 0.5 + rng.random(),
+                                                "dread": rng.random() * 0.5, "defend": 0.5, "border": 0.1,
+                                                "curse": 0.3 * (k % 3 == 0)}, family=f"f{k % 5}") for k in range(9)]
+        n = 20000
+        P = np.where(rng.random((9, n)) < 0.25, np.exp(rng.normal(5, 3, (9, n))), 0.0)
+        q, crowd = rng.random((9, n)), rng.random((9, n)) * 3
+        blame, want_blame = np.zeros_like(P), np.zeros_like(P)
+        loss, press = conflict(P, q, crowd, races, blame=blame)
+        want_loss, want_press = conflict_full(P, q, crowd, races, blame=want_blame)
+        self.assertTrue(np.array_equal(loss, want_loss))
+        self.assertTrue(np.array_equal(press, want_press))
+        self.assertTrue(np.array_equal(blame, want_blame))
+
+    def test_conflict_matches_full_arrays(self):
+        from tests.helpers import make_race
+        from worldhistory.conflict import conflict
+        rng = np.random.default_rng(5)
+        races = [make_race(id=k, conflict={"aggression": a, "power": p, "dread": d, "defend": 0.5, "border": 0.1,
+                                          "curse": c}, family=f)
+                 for k, a, p, d, c, f in (("a", 0.4, 1.0, 0.2, 0.0, "x"), ("b", 0.9, 1.5, 0.0, 0.5, "y"),
+                                          ("c", 0.2, 0.7, 0.6, 0.0, "z"))]
+        n = 5000
+        P = np.where(rng.random((3, n)) < 0.3, rng.random((3, n)) * 100, 0.0)
+        P[2] = 0.0                                       # an absent race
+        q, crowd = rng.random((3, n)), rng.random((3, n)) * 3
+        blame = np.zeros_like(P)
+        loss, press = conflict(P, q, crowd, races, blame=blame)
+        want_blame = np.zeros_like(P)
+        want_loss, want_press = conflict_full(P, q, crowd, races, blame=want_blame)
+        self.assertTrue(np.array_equal(loss, want_loss))
+        self.assertTrue(np.array_equal(press, want_press))
+        self.assertTrue(np.array_equal(blame, want_blame))
+
+
+def conflict_full(P, q_eff, crowd, races, core_q=0.6, blame=None):
+    R = len(races)
+    tot = P.sum(0)
+    share = np.where(tot > 0, P / np.maximum(tot, 1e-12), 0.0)
+    fam = [r.family for r in races]
+    C = [r.conflict for r in races]
+    loss, press = np.zeros_like(P), np.zeros_like(P)
+    for i in range(R):
+        ci = C[i]
+        loss[i] += ci["internal"] * np.minimum(crowd[i], 2.0) * P[i]
+        for j in range(R):
+            if fam[j] == fam[i]:
+                continue
+            cj = C[j]
+            attacked = cj["aggression"] * (1 - ci["dread"]) * cj["power"] / ci["power"] * share[j]
+            attacking = ci["aggression"] * cj["defend"] * (q_eff[j] >= core_q) * cj["power"] / ci["power"] * share[j]
+            press[i] += attacked + cj["dread"] * share[j]
+            loss[i] += ci["border"] * (attacked + attacking) * P[i]
+            if blame is not None and ci["curse"] > 0:
+                blame[j] += ci["curse"] * ci["border"] * attacked
+    return np.minimum(loss, 0.9 * P), press
+
+
+def prospective_full(st, world):
+    """The original prospective: neighbour sums over every cell."""
+    for r in range(st.P.shape[0]):
+        pos = st.P[r] > 0
+        if not pos.any():
+            continue
+        P = np.where(pos, st.P[r], 0.0)
+        wsum = world.nb_sum(P)
+        empty = np.flatnonzero((st.P[r] <= 0) & (wsum > 0))
+        if not len(empty):
+            continue
+        st.O[r][:, empty] = world.nb_sum(P * st.O[r])[:, empty] / np.maximum(wsum[empty], 1e-12)
+
+
+def step_tech_full(st, world, races, tp, steps=1.0):
+    """The original step_tech: neighbour sums over every cell."""
+    from worldhistory.config import DOMAINS
+    from worldhistory.tech import neigh_pop
+    for r, race in enumerate(races):
+        P = st.P[r]
+        occ = P >= 1
+        if not occ.any():
+            continue
+        oc = np.flatnonzero(occ)
+        N = neigh_pop(world, P, tp["passes"])[oc]
+        s = N / (N + tp["n_half"])
+        T = st.T[r]
+        for d, name in enumerate(DOMAINS):
+            Td = T[d, oc]
+            gain = tp["rate"] * steps * race.tech.get(name, 1.0) * s * (1 - Td)
+            loss = np.where(N < tp["loss_below"], tp["loss_rate"] * steps * Td, 0.0)
+            T[d, oc] = np.clip(Td + gain - loss, 0, 1)
+        if tp["diffuse"] > 0:
+            PT = P * T
+            Pc, Tc = P[oc], T[:, oc]
+            M = (PT[:, oc] + world.nb_sum(PT)[:, oc]) / np.maximum(Pc + world.nb_sum(P)[oc], 1e-12)
+            T[:, oc] = Tc + tp["diffuse"] * np.clip(M - Tc, 0, None)
+
+
+def sparse_state(n, seed, frac):
+    """Groups on a few patches (as in a run: most cells empty), thin and fractional groups among them."""
+    st, rng = random_state(n, seed)
+    st.P[:] = np.where(rng.random((2, n)) < frac, rng.random((2, n)) * 3e3, 0.0)
+    st.P[0, :5] = rng.random(5) * 0.5
+    return st
+
+
+class NeighbourhoodTest(unittest.TestCase):
+    def setUp(self):
+        from tests.helpers import globe_world
+        self.w = globe_world(2)                                         # 5882 cells
+
+    def test_adjacency_is_symmetric(self):
+        a = self.w._adj
+        self.assertEqual((a != a.T).nnz, 0)
+
+    def test_nb_sum_at_matches_full_sum(self):
+        rng = np.random.default_rng(3)
+        x = rng.normal(size=(3, self.w.n))
+        rows = np.sort(rng.choice(self.w.n, 400, replace=False))
+        self.assertTrue(np.array_equal(self.w.nb_sum_at(x, rows), self.w.nb_sum(x)[:, rows]))
+        self.assertTrue(np.array_equal(self.w.nb_sum_at(x[0], rows), self.w.nb_sum(x[0])[rows]))
+
+    def test_prospective_matches_full(self):
+        from worldhistory.adaptation import prospective
+        for seed, frac in ((0, 0.02), (1, 0.2), (2, 0.0005), (3, 0.9)):
+            a, b = sparse_state(self.w.n, seed, frac), sparse_state(self.w.n, seed, frac)
+            prospective(a, self.w)
+            prospective_full(b, self.w)
+            with self.subTest(seed=seed):
+                self.assertTrue(np.array_equal(a.O, b.O))
+                self.assertFalse(np.array_equal(a.O, sparse_state(self.w.n, seed, frac).O))   # it did change O
+
+    def test_step_tech_matches_full(self):
+        from tests.helpers import make_history, make_race
+        from worldhistory.tech import step_tech
+        tp = make_history().tech
+        races = [make_race("a"), make_race("b", tech={"farming": 2.0})]
+        for seed, frac in ((0, 0.02), (1, 0.3), (2, 0.9)):
+            for passes in (1, 2):
+                a, b = sparse_state(self.w.n, seed, frac), sparse_state(self.w.n, seed, frac)
+                for _ in range(3):
+                    step_tech(a, self.w, races, {**tp, "passes": passes}, 1.5)
+                    step_tech_full(b, self.w, races, {**tp, "passes": passes}, 1.5)
+                with self.subTest(seed=seed, passes=passes):
+                    self.assertTrue(np.array_equal(a.T, b.T))
+
+
+def inherit_full(st, world, races, rules, births, regions, ok):
+    """The original inherit: neighbour sums and gates over every cell."""
+    from worldhistory.lineage import _born, convert, native, progress, smoothstep
+    for rule in _born(rules):
+        s, p = rule["s"], rule["p"]
+        spec = races[s].emerge
+        Ps, Pp = st.P[s], st.P[p]
+        near_s = Ps + world.nb_sum(Ps)
+        near = near_s + Pp + world.nb_sum(Pp)
+        share = np.where(near > 0, near_s / np.maximum(near, 1e-12), 0.0)
+        where = (Pp > 0) & (births[p] > 0) & (near_s > 0) & ok[s]
+        if spec["spread"] == "region" and spec["region"]:
+            where &= regions(spec["region"])
+        if spec["mode"] == "ritual":
+            gate = np.ones(world.n)
+        else:
+            a = progress(st.O[p], races[p], races[s])
+            gate = smoothstep((a - spec["mix_min"]) / max(spec["birth_sure"] - spec["mix_min"], 1e-12))
+        amt = np.where(where, spec["dominance"] * gate * share * births[p], 0.0)
+        cells = np.flatnonzero(amt > 0)
+        if len(cells):
+            convert(st, p, s, cells, amt[cells], O_new=np.repeat(native(races[s])[:, None], len(cells), 1))
+
+
+class InheritFastTest(unittest.TestCase):
+    def test_inherit_matches_full(self):
+        from tests.helpers import globe_world, make_history, make_race
+        from tests.test_lineage import D, PARENT_TOL, SUB_TOL
+        from worldhistory.config import link
+        from worldhistory.habitat import native_optima
+        from worldhistory.lineage import inherit, init_rules
+        from worldhistory.regions import Regions
+        w = globe_world(2)
+        for seed, emerge in enumerate(({}, {"mode": "ritual"}, {"spread": "region", "region": "north"})):
+            races = link(make_history(regions={"north": {"kind": "box", "lat": [0, 90], "lon": [-180, 180]}}), [make_race("p", tolerance=PARENT_TOL),
+                                          make_race("s", tolerance=SUB_TOL, emerge={"parent": "p", **emerge})])
+            regs = Regions(w, {"north": {"kind": "box", "lat": [0, 90], "lon": [-180, 180]}})
+            out = []
+            for fn in (inherit, inherit_full):
+                rng = np.random.default_rng(seed)
+                st = new_state(2, w.n, native_optima(races), seed=0)
+                st.P[:] = np.where(rng.random((2, w.n)) < 0.1, rng.random((2, w.n)) * 500, 0.0)
+                st.O[0, D] = rng.uniform(0, 2500, w.n).astype(st.O.dtype)
+                births = np.where(rng.random((2, w.n)) < 0.7, rng.random((2, w.n)) * 20, 0.0)
+                ok = rng.random((2, w.n)) < 0.8
+                rules = init_rules(races)
+                rules[0]["origin"] = 0
+                fn(st, w, races, rules, births, regs, ok)
+                out.append(st)
+            with self.subTest(emerge=emerge):
+                self.assertTrue(np.array_equal(out[0].P, out[1].P))
+                self.assertTrue(np.array_equal(out[0].O, out[1].O))
+                self.assertGreater(out[0].P[1].sum(), 0)
+
+
+class ChangeCellsTest(unittest.TestCase):
+    def test_compiled_scan_matches_numpy(self):
+        if S._change_cells_jit is None:
+            self.skipTest("numba not installed")
+        rng = np.random.default_rng(11)
+        n = 20000
+        for seed in range(5):
+            P = np.where(rng.random(n) < 0.3, rng.random(n) * 1e3, 0.0)
+            k = rng.choice(n, 400, replace=False)
+            P[k[:100]] = -0.0
+            P[k[100:200]] = rng.random(100) * 1e-12
+            P[k[200:250]] = -rng.random(50)
+            P[k[250:270]] = np.nan
+            P[k[270:280]] = -np.nan
+            P[k[280:290]] = 1e-12
+            C = np.where(rng.random((1, n)) < 0.05, rng.random((1, n)), 0.0)
+            C[0, k[290:300]] = np.nan
+            C[0, k[300:310]] = -0.0
+            src, dst = rng.integers(0, n, 300), rng.integers(0, n, 300)
+            for wide in ([C], []):
+                with self.subTest(seed=seed, wide=len(wide)):
+                    self.assertTrue(np.array_equal(S.change_cells(P, wide, src, dst), S._change_cells_np(P, wide, src, dst)))
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/test_habitat.py b/tests/test_habitat.py
new file mode 100644
index 0000000..26f1d95
--- /dev/null
+++ b/tests/test_habitat.py
@@ -0,0 +1,104 @@
+import unittest
+
+import numpy as np
+
+from tests.helpers import line_world, make_race
+from worldhistory.config import CONDITIONS
+from worldhistory.habitat import comfort, environment, fitness, intervals, native_optima, quality, suitability
+
+DEPTH = CONDITIONS.index("depth")
+GRAV = CONDITIONS.index("gravity")
+
+
+class HabitatTest(unittest.TestCase):
+    def test_intervals_span_to_neighbour_midpoints(self):
+        w = line_world(4)
+        lo, hi = intervals(w, np.array([0.0, 100.0, 300.0, 300.0]))
+        np.testing.assert_array_equal(lo, [0, 50, 200, 300])
+        np.testing.assert_array_equal(hi, [50, 200, 300, 300])
+
+    def test_environment_depth_and_temperature(self):
+        w = line_world(3, ocean=[False, True, True], elevation_m=[10, -200, -1000], T_mean=[20, 15, 15],
+                       bottom_temp_c=[0, 8, 3])
+        env = environment(w)
+        np.testing.assert_array_equal(env.x["depth"], [0, 200, 1000])
+        np.testing.assert_array_equal(env.x["temperature"], [20, 8, 3])
+        np.testing.assert_array_equal(env.hi["depth"], [100, 600, 1000])
+
+    def test_suitability_terms_gates_realm(self):
+        w = line_world(4, ocean=[False, False, False, True], holdridge=[20, 23, 20, 0], T_mean=[20, 20, 0, 20])
+        r = make_race(habitat={"realm": "land",
+                               "terms": [{"p": "land", "w": 0.25}, {"p": "forest", "w": 0.5}],
+                               "gates": [{"require": "warm:0:20", "realm": "land"}]})
+        np.testing.assert_allclose(suitability(w, r), [0.75, 0.25, 0.0, 0.0])
+        pen = make_race(habitat={"realm": "both", "terms": [{"p": "land", "w": 1}, {"p": "sea", "w": 1}],
+                                 "gates": [{"penalty": "sea", "factor": 0.5}]})
+        np.testing.assert_allclose(suitability(w, pen), [1, 1, 1, 0.5])
+
+    def test_quality_normalised_to_p99(self):
+        s = np.r_[np.linspace(0, 1, 101), 50.0]         # one outlier must not squash everyone else
+        q = quality(s, make_race(habitat={"qmax": 0.7}), np.ones_like(s))
+        self.assertAlmostEqual(float(q[100]), 0.7, places=6)
+        self.assertLessEqual(float(q.max()), 0.7)
+
+    def test_quality_all_zero(self):
+        q = quality(np.zeros(5), make_race(), np.ones(5))
+        np.testing.assert_array_equal(q, 0)
+
+    def test_fitness_uses_distance_to_interval_and_reach(self):
+        w = line_world(3, ocean=[False, True, True], elevation_m=[0, -100, -400])
+        env = environment(w)
+        r = make_race(tolerance={"depth": {"optimum": 0, "width": 50, "lo": 0, "hi": 5000}})
+        O = native_optima([r])[0]
+        O = np.repeat(O[:, None], 3, 1)
+        f = fitness(env, O, r)
+        self.assertAlmostEqual(float(f[0]), 1.0)                                  # interval [0, 50] holds 0
+        self.assertAlmostEqual(float(f[1]), np.exp(-0.5), places=6)               # interval [50, 250]: d = 50 = 1 width
+        self.assertAlmostEqual(float(fitness(env, O, r, reach=1.0)[1]), np.exp(-0.125), places=6)
+
+    def test_unlisted_condition_has_no_effect(self):
+        w = line_world(2, gravity_g=[1.0, 0.35])
+        r = make_race()
+        O = np.repeat(native_optima([r])[0][:, None], 2, 1)
+        np.testing.assert_array_equal(fitness(environment(w), O, r), [1, 1])
+
+    def test_comfort_piecewise(self):
+        r = make_race(tolerance={"gravity": {"optimum": 1.0, "width": 0.2, "lo": 0.3, "hi": 1.3,
+                                             "comfort": [[1.0, 1.0], [0.35, 1.2]]}})
+        O = np.zeros((len(CONDITIONS), 3))
+        O[GRAV] = [1.0, 0.675, 0.35]
+        np.testing.assert_allclose(comfort(O, r), [1.0, 1.1, 1.2])
+
+    def test_lopsided_widths(self):
+        w = line_world(5, gravity_g=[1.05, 1.05, 0.75, 0.45, 0.45])
+        r = make_race(tolerance={"gravity": {"optimum": 0.75, "width_lo": 0.3, "width_hi": 0.15, "lo": 0.3, "hi": 1.0}})
+        O = np.repeat(native_optima([r])[0][:, None], 5, 1)
+        f = fitness(environment(w), O, r)
+        # cell 0 spans [1.05, 1.05]: 0.3 above = 2 upper widths; cell 4 spans [0.45, 0.45]: 0.3 below = 1 lower width
+        self.assertAlmostEqual(float(f[0]), np.exp(-2.0), places=6)
+        self.assertAlmostEqual(float(f[4]), np.exp(-0.5), places=6)
+
+    def test_width_sets_both_sides(self):
+        r = make_race(tolerance={"depth": {"optimum": 0, "width": 50, "lo": 0, "hi": 5000}})
+        t = r.tolerance["depth"]
+        self.assertEqual((t.width_lo, t.width_hi), (50, 50))
+
+    def test_rain_is_log_and_ignored_at_sea(self):
+        w = line_world(3, ocean=[False, False, True], P_ann=[100.0, 1000.0, 1000.0], elevation_m=[10, 10, -100])
+        env = environment(w)
+        np.testing.assert_allclose(env.x["rain"][:2], [2.0, 3.0])
+        self.assertEqual(env.lo["rain"][2], -np.inf)
+        self.assertEqual(env.hi["rain"][2], np.inf)
+        r = make_race(habitat={"realm": "both", "terms": [{"p": "land", "w": 1}]},
+                      tolerance={"rain": {"optimum": 2.0, "width": 0.3, "lo": 1.5, "hi": 3.5}})
+        O = np.repeat(native_optima([r])[0][:, None], 3, 1)
+        f = fitness(env, O, r)
+        self.assertEqual(float(f[2]), 1.0)                     # no rain effect in the sea
+        self.assertLess(float(f[1]), float(f[0]))
+
+    def test_default_strain_off_native(self):
+        # adapted away from the native optimum: livable, not good — the peak falls to 0.35 at the range edge
+        r = make_race(tolerance={"gravity": {"optimum": 1.0, "width": 0.2, "lo": 0.3, "hi": 1.3}})
+        O = np.zeros((len(CONDITIONS), 4))
+        O[GRAV] = [1.0, 0.3, 0.65, 1.3]
+        np.testing.assert_allclose(comfort(O, r), [1.0, 0.35, 0.35 ** 0.25, 0.35], rtol=1e-9)
diff --git a/tests/test_hazards.py b/tests/test_hazards.py
new file mode 100644
index 0000000..f3f366c
--- /dev/null
+++ b/tests/test_hazards.py
@@ -0,0 +1,106 @@
+import unittest
+
+import numpy as np
+
+from tests.helpers import globe_world, line_world, make_history, make_race
+from worldhistory.hazards import make_hazard
+from worldhistory.regions import Regions
+from worldhistory.state import new_state
+
+
+def hist(**hz):
+    return make_history(regions={"north": {"kind": "box", "lat": [30, 90], "lon": [-180, 180]}}, hazard=[hz])
+
+
+class HazardTest(unittest.TestCase):
+    def setUp(self):
+        self.w = globe_world(1)
+        self.st = new_state(1, self.w.n, np.zeros((1, 6)), seed=0)
+        self.st.P[:] = 100.0
+
+    def test_static_in_window_and_region(self):
+        h = hist(kind="static", region="north", mortality=0.1, years=[10, 30])
+        hz = make_hazard(h.hazards[0])
+        reg = Regions(self.w, h.regions)
+        north = reg("north")
+        hz.step(self.st, self.w, [make_race()], reg, 0, np.random.default_rng(0))
+        self.assertTrue((self.st.P == 100).all())
+        hz.step(self.st, self.w, [make_race()], reg, 10, np.random.default_rng(0))
+        np.testing.assert_allclose(self.st.P[0, north], 90.0)
+        np.testing.assert_allclose(self.st.P[0, ~north], 100.0)
+
+    def test_static_push_drives_people_off_land(self):
+        # megafauna: on land in the region, a share flees each step (not killed); drift is told to avoid those cells
+        w = line_world(4, ocean=[False, False, True, True])
+        h = make_history(regions={"coast": {"kind": "box", "lat": [-1, 1], "lon": [-1, 99]}},
+                         hazard=[{"kind": "static", "name": "beasts", "region": "coast", "mortality": 0.0,
+                                  "years": [0, 100], "push": 0.2}])
+        st = new_state(1, 4, np.zeros((1, 6)), seed=0)
+        st.P[:] = 100.0
+        hz = make_hazard(h.hazards[0])
+        hz.step(st, w, [make_race()], Regions(w, h.regions), 0, np.random.default_rng(0))
+        np.testing.assert_allclose(st.P[0], 100.0)
+        np.testing.assert_allclose(st.push[0], [20, 20, 0, 0])
+        np.testing.assert_allclose(hz.danger, [1, 1, 0, 0])
+        hz.step(st, w, [make_race()], Regions(w, h.regions), 200, np.random.default_rng(0))
+        np.testing.assert_allclose(hz.danger, 0)                          # outside its years: no danger
+
+    def test_roaming_local_stay_home_and_raid(self):
+        h = hist(kind="roaming", start=20, home_region="north", home_near="a", local=10, wanderers=0,
+                 radius_km=600, raid_chance=0.5, growth=0.0, infight=0.0, raid_death=0.0)
+        hz = make_hazard(h.hazards[0])
+        reg = Regions(self.w, h.regions)
+        rng = np.random.default_rng(3)
+        hz.step(self.st, self.w, [make_race()], reg, 10, rng)
+        self.assertTrue((self.st.P == 100).all())                    # not started yet
+        stats = [hz.step(self.st, self.w, [make_race()], reg, t, rng) for t in range(20, 200, 10)]
+        self.assertTrue(reg("north")[hz.home].all())
+        self.assertEqual(len(hz.home), 10)
+        self.assertGreater(sum(s["raids"] for s in stats), 0)
+        self.assertTrue((self.st.P[0, reg("north")] < 100).any())
+        np.testing.assert_array_equal(self.st.P[0, self.w.lat < 20], 100.0)   # north (lat ≥ 30) + 600 km only
+
+    def test_roaming_wanderer_leaves_home_region(self):
+        h = hist(kind="roaming", start=0, home_region="north", local=0, wanderers=1, wander_stop=1.0,
+                 wander_step_km=1500.0, growth=0.0, infight=0.0)
+        hz = make_hazard(h.hazards[0])
+        reg = Regions(self.w, h.regions)
+        rng = np.random.default_rng(0)
+        for t in range(0, 400, 10):
+            hz.step(self.st, self.w, [make_race()], reg, t, rng)
+        self.assertTrue((self.st.P[0, ~reg("north")] < 100).any())
+
+    def test_roaming_growth_and_deaths(self):
+        h = hist(kind="roaming", start=0, home_region="north", local=100, wanderers=0, growth=0.1,
+                 clutch=[2, 4], infight=0.05, raid_chance=0.0)
+        hz = make_hazard(h.hazards[0])
+        reg = Regions(self.w, h.regions)
+        rng = np.random.default_rng(0)
+        births = deaths = 0
+        for t in range(0, 500, 10):
+            s = hz.step(self.st, self.w, [make_race()], reg, t, rng)
+            births += s["births"]
+            deaths += s["deaths"]
+        self.assertGreater(births, 0)
+        self.assertGreater(deaths, 0)
+        self.assertEqual(len(hz.home), 100 + births - deaths)
+
+    def test_roaming_growth_levels_off_at_cap(self):
+        # growth 0.1/step uncapped would give 20·1.1^100 ≈ 275k units; logistic cap holds it near 60
+        h = hist(kind="roaming", start=0, home_region="north", local=20, wanderers=0, growth=0.1,
+                 clutch=[2, 4], infight=0.0, raid_chance=0.0, cap=60)
+        hz = make_hazard(h.hazards[0])
+        reg = Regions(self.w, h.regions)
+        rng = np.random.default_rng(0)
+        for t in range(0, 1000, 10):
+            hz.step(self.st, self.w, [make_race()], reg, t, rng)
+        self.assertGreater(len(hz.home), 40)
+        self.assertLessEqual(len(hz.home), 60 + 4)
+
+    def test_roaming_fallback_home(self):
+        self.st.P[:] = 0.0
+        h = hist(kind="roaming", start=0, home_region="north", home_near="a", local=5, wanderers=0)
+        hz = make_hazard(h.hazards[0])
+        reg = Regions(self.w, h.regions)
+        hz.step(self.st, self.w, [make_race()], reg, 0, np.random.default_rng(0))
+        self.assertTrue(reg("north")[hz.home].all())
diff --git a/tests/test_lineage.py b/tests/test_lineage.py
new file mode 100644
index 0000000..403bccd
--- /dev/null
+++ b/tests/test_lineage.py
@@ -0,0 +1,516 @@
+import unittest
+
+import numpy as np
+
+from tests.helpers import line_world, make_race
+from worldhistory.config import CONDITIONS, link
+from worldhistory.habitat import native_optima
+from worldhistory.lineage import birth_step, chance, init_rules, native, progress, return_step
+from worldhistory.regions import Regions
+from worldhistory.state import new_state
+from tests.helpers import make_history
+
+D, G = CONDITIONS.index("depth"), CONDITIONS.index("gravity")
+PARENT_TOL = {"depth": {"optimum": 0, "width_lo": 60, "width_hi": 100, "lo": 0, "hi": 6000, "rate": 2.0},
+              "gravity": {"optimum": 1.0, "width": 0.5, "lo": 0.3, "hi": 1.3, "rate": 0.01}}
+SUB_TOL = {"depth": {"optimum": 2000, "width_lo": 400, "width_hi": 1500, "lo": 500, "hi": 6000, "rate": 2.0},
+           "gravity": {"optimum": 0.5, "width": 0.5, "lo": 0.3, "hi": 1.3, "rate": 0.01}}
+
+
+def pair(n=6, emerge=None, **fields):
+    w = line_world(n, **fields)
+    races = link(make_history(), [make_race("p", tolerance=PARENT_TOL),
+                                  make_race("s", tolerance=SUB_TOL, emerge={"parent": "p", **(emerge or {})})])
+    st = new_state(2, n, native_optima(races), seed=0)
+    return w, races, st
+
+
+class ProgressTest(unittest.TestCase):
+    def test_weighted_mean_over_distinct_axes(self):
+        w, (p, s), st = pair(1)
+        O = native(p)[:, None].copy()
+        O[D, 0] = 1000.0                      # halfway on depth (weight 2000/100 = 20)
+        a = progress(O, p, s)                 # gravity untouched (weight 0.5/0.5 = 1)
+        self.assertAlmostEqual(float(a[0]), 20 * 0.5 / 21, places=6)
+
+    def test_overshoot_clips_and_equal_axes_ignored(self):
+        w, (p, s), st = pair(1)
+        O = native(p)[:, None].copy()
+        O[D, 0], O[G, 0] = 3000.0, 0.2       # both past the sub-species' optimum
+        self.assertAlmostEqual(float(progress(O, p, s)[0]), 1.0)
+        self.assertEqual(float(progress(O, p, p)[0]), 0.0)                # no distinct axis → 0
+
+    def test_back_toward_parent(self):
+        w, (p, s), st = pair(1)
+        O = native(s)[:, None].copy()
+        O[D, 0] = 1000.0                      # Havs risen halfway; toward the parent the width is s's lower side (400)
+        a = progress(O, s, p)                 # weights: depth 2000/400 = 5, gravity 1
+        self.assertAlmostEqual(float(a[0]), 5 * 0.5 / 6, places=6)
+
+
+class ChanceTest(unittest.TestCase):
+    def test_floor_ceiling_and_monotone(self):
+        a = np.array([0.0, 0.39, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0])
+        c = chance(a, 0.4, 0.8, 0.2)
+        np.testing.assert_allclose(c[[0, 1, 2]], 0.0)
+        np.testing.assert_allclose(c[[6, 7]], 0.2)
+        self.assertTrue(np.all(np.diff(c) >= 0))
+        self.assertAlmostEqual(float(c[4]), 0.1)                         # smoothstep(0.5) = 0.5
+
+    def test_equal_bounds_is_a_step(self):
+        np.testing.assert_allclose(chance(np.array([0.0, 0.5]), 0.5, 0.5, 1.0), [0.0, 1.0])
+
+
+class BirthTest(unittest.TestCase):
+    def deep(self, st, cells, depth=2000.0):
+        st.O[0, D, cells] = depth
+
+    def test_no_birth_below_min(self):
+        w, races, st = pair()
+        st.P[0] = 100
+        self.deep(st, slice(None), 500.0)                                # a = 20·0.25/21 < 0.4
+        rules = init_rules(races)
+        ok = np.ones((2, 6), bool)
+        for s in range(200):
+            self.assertEqual(birth_step(st, w, races, rules, Regions(w, {}), np.random.default_rng(s), ok), [])
+
+    def test_birth_resets_to_native_and_is_once(self):
+        w, races, st = pair(emerge={"birth_rate": 1.0})
+        st.P[0] = [100, 100, 0, 0, 100, 100]
+        self.deep(st, slice(None))
+        st.T[0, 0, :] = 0.5
+        rules = init_rules(races)
+        ok = np.ones((2, 6), bool)
+        new = birth_step(st, w, races, rules, Regions(w, {}), np.random.default_rng(0), ok)
+        self.assertEqual(len(new), 1)
+        c = new[0]["cell"]
+        self.assertEqual(float(st.P[0, c]), 0.0)                         # converted wholly
+        self.assertAlmostEqual(float(st.P[1, c]), 100.0)
+        np.testing.assert_allclose(st.O[1, :, c], native(races[1]))     # baseline shift = 1
+        self.assertAlmostEqual(float(st.T[1, 0, c]), 0.5)                # tech kept
+        self.assertGreater(new[0]["a"], 0.8)
+        for s in range(20):
+            self.assertEqual(birth_step(st, w, races, rules, Regions(w, {}), np.random.default_rng(s + 1), ok), [])
+
+    def test_birth_needs_native_fit(self):
+        # Review Focus 1: where the sub-species' own curves do not fit, no birth
+        w, races, st = pair(emerge={"birth_rate": 1.0})
+        st.P[0] = 100
+        self.deep(st, slice(None))
+        ok = np.ones((2, 6), bool)
+        ok[1] = False
+        rules = init_rules(races)
+        for s in range(50):
+            self.assertEqual(birth_step(st, w, races, rules, Regions(w, {}), np.random.default_rng(s), ok), [])
+
+    def test_origin_rules_still_apply(self):
+        w, races, st = pair(emerge={"birth_rate": 1.0, "settlement_density": 0.001})
+        st.P[0] = [10, 20, 90, 30, 10, 10]
+        self.deep(st, slice(None))
+        rules = init_rules(races)
+        birth_step(st, w, races, rules, Regions(w, {}), np.random.default_rng(0), np.ones((2, 6), bool))
+        self.assertEqual(rules[0]["origin"], 2)                          # the densest settlement
+
+    def test_small_groups_rarely_found(self):
+        counts = []
+        for P in (1.0, 1000.0):
+            n = 0
+            for s in range(300):
+                w = line_world(1)
+                races = link(make_history(), [make_race("p", tolerance=PARENT_TOL, demography={"founder": 50.0}),
+                                              make_race("s", tolerance=SUB_TOL, emerge={"parent": "p"})])
+                st = new_state(2, 1, native_optima(races), seed=0)
+                st.P[0] = P
+                st.O[0, D] = 2000.0
+                if birth_step(st, w, races, init_rules(races), Regions(w, {}), np.random.default_rng(s),
+                              np.ones((2, 1), bool)):
+                    n += 1
+            counts.append(n)
+        self.assertLess(counts[0], counts[1] / 5)                        # 1/51 vs 0.95 of 0.2 per step
+
+
+class ReturnTest(unittest.TestCase):
+    def test_return_births_parents_many_times_capped(self):
+        # Review Focus 3: returns may fire in several cells, never convert more than the group, never log emergence
+        w, races, st = pair(emerge={"birth_rate": 1.0})
+        rules = init_rules(races)
+        rules[0]["origin"], rules[0]["year"] = 0, 0                       # already born
+        st.P[1] = [50, 50, 50, 0, 0, 0]
+        st.O[1, :, :] = native(races[1])[:, None]
+        st.O[1, D, :3] = 0.0                                             # risen all the way back
+        ok = np.ones((2, 6), bool)
+        rets = return_step(st, w, races, rules, np.random.default_rng(0), ok)
+        self.assertEqual(sorted(r["cell"] for r in rets), [0, 1, 2])
+        np.testing.assert_allclose(st.P[1, :3], 0.0)
+        np.testing.assert_allclose(st.P[0, :3], 50.0)
+        np.testing.assert_allclose(st.O[0, :, 0], native(races[0]))
+        self.assertTrue((st.P >= 0).all())
+
+    def test_no_return_before_birth_or_for_ritual(self):
+        w, races, st = pair(emerge={"birth_rate": 1.0})
+        rules = init_rules(races)
+        st.P[1] = 50
+        st.O[1, D] = 0.0
+        self.assertEqual(return_step(st, w, races, rules, np.random.default_rng(0), np.ones((2, 6), bool)), [])
+
+
+class RitualTest(unittest.TestCase):
+    def test_ritual_birth_ignores_progress(self):
+        w, races, st = pair(emerge={"mode": "ritual", "trigger": 1.0})
+        st.P[0] = 100                                                    # optima untouched: a = 0
+        rules = init_rules(races)
+        new = birth_step(st, w, races, rules, Regions(w, {}), np.random.default_rng(0), np.ones((2, 6), bool))
+        self.assertEqual(len(new), 1)
+
+
+from worldhistory.lineage import displacement, inherit, rituals
+
+
+class DisplacementTest(unittest.TestCase):
+    def test_less_fit_race_squeezed(self):
+        w, races, st = pair(3)
+        rules = init_rules(races)
+        rules[0]["origin"] = 0
+        K = np.full((2, 3), 100.0)
+        P = np.array([[30.0, 50.0, 10.0], [10.0, 0.0, 30.0]])
+        fit = np.array([[0.2, 0.9, 0.9], [0.9, 0.1, 0.1]])
+        displacement(K, P, fit, rules, races)
+        np.testing.assert_allclose(K[0], [100 * (1 - 0.9 * 10 / 40), 100.0, 100.0])      # parent less fit in cell 0
+        np.testing.assert_allclose(K[1], [100.0, 100.0, 100 * (1 - 0.9 * 10 / 40)])      # sub less fit in cell 2
+        self.assertTrue((K >= 0).all())
+
+    def test_no_effect_before_birth_or_on_strangers(self):
+        # Review Focus 5
+        w, races, st = pair(3)
+        rules = init_rules(races)
+        K = np.full((2, 3), 100.0)
+        displacement(K, np.ones((2, 3)), np.array([[0.1] * 3, [0.9] * 3]), rules, races)
+        np.testing.assert_allclose(K, 100.0)
+
+
+class InheritTest(unittest.TestCase):
+    def test_adapted_neighbours_bear_the_dominant_race(self):
+        w, races, st = pair(4)
+        rules = init_rules(races)
+        rules[0]["origin"] = 0
+        st.P[1] = [100, 0, 0, 0]
+        st.P[0] = [0, 100, 100, 100]
+        st.O[0, D, 1] = 2000.0                    # adapted (a = 20/21 ≥ sure)
+        st.O[0, D, 2] = 2000.0                    # adapted but not in contact
+        births = np.zeros((2, 4))
+        births[0] = 10.0
+        inherit(st, w, races, rules, births, Regions(w, {}), np.ones((2, w.n), bool))
+        share = 100 / 300                         # around cell 1 (cells 0..2): 100 sub-species, 200 parents
+        self.assertAlmostEqual(float(st.P[1, 1]), 10 * share, places=6)
+        np.testing.assert_allclose(st.O[1, :, 1], native(races[1]))
+        self.assertEqual(float(st.P[1, 3]), 0.0)
+
+    def test_no_dominant_births_where_native_curves_do_not_fit(self):
+        w, races, st = pair(2)
+        rules = init_rules(races)
+        rules[0]["origin"] = 0
+        st.P[1] = [100, 0]
+        st.P[0] = [0, 100]
+        st.O[0, D, 1] = 2000.0
+        births = np.zeros((2, 2))
+        births[0] = 10.0
+        ok = np.ones((2, 2), bool)
+        ok[1, 1] = False
+        inherit(st, w, races, rules, births, Regions(w, {}), ok)
+        self.assertEqual(float(st.P[1, 1]), 0.0)
+
+    def test_holdouts_refuse_to_mix(self):
+        w, races, st = pair(2)
+        rules = init_rules(races)
+        rules[0]["origin"] = 0
+        st.P[1] = [100, 0]
+        st.P[0] = [0, 100]                        # optima native: a = 0 < mix_min
+        births = np.zeros((2, 2))
+        births[0] = 10.0
+        inherit(st, w, races, rules, births, Regions(w, {}), np.ones((2, w.n), bool))
+        self.assertEqual(float(st.P[1, 1]), 0.0)
+
+    def test_ritual_race_is_plainly_dominant(self):
+        w, races, st = pair(2, emerge={"mode": "ritual"})
+        rules = init_rules(races)
+        rules[0]["origin"] = 0
+        st.P[1] = [100, 0]
+        st.P[0] = [0, 100]                        # a = 0 but no gate in ritual mode
+        births = np.zeros((2, 2))
+        births[0] = 10.0
+        inherit(st, w, races, rules, births, Regions(w, {}), np.ones((2, w.n), bool))
+        self.assertAlmostEqual(float(st.P[1, 1]), 10 * 100 / 200, places=6)
+
+
+class RitualConvertTest(unittest.TestCase):
+    def test_rituals_cost_hundreds_and_scale_sublinearly(self):
+        w, races, st = pair(2, emerge={"mode": "ritual", "ritual_rate": 0.5, "ritual_cost": 300,
+                                       "ritual_converts": 30})
+        rules = init_rules(races)
+        rules[0]["origin"] = 0
+        st.P[1] = [400, 0]
+        st.P[0] = [0, 100000]
+        out = rituals(st, w, races, rules, np.random.default_rng(0))
+        n = sum(r["converts"] for r in out) / 30
+        self.assertGreater(n, 0)
+        self.assertAlmostEqual(float(st.P[0, 1]), 100000 - 330 * n)
+        self.assertAlmostEqual(float(st.P[1].sum()), 400 + 30 * n)
+
+    def test_rituals_never_overdraw(self):
+        # Review Focus 4
+        w, races, st = pair(2, emerge={"mode": "ritual", "ritual_rate": 50.0})
+        rules = init_rules(races)
+        rules[0]["origin"] = 0
+        st.P[1] = [10000, 0]
+        st.P[0] = [0, 500]
+        rituals(st, w, races, rules, np.random.default_rng(0))
+        self.assertTrue((st.P >= 0).all())
+
+
+class SummaryTest(unittest.TestCase):
+    def test_events_aggregate_per_step(self):
+        from worldhistory.lineage import summarise
+        ev = [{"race": "p", "from": "s", "cell": 1, "year": 10, "people": 5.0},
+              {"race": "p", "from": "s", "cell": 2, "year": 10, "people": 7.0}]
+        self.assertEqual(summarise(ev, ("people",)), [{"race": "p", "from": "s", "year": 10, "groups": 2, "people": 12.0}])
+        self.assertEqual(summarise([], ("people",)), [])
+
+    def test_split_by_source_and_dicts_summed(self):
+        from worldhistory.lineage import summarise
+        ev = [{"race": "p", "from": "s", "cell": 1, "year": 10, "people": 5.0},
+              {"race": "p", "from": "t", "cell": 2, "year": 10, "people": 7.0},
+              {"race": "p", "from": "s", "cell": 3, "year": 10, "people": 1.0}]
+        self.assertEqual(summarise(ev, ("people",)),
+                         [{"race": "p", "from": "s", "year": 10, "groups": 2, "people": 6.0},
+                          {"race": "p", "from": "t", "year": 10, "groups": 1, "people": 7.0}])
+        ev = [{"race": "b", "cell": 1, "year": 10, "victims": 4.0, "dead": {"p": 1.0, "o": 3.0}},
+              {"race": "b", "cell": 2, "year": 10, "victims": 2.0, "dead": {"p": 2.0}}]
+        self.assertEqual(summarise(ev, ("victims", "dead")),
+                         [{"race": "b", "year": 10, "groups": 2, "victims": 6.0, "dead": {"p": 3.0, "o": 3.0}}])
+        ev.append({"race": "b", "cell": -1, "year": 10, "victims": 0.0, "backlash": 9.0})    # no "dead" field
+        self.assertEqual(summarise(ev, ("victims", "dead", "backlash"))[0]["dead"], {"p": 3.0, "o": 3.0})
+
+    def test_backlash_record_is_summed_not_counted(self):
+        from worldhistory.lineage import summarise
+        ev = [{"race": "b", "cell": 3, "year": 10, "victims": 300.0, "converts": 30.0},
+              {"race": "b", "cell": -1, "year": 10, "victims": 0.0, "converts": 0.0, "backlash": 70.0}]
+        self.assertEqual(summarise(ev, ("victims", "converts", "backlash")),
+                         [{"race": "b", "year": 10, "groups": 1, "victims": 300.0, "converts": 30.0, "backlash": 70.0}])
+
+def trio(n=3, emerge=None):
+    """Parent p, sub-species s, and an unrelated race o."""
+    w = line_world(n)
+    races = link(make_history(), [make_race("p", tolerance=PARENT_TOL), make_race("o"),
+                                  make_race("s", tolerance=SUB_TOL, emerge={"parent": "p", **(emerge or {})})])
+    st = new_state(3, n, native_optima(races), seed=0)
+    return w, races, st
+
+
+class AfterTest(unittest.TestCase):
+    def test_no_birth_before_the_earliest_year(self):
+        for mode in ("ritual", "adapt"):
+            em = {"mode": mode, "trigger": 1.0, "after": 100} if mode == "ritual" else \
+                 {"after": 100, "birth_min": 0.0, "birth_sure": 0.0, "birth_rate": 1.0}
+            w, races, st = pair(emerge=em)
+            st.P[0] = 1e6
+            rules = init_rules(races)
+            ok = np.ones((2, 6), bool)
+            st.t = 90
+            self.assertEqual(birth_step(st, w, races, rules, Regions(w, {}), np.random.default_rng(0), ok), [], mode)
+            st.t = 100
+            self.assertEqual(len(birth_step(st, w, races, rules, Regions(w, {}), np.random.default_rng(0), ok)), 1,
+                             mode)
+
+
+class CrisisTest(unittest.TestCase):
+    def test_ritual_birth_only_where_parents_were_almost_eradicated(self):
+        # cell 1 falls from 10,000 to 500 (−95 %), cell 3 from 10,000 to 5,000 (−50 %), cell 4 tiny all along
+        w, races, st = pair(emerge={"mode": "ritual", "trigger": 1.0, "crisis_drop": 0.9, "crisis_min": 1000,
+                                    "crisis_years": 200})
+        rules = init_rules(races)
+        ok = np.ones((2, 6), bool)
+        rng = np.random.default_rng(0)
+        st.P[0] = [0, 10000, 0, 10000, 50, 0]
+        st.t = 0
+        self.assertEqual(birth_step(st, w, races, rules, Regions(w, {}), rng, ok), [])   # peaks recorded, no crisis
+        st.P[0] = [0, 500, 0, 5000, 5, 0]
+        st.t = 10
+        new = birth_step(st, w, races, rules, Regions(w, {}), rng, ok)
+        self.assertEqual([e["cell"] for e in new], [1])
+
+    def _off_habitat(self, km):
+        # crisis in cell 1 (10,000 → 500), where s cannot live; s can live only in cell 3, 200 km away
+        w, races, st = pair(emerge={"mode": "ritual", "trigger": 1.0, "crisis_drop": 0.9, "crisis_min": 1000,
+                                    "convert": 0.5, **({"crisis_km": km} if km is not None else {})})
+        rules = init_rules(races)
+        ok = np.ones((2, 6), bool)
+        ok[1] = [False, False, False, True, False, False]
+        rng = np.random.default_rng(0)
+        st.P[0] = [0, 10000, 0, 0, 0, 0]
+        st.t = 0
+        birth_step(st, w, races, rules, Regions(w, {}), rng, ok)
+        st.P[0] = [0, 500, 0, 0, 0, 0]
+        st.t = 10
+        return st, birth_step(st, w, races, rules, Regions(w, {}), rng, ok)
+
+    def test_crisis_off_habitat_births_in_nearest_habitat_within_range(self):
+        st, new = self._off_habitat(250)
+        self.assertEqual([e["cell"] for e in new], [3])
+        self.assertEqual(new[0]["crisis"], 1)
+        np.testing.assert_allclose(st.P[1], [0, 0, 0, 250, 0, 0])      # half the 500 survivors, moved to cell 3
+        np.testing.assert_allclose(st.P[0], [0, 250, 0, 0, 0, 0])
+
+    def test_crisis_off_habitat_no_birth_out_of_range_or_by_default(self):
+        for km in (150, None):
+            st, new = self._off_habitat(km)
+            self.assertEqual(new, [], km)
+            np.testing.assert_allclose(st.P[1], 0)
+
+    def test_no_crisis_no_birth(self):
+        w, races, st = pair(emerge={"mode": "ritual", "trigger": 1.0, "crisis_drop": 0.9})
+        st.P[0] = 10000
+        rules = init_rules(races)
+        self.assertEqual(birth_step(st, w, races, rules, Regions(w, {}), np.random.default_rng(0),
+                                    np.ones((2, 6), bool)), [])
+
+
+class VictimsTest(unittest.TestCase):
+    def test_victims_from_every_nearby_race_converts_from_parents(self):
+        # Blod (s) in cell 0; parents 30,000 in cell 1; other race 90,000 in cell 1 → victims split 1:3
+        w, races, st = trio(emerge={"mode": "ritual", "ritual_rate": 0.5, "ritual_cost": 300,
+                                    "ritual_converts": 30, "victims": "all"})
+        rules = init_rules(races)
+        rules[0]["origin"] = 0
+        p, o, s = 0, 1, 2
+        st.P[s] = [400, 0, 0]
+        st.P[p] = [0, 30000, 0]
+        st.P[o] = [0, 90000, 0]
+        out = rituals(st, w, races, rules, np.random.default_rng(0))
+        k = sum(r["converts"] for r in out) / 30
+        self.assertGreater(k, 0)
+        self.assertAlmostEqual(float(st.P[p, 1]), 30000 - 30 * k - 300 * k * 0.25)
+        self.assertAlmostEqual(float(st.P[o, 1]), 90000 - 300 * k * 0.75)
+        self.assertAlmostEqual(float(st.P[s].sum()), 400 + 30 * k)
+        self.assertTrue((st.P >= 0).all())
+        dead = {}
+        for r in out:
+            for key, x in r["dead"].items():
+                dead[key] = dead.get(key, 0.0) + x
+        self.assertAlmostEqual(dead["p"], 300 * k * 0.25)
+        self.assertAlmostEqual(dead["o"], 300 * k * 0.75)
+
+
+class BacklashTest(unittest.TestCase):
+    def test_smackdown_once_then_calmer(self):
+        w, races, st = pair(2, emerge={"mode": "ritual", "ritual_rate": 0.5, "ritual_cost": 300,
+                                       "ritual_converts": 30, "backlash_victims": 600, "backlash_loss": 0.7,
+                                       "backlash_calm": 0.1})
+        rules = init_rules(races)
+        rules[0]["origin"] = 0
+        st.P[1] = [400, 0]
+        st.P[0] = [0, 100000]
+        rng = np.random.default_rng(0)
+        rituals(st, w, races, rules, rng)
+        v = rules[0]["victims"]
+        self.assertGreaterEqual(v, 600)                                  # ≥ 2 rituals: the smackdown fires
+        self.assertTrue(rules[0]["backlash"])
+        blod = 400 + v / 300 * 30
+        self.assertAlmostEqual(float(st.P[1].sum()), blod * 0.3)
+        self.assertAlmostEqual(rules[0]["rate"], 0.05)
+        before = st.P[1].sum()
+        rules[0]["victims"] += 1e6
+        rituals(st, w, races, rules, rng)
+        self.assertGreaterEqual(float(st.P[1].sum()), before)            # no second smackdown
+
+
+class ReturnRateTest(unittest.TestCase):
+    def test_return_rate_sets_the_chance(self):
+        w, races, st = pair(emerge={"birth_rate": 1.0, "return_rate": 0.0})
+        rules = init_rules(races)
+        rules[0]["origin"] = 0
+        st.P[1] = 1e6
+        st.O[1, :, :] = native(races[1])[:, None]
+        st.O[1, D, :] = 0.0
+        self.assertEqual(return_step(st, w, races, rules, np.random.default_rng(0), np.ones((2, 6), bool)), [])
+
+
+from worldhistory.lineage import raids
+
+RAID = {"mode": "ritual", "raid_rate": 5.0, "raid_size": 100, "raid_min": 1000, "raid_km": 250}
+
+
+class RaidTest(unittest.TestCase):
+    def setUp(self):
+        self.w, self.races, self.st = pair(6, emerge=RAID)
+        self.rules = init_rules(self.races)
+        self.rules[0]["origin"] = 0
+        self.st.P[1] = [5000, 0, 0, 0, 0, 0]
+        self.st.P[0] = [0, 0, 800, 0, 900, 0]         # cell 1 empty, cell 2 in reach (200 km), cell 4 out (400 km)
+        self.ok = np.ones((2, 6), bool)
+
+    def test_cells_raid_occupied_non_blod_land_in_reach(self):
+        out = raids(self.st, self.w, self.races, self.rules, np.random.default_rng(0), self.ok)
+        n = sum(r["raids"] for r in out)
+        self.assertGreater(n, 0)
+        np.testing.assert_allclose(self.st.P[1], [5000 - 100 * n, 0, 100 * n, 0, 0, 0])
+        np.testing.assert_allclose(self.st.P[0], [0, 0, 800, 0, 900, 0])     # the killing is the rituals' job
+
+    def test_no_raids_after_the_smackdown(self):
+        self.rules[0]["backlash"] = True
+        self.assertEqual(raids(self.st, self.w, self.races, self.rules, np.random.default_rng(0), self.ok), [])
+        self.assertEqual(float(self.st.P[1, 0]), 5000.0)
+
+    def test_small_groups_do_not_raid(self):
+        self.st.P[1, 0] = 900
+        self.assertEqual(raids(self.st, self.w, self.races, self.rules, np.random.default_rng(0), self.ok), [])
+
+    def test_no_raids_where_blod_cannot_live(self):
+        self.ok[1, 2] = False
+        self.assertEqual(raids(self.st, self.w, self.races, self.rules, np.random.default_rng(0), self.ok), [])
+
+
+class FrenzyTest(unittest.TestCase):
+    def test_linear_rituals_scale_with_cultists(self):
+        # ritual_power 1: rituals ≈ rate·P per step; 10,000 Blod at 1/15 → ≈ 667 rituals, 200k dead
+        w, races, st = trio(emerge={"mode": "ritual", "ritual_rate": 1 / 15, "ritual_power": 1.0,
+                                    "ritual_cost": 300, "ritual_converts": 30, "victims": "all"})
+        rules = init_rules(races)
+        rules[0]["origin"] = 0
+        st.P[2] = [10000, 0, 0]
+        st.P[0] = [0, 1e7, 0]
+        st.P[1] = [0, 1e7, 0]
+        out = rituals(st, w, races, rules, np.random.default_rng(0))
+        dead = sum(r["victims"] for r in out)
+        self.assertAlmostEqual(dead / 200000, 1.0, delta=0.08)                 # Poisson(667): sd ≈ 4 %
+
+
+class SmackdownTest(unittest.TestCase):
+    def rig(self, **em):
+        w, races, st = trio(5, emerge={"mode": "ritual", "ritual_rate": 0.0, "backlash_loss": 0.7, **em})
+        rules = init_rules(races)
+        rules[0]["origin"], rules[0]["year"] = 0, 100
+        return w, races, st, rules
+
+    def test_time_limit_fires_without_the_body_count(self):
+        w, races, st, rules = self.rig(backlash_victims=3e6, backlash_years=50)
+        st.P[2] = [1000, 0, 0, 0, 0]
+        st.t = 140
+        rituals(st, w, races, rules, np.random.default_rng(0))
+        self.assertFalse(rules[0]["backlash"])
+        st.t = 150
+        rituals(st, w, races, rules, np.random.default_rng(0))
+        self.assertTrue(rules[0]["backlash"])
+        self.assertAlmostEqual(float(st.P[2].sum()), 300.0)
+
+    def test_small_remote_cells_survive_big_exposed_ones_pay(self):
+        # cell 0: 10,000 Blod beside 50,000 enemies; cell 4: 100 Blod alone → survivors 30 % of 10,100 = 3,030
+        w, races, st, rules = self.rig(backlash_years=1)
+        st.P[2] = [10000, 0, 0, 0, 100]
+        st.P[0] = [0, 50000, 0, 0, 0]
+        st.t = 101
+        rituals(st, w, races, rules, np.random.default_rng(0))
+        self.assertAlmostEqual(float(st.P[2].sum()), 3030.0, places=6)
+        # kill share ∝ exposure (own + nearby enemies): 60,000 vs 100 → cell 4 loses 100·100·κ ≈ 0.12
+        kappa = 7070 / (10000 * 60000 + 100 * 100)
+        self.assertAlmostEqual(float(st.P[2, 4]), 100 - 100 * 100 * kappa, places=6)
+        self.assertGreater(float(st.P[2, 4]), 99.8)
diff --git a/tests/test_migration.py b/tests/test_migration.py
new file mode 100644
index 0000000..160f75a
--- /dev/null
+++ b/tests/test_migration.py
@@ -0,0 +1,200 @@
+import unittest
+
+import numpy as np
+
+from tests.helpers import line_world, make_history, make_race
+from worldhistory.migration import drift_and_bud, long_jumps, travel_cost
+from worldhistory.config import CONDITIONS
+from worldhistory.state import new_state
+
+D_DEPTH = CONDITIONS.index("depth")
+
+MP = make_history().migration
+
+
+def setup(n, pop, race, K=1e6):
+    st = new_state(1, n, np.zeros((1, 6)), seed=0)
+    st.P[0, 0] = pop
+    return st, np.full((1, n), K), np.ones((1, n), bool)
+
+
+class MigrationTest(unittest.TestCase):
+    def test_travel_cost(self):
+        w = line_world(4, ocean=[False, False, False, True], landform=[1, 3, 1, 0], river=[True, False, False, False],
+                       holdridge=[20, 20, 23, 0])
+        r = make_race(travel={"mountain": 2.0, "desert": 1.0})
+        np.testing.assert_allclose(travel_cost(w, r, np.zeros(4)), [-1, 2, 0, 0])   # cell 2: desert +1, coast -1
+        np.testing.assert_allclose(travel_cost(w, r, np.full(4, 0.5)), [-1.5, 1.5, -0.5, -0.5])
+
+    def test_budding_needs_founder_group(self):
+        w = line_world(3)
+        r = make_race(demography={"founder": 50.0, "mobility": 0.1})
+        st, K, OK = setup(3, 60.0, r, K=100.0)          # crowded, but below 2 founders
+        drift_and_bud(st, w, [r], K, OK, np.zeros((1, 3)), np.zeros((1, 3)), np.ones((1, 3)),
+                      np.random.default_rng(0), MP)
+        self.assertEqual(float(st.P[0, 1]), 0.0)
+        st, K, OK = setup(3, 120.0, r, K=100.0)
+        drift_and_bud(st, w, [r], K, OK, np.zeros((1, 3)), np.zeros((1, 3)), np.ones((1, 3)),
+                      np.random.default_rng(0), MP)
+        self.assertIn(float(st.P[0, 1]), (0.0, 20.0))     # a founder group (capped at P/6 = 20) or nothing
+        self.assertAlmostEqual(float(st.P[0].sum()), 120.0)
+
+    def test_river_colonised_before_mountain(self):
+        r = make_race(demography={"founder": 5.0, "mobility": 0.1}, travel={"mountain": 2.0})
+        hits = {0: 0, 2: 0}
+        for seed in range(60):
+            st = new_state(1, 3, np.zeros((1, 6)), seed=0)
+            st.P[0, 1] = 100.0
+            w2 = line_world(3, landform=[1, 1, 3], river=[True, True, False])   # cell 0 river plain, cell 2 mountain
+            tc = travel_cost(w2, r, np.zeros(3))[None]
+            drift_and_bud(st, w2, [r], np.full((1, 3), 100.0), np.ones((1, 3), bool), tc, np.zeros((1, 3)),
+                          np.ones((1, 3)), np.random.default_rng(seed), MP)
+            hits[0] += st.P[0, 0] > 0
+            hits[2] += st.P[0, 2] > 0
+        self.assertGreater(hits[0], hits[2])
+
+    def test_drift_toward_headroom(self):
+        w = line_world(3)
+        r = make_race(demography={"mobility": 0.2})
+        st = new_state(1, 3, np.zeros((1, 6)), seed=0)
+        st.P[0] = [100.0, 100.0, 100.0]
+        K = np.array([[1000.0, 1000.0, 5000.0]])
+        drift_and_bud(st, w, [r], K, np.ones((1, 3), bool), np.zeros((1, 3)), np.zeros((1, 3)),
+                      np.ones((1, 3)), np.random.default_rng(0), MP)
+        self.assertGreater(float(st.P[0, 2]), float(st.P[0, 0]))  # cell 1 sends more toward the roomier side
+        self.assertAlmostEqual(float(st.P[0].sum()), 300.0)
+
+    def test_never_into_uninhabitable(self):
+        w = line_world(3)
+        r = make_race(demography={"mobility": 0.2, "founder": 1.0})
+        st = new_state(1, 3, np.zeros((1, 6)), seed=0)
+        st.P[0] = [0.0, 100.0, 100.0]
+        drift_and_bud(st, w, [r], np.full((1, 3), 100.0), np.array([[False, True, True]]), np.zeros((1, 3)),
+                      np.zeros((1, 3)), np.ones((1, 3)), np.random.default_rng(0), MP)
+        self.assertEqual(float(st.P[0, 0]), 0.0)
+        self.assertAlmostEqual(float(st.P[0].sum()), 200.0)
+
+    def test_drift_keeps_viable_source_viable(self):
+        # uncrowded viable group (22, Allee 20) beside a small one: drift may not take it below the Allee threshold
+        w = line_world(2)
+        r = make_race(demography={"mobility": 0.5, "allee": 20.0, "founder": 20.0})
+        st = new_state(1, 2, np.zeros((1, 6)), seed=0)
+        st.P[0] = [22.0, 5.0]
+        drift_and_bud(st, w, [r], np.full((1, 2), 1e5), np.ones((1, 2), bool), np.zeros((1, 2)),
+                      np.zeros((1, 2)), np.ones((1, 2)), np.random.default_rng(0), MP)
+        self.assertGreaterEqual(float(st.P[0, 0]), 20.0)
+        self.assertAlmostEqual(float(st.P[0].sum()), 27.0)
+
+    def test_small_group_joins_viable_neighbour(self):
+        # below the Allee threshold, a group gathers into its viable neighbour rather than dwindling alone
+        w = line_world(3)
+        r = make_race(demography={"mobility": 0.1, "allee": 20.0, "founder": 20.0})
+        st = new_state(1, 3, np.zeros((1, 6)), seed=0)
+        st.P[0] = [8.0, 30.0, 0.0]
+        drift_and_bud(st, w, [r], np.full((1, 3), 1e5), np.ones((1, 3), bool), np.zeros((1, 3)),
+                      np.zeros((1, 3)), np.ones((1, 3)), np.random.default_rng(0), MP)
+        np.testing.assert_allclose(st.P[0], [0.0, 38.0, 0.0])
+
+    def test_push_flees_into_empty_land(self):
+        # people driven out (push) settle an empty livable neighbour even when nobody is crowded
+        w = line_world(2)
+        r = make_race(demography={"mobility": 0.0, "allee": 20.0, "founder": 20.0})
+        st = new_state(1, 2, np.zeros((1, 6)), seed=0)
+        st.P[0] = [120.0, 0.0]
+        st.push[0, 0] = 30.0
+        drift_and_bud(st, w, [r], np.full((1, 2), 1e5), np.ones((1, 2), bool), np.zeros((1, 2)),
+                      np.zeros((1, 2)), np.ones((1, 2)), np.random.default_rng(0), MP)
+        np.testing.assert_allclose(st.P[0], [100.0, 20.0])                  # capped at P/6 per neighbour
+
+    def test_roaming_groups_travel_together(self):
+        # roam = chance per step that a whole group (a caravan) moves on, together, to one neighbouring cell
+        w = line_world(3)
+        r = make_race(demography={"mobility": 0.0, "allee": 2.0, "founder": 2.0, "roam": 1.0})
+        ends = set()
+        for seed in range(20):
+            st = new_state(1, 3, np.zeros((1, 6)), seed=0)
+            st.P[0] = [0.0, 120.0, 0.0]
+            drift_and_bud(st, w, [r], np.full((1, 3), 1e5), np.ones((1, 3), bool), np.zeros((1, 3)),
+                          np.zeros((1, 3)), np.ones((1, 3)), np.random.default_rng(seed), MP)
+            self.assertEqual(float(st.P[0, 1]), 0.0)
+            self.assertIn(float(st.P[0].max()), (120.0,))
+            ends.add(int(st.P[0].argmax()))
+        self.assertEqual(ends, {0, 2})                                     # either way, never split
+
+    def test_long_jump_crosses_chasm_on_migrants_adaptation(self):
+        # deep-adapted swimmers skip an abyss they cannot live in and land on ground that suits *them*, even though
+        # the landing cells' stored optimum (native, 0 m) would call them unfit
+        from worldhistory.habitat import environment, suitability
+        depth = np.r_[np.full(5, 2000.0), np.full(35, 6000.0), np.full(40, 2000.0)]
+        w = line_world(80, spacing_km=10.0, ocean=True, elevation_m=-depth)
+        tol = {"depth": {"optimum": 0, "width": 60, "lo": 0, "hi": 7000, "rate": 2.0}}
+        r = make_race(demography={"founder": 2.0, "allee": 1.0, "jump": 1.0}, tolerance=tol,
+                      habitat={"realm": "sea", "terms": [{"p": "sea", "w": 1.0}]})
+        env, suit = environment(w), suitability(w, r)[None]
+        landed = set()
+        for seed in range(200):
+            st = new_state(1, w.n, np.zeros((1, 6)), seed=0)
+            st.P[0, 0] = 1e6
+            st.O[0, D_DEPTH, :5] = 2000.0
+            long_jumps(st, w, [r], np.full((1, w.n), 1e6), np.ones((1, w.n), bool), np.random.default_rng(seed),
+                       {**MP, "jump_rate": 1.0}, np.ones((1, w.n)), env=env, suit=suit, f_min=0.1)
+            landed |= set(np.flatnonzero(st.P[0, 1:] > 0) + 1)
+        self.assertTrue(any(c >= 40 for c in landed))                    # across the chasm
+        self.assertFalse(any(5 <= c < 40 for c in landed))               # never into the abyss
+
+    def test_sea_path_jumps_never_cross_land(self):
+        # swimmers (jump_path = sea) cannot hop over a strip of land; without it they can
+        ocean = np.r_[np.ones(5, bool), np.zeros(5, bool), np.ones(70, bool)]
+        w = line_world(80, spacing_km=10.0, ocean=ocean, elevation_m=np.where(ocean, -100.0, 10.0))
+        from worldhistory.habitat import environment, suitability
+        for path, beyond in (("any", True), ("sea", False)):
+            r = make_race(demography={"founder": 2.0, "allee": 1.0, "jump_path": path},
+                          habitat={"realm": "sea", "terms": [{"p": "sea", "w": 1.0}]})
+            env, suit = environment(w), suitability(w, r)[None]
+            landed = set()
+            for seed in range(200):
+                st = new_state(1, w.n, np.zeros((1, 6)), seed=0)
+                st.P[0, 0] = 1e6
+                long_jumps(st, w, [r], np.full((1, w.n), 1e6), np.ones((1, w.n), bool), np.random.default_rng(seed),
+                           {**MP, "jump_rate": 1.0}, np.ones((1, w.n)), env=env, suit=suit, f_min=0.1)
+                landed |= set(np.flatnonzero(st.P[0, 1:] > 0) + 1)
+            self.assertEqual(any(c >= 10 for c in landed), beyond, path)
+
+    def test_large_group_splits(self):
+        # clans above group_max break up: half leaves as a new group into empty land
+        w = line_world(3)
+        r = make_race(demography={"mobility": 0.0, "allee": 20.0, "founder": 20.0, "group_max": 500.0})
+        st = new_state(1, 3, np.zeros((1, 6)), seed=0)
+        st.P[0] = [0.0, 1200.0, 0.0]
+        drift_and_bud(st, w, [r], np.full((1, 3), 1e5), np.ones((1, 3), bool), np.zeros((1, 3)),
+                      np.zeros((1, 3)), np.ones((1, 3)), np.random.default_rng(0), MP)
+        np.testing.assert_allclose(st.P[0], [200.0, 800.0, 200.0])          # 600 leave, capped P/6 per side
+
+    def test_push_leaves(self):
+        w = line_world(2)
+        r = make_race(demography={"mobility": 0.0})
+        st = new_state(1, 2, np.zeros((1, 6)), seed=0)
+        st.P[0] = [100.0, 100.0]
+        st.push[0, 0] = 10.0
+        drift_and_bud(st, w, [r], np.full((1, 2), 1000.0), np.ones((1, 2), bool), np.zeros((1, 2)),
+                      np.zeros((1, 2)), np.ones((1, 2)), np.random.default_rng(0), MP)
+        np.testing.assert_allclose(st.P[0], [90.0, 110.0])
+        self.assertEqual(float(st.push.sum()), 0.0)
+
+    def test_long_jump_distances(self):
+        w = line_world(4000, spacing_km=10.0)
+        r = make_race(demography={"founder": 2.0, "allee": 1.0})
+        st = new_state(1, w.n, np.zeros((1, 6)), seed=0)
+        st.P[0, 0] = 1e6
+        d = []
+        for seed in range(300):
+            st.P[0] = 0
+            st.P[0, 0] = 1e6
+            long_jumps(st, w, [r], np.full((1, w.n), 1e6), np.ones((1, w.n), bool), np.random.default_rng(seed),
+                       {**MP, "jump_rate": 1.0}, np.ones((1, w.n)))
+            hit = np.flatnonzero(st.P[0, 1:]) + 1
+            d += list(np.minimum(hit, w.n - hit) * 10.0)          # the line circles the globe: wrap
+        d = np.array(d)
+        self.assertGreater(len(d), 100)          # a line world: only jumps along the line land near it
+        self.assertTrue((d <= MP["jump_cap_km"] + 10).all())
+        self.assertGreater(float(np.median(d)), MP["jump_xm_km"] * 0.8)
diff --git a/tests/test_predicates.py b/tests/test_predicates.py
new file mode 100644
index 0000000..542b67f
--- /dev/null
+++ b/tests/test_predicates.py
@@ -0,0 +1,80 @@
+import unittest
+
+import numpy as np
+
+from tests.helpers import fields, line_world
+from worldhistory.predicates import check, depth, evaluate, parse
+from worldhistory.regions import Regions, region_mask
+
+
+class PredicateTest(unittest.TestCase):
+    def setUp(self):
+        # cells: 0 forest land, 1 desert land by a river, 2 mountain land, 3 shelf sea, 4 deep sea with a vent
+        self.w = line_world(5, ocean=[False, False, False, True, True], holdridge=[20, 23, 12, 0, 0],
+                            landform=[1, 1, 3, 0, 0], river=[False, True, False, False, False],
+                            strahler=[0, 4, 0, 0, 0], elevation_m=[100, 50, 2000, -150, -3000],
+                            vent_potential=[0, 0, 0, 0, 0.8], T_mean=[10, 30, 0, 12, 5], lithology=[1, 3, 4, 0, 0])
+
+    def test_parse(self):
+        self.assertEqual(parse("land"), ("land", []))
+        self.assertEqual(parse("warm:5:20"), ("warm", [5.0, 20.0]))
+        self.assertEqual(parse("lithology:3,4"), ("lithology", [[3, 4]]))
+        self.assertEqual(parse("above:o2_fraction:0.3"), ("above", ["o2_fraction", 0.3]))
+
+    def test_basic(self):
+        e = lambda p: evaluate(p, self.w).tolist()
+        self.assertEqual(e("land"), [1, 1, 1, 0, 0])
+        self.assertEqual(e("forest"), [1, 0, 0, 0, 0])
+        self.assertEqual(e("desert"), [0, 1, 0, 0, 0])
+        self.assertEqual(e("mountain"), [0, 0, 1, 0, 0])
+        self.assertEqual(e("coast"), [0, 0, 1, 0, 0])
+        self.assertEqual(e("shelf"), [0, 0, 0, 1, 0])
+        self.assertEqual(e("deep_sea"), [0, 0, 0, 0, 1])
+        self.assertAlmostEqual(e("water")[1], 1.0)              # 0.4 + 0.15*4 = 1.0
+        self.assertAlmostEqual(e("vent")[4], 0.8)
+        self.assertEqual(e("lithology:3,4"), [0, 1, 1, 0, 0])
+
+    def test_ranges_and_products(self):
+        self.assertEqual(evaluate("warm:5:25", self.w).tolist(), [0.25, 1.0, 0.0, 0.35, 0.0])
+        self.assertEqual(evaluate(["land", "above:T_mean:20"], self.w).tolist(), [0, 1, 0, 0, 0])
+        np.testing.assert_array_equal(depth(self.w.fields), [0, 0, 0, 150, 3000])
+
+    def test_check(self):
+        check(["land", "warm:1:2", "lithology:1"])
+        for bad in ("nosuch", "warm:1", "above:nofield:1"):
+            with self.assertRaises(ValueError):
+                check(bad, field_names=set(fields(1)))
+
+
+class RegionTest(unittest.TestCase):
+    def test_kinds(self):
+        w = line_world(4, ocean=[False, False, False, True], plate=[1, 1, 2, 2], gravity_g=1.0)
+        w.eras["late"] = fields(4, ocean=[False, True, False, True], gravity_g=[1.0, 0.35, 0.35, 0.35])
+        self.assertEqual(region_mask(w, {"kind": "all"}).sum(), 4)
+        self.assertEqual(region_mask(w, {"kind": "plate", "ids": [2]}).tolist(), [False, False, True, True])
+        self.assertEqual(region_mask(w, {"kind": "changed", "era": "late", "fields": ["gravity_g"]}).tolist(),
+                         [False, True, True, True])
+        self.assertEqual(region_mask(w, {"kind": "changed", "era": "late", "fields": ["gravity_g"], "land": "base"})
+                         .tolist(), [False, True, True, False])
+        box = region_mask(w, {"kind": "box", "lat": [-1, 1], "lon": [0.5, 10]})
+        self.assertEqual(box.tolist(), [False, True, True, True])
+        self.assertEqual(region_mask(w, {"kind": "predicate", "p": "land"}).tolist(), [True, True, True, False])
+
+    def test_box_with_predicate_filter(self):
+        w = line_world(4, ocean=[False, False, True, True])
+        m = region_mask(w, {"kind": "box", "lat": [-1, 1], "lon": [0.5, 10], "p": "land"})
+        self.assertEqual(m.tolist(), [False, True, False, False])
+
+    def test_exclude_other_region(self):
+        w = line_world(4)
+        r = Regions(w, {"east": {"kind": "box", "lat": [-1, 1], "lon": [0.5, 10]},
+                        "land_not_east": {"kind": "predicate", "p": "land", "exclude": "east"}})
+        self.assertEqual(r("land_not_east").tolist(), [True, False, False, False])
+
+    def test_registry(self):
+        w = line_world(3)
+        r = Regions(w, {"west": {"kind": "box", "lat": [-1, 1], "lon": [-1, 0.5]}})
+        self.assertEqual(r("west").tolist(), [True, False, False])
+        self.assertTrue(r("all").all())
+        with self.assertRaises(KeyError):
+            r("nowhere")
diff --git a/tests/test_state.py b/tests/test_state.py
new file mode 100644
index 0000000..6ed8351
--- /dev/null
+++ b/tests/test_state.py
@@ -0,0 +1,63 @@
+import unittest
+
+import numpy as np
+
+from worldhistory.state import convert, move, new_state
+
+
+def st3():
+    st = new_state(2, 3, np.array([[0.0, 1.0, 0, 0, 0, 0], [5.0, 1.0, 0, 0, 0, 0]]), seed=0)
+    st.P[0] = [10.0, 0.0, 30.0]
+    st.O[0, 0] = [100.0, 999.0, 300.0]       # cell 1 is empty with a stale value
+    return st
+
+
+class StateTest(unittest.TestCase):
+    def test_natives_need_one_column_per_condition(self):
+        from worldhistory.config import CONDITIONS
+        with self.assertRaisesRegex(ValueError, "conditions"):
+            new_state(1, 3, np.zeros((1, len(CONDITIONS) - 1)), seed=0)
+
+    def test_shapes(self):
+        st = new_state(2, 3, np.zeros((2, 6)), seed=0)
+        self.assertEqual(st.P.shape, (2, 3))
+        self.assertEqual(st.O.shape, (2, 6, 3))
+        self.assertFalse(hasattr(st, "E"))
+        self.assertEqual(st.T.shape, (2, 6, 3))
+
+    def test_move_conserves_and_mixes(self):
+        st = st3()
+        move(st, 0, [0], [2], [10.0])
+        np.testing.assert_allclose(st.P[0], [0, 0, 40])
+        self.assertAlmostEqual(float(st.O[0, 0, 2]), (30 * 300 + 10 * 100) / 40)
+
+    def test_move_into_empty_takes_arrivals_attrs(self):
+        st = st3()
+        move(st, 0, [0], [1], [4.0])
+        self.assertAlmostEqual(float(st.O[0, 0, 1]), 100.0)
+
+    def test_move_caps_at_available(self):
+        st = st3()
+        move(st, 0, [0, 0], [1, 2], [15.0, 5.0])          # asks 20 of 10: scaled to 7.5 + 2.5
+        np.testing.assert_allclose(st.P[0], [0, 7.5, 32.5])
+        self.assertAlmostEqual(float(st.P[0].sum()), 40.0)
+
+    def test_convert_between_slots(self):
+        st = st3()
+        st.P[1, 2] = 10.0
+        st.O[1, 0, 2] = 0.0
+        convert(st, 0, 1, np.array([2]), np.array([10.0]))
+        self.assertAlmostEqual(float(st.P[0, 2]), 20.0)
+        self.assertAlmostEqual(float(st.P[1, 2]), 20.0)
+        self.assertAlmostEqual(float(st.O[1, 0, 2]), 150.0)     # (10*0 + 10*300) / 20
+
+    def test_convert_with_new_optima(self):
+        st = new_state(2, 3, np.zeros((2, 6)), seed=0)
+        st.P[0, 2] = 30.0
+        st.P[1, 2] = 10.0
+        st.O[1, 0, 2] = 100.0
+        O_new = np.zeros((6, 1))
+        O_new[0, 0] = 500.0
+        convert(st, 0, 1, [2], [30.0], O_new=O_new)
+        self.assertAlmostEqual(float(st.O[1, 0, 2]), (10 * 100 + 30 * 500) / 40)
+        self.assertEqual(float(st.O[0, 0, 2]), 0.0)                  # the parents' optima are untouched
diff --git a/tests/test_tech.py b/tests/test_tech.py
new file mode 100644
index 0000000..4339c4a
--- /dev/null
+++ b/tests/test_tech.py
@@ -0,0 +1,64 @@
+import unittest
+
+import numpy as np
+
+from tests.helpers import line_world, make_history, make_race
+from worldhistory.state import new_state
+from worldhistory.tech import DI, effective, reach, step_land, step_tech
+
+TP = make_history().tech
+
+
+def state(n, pops):
+    st = new_state(1, n, np.zeros((1, 6)), seed=0)
+    st.P[0] = pops
+    return st
+
+
+class TechTest(unittest.TestCase):
+    def test_progress_grows_with_population(self):
+        w = line_world(7)
+        st = state(7, [0, 0, 10000, 0, 0, 20, 0])
+        step_tech(st, w, [make_race()], {**TP, "diffuse": 0.0, "loss_below": 0.0})
+        big, small = st.T[0, DI["farming"], 2], st.T[0, DI["farming"], 5]
+        # s = N/(N+n_half): 10000/15000 vs 20/5020; gain = rate*s
+        self.assertAlmostEqual(float(big), TP["rate"] * 10000 / 15000, places=6)
+        self.assertAlmostEqual(float(small), TP["rate"] * 20 / 5020, places=6)
+        self.assertEqual(float(st.T[0, 0, 0]), 0.0)                    # empty cells do not learn
+
+    def test_priority_scales(self):
+        w = line_world(1)
+        st = state(1, [10000])
+        step_tech(st, w, [make_race(tech={"travel": 0.0, "farming": 2.0})], {**TP, "diffuse": 0.0})
+        self.assertEqual(float(st.T[0, DI["travel"], 0]), 0.0)
+        self.assertAlmostEqual(float(st.T[0, DI["farming"], 0]), 2 * TP["rate"] * 10000 / 15000, places=6)
+
+    def test_small_isolated_groups_lose_tech(self):
+        w = line_world(1)
+        st = state(1, [50])
+        st.T[0, :, 0] = 0.8
+        step_tech(st, w, [make_race()], {**TP, "diffuse": 0.0, "rate": 0.0})
+        self.assertAlmostEqual(float(st.T[0, 0, 0]), 0.8 * (1 - TP["loss_rate"]), places=6)
+
+    def test_diffusion_only_upward(self):
+        w = line_world(2)
+        st = state(2, [1000, 1000])
+        st.T[0, 0] = [1.0, 0.0]
+        step_tech(st, w, [make_race()], {**TP, "rate": 0.0, "loss_below": 0.0, "diffuse": 0.5})
+        self.assertAlmostEqual(float(st.T[0, 0, 0]), 1.0)
+        self.assertAlmostEqual(float(st.T[0, 0, 1]), 0.25)            # 0.5 * (mean 0.5 - 0)
+
+    def test_magic_folds_in_and_reach_capped(self):
+        T = np.zeros((6, 1))
+        T[DI["magic"]] = 1.0
+        eff = effective(T, TP)
+        self.assertAlmostEqual(float(eff["farming"][0]), TP["magic_share"])
+        self.assertAlmostEqual(float(reach(eff, make_race(tech={"reach_cap": 0.1}), TP)[0]), 0.1)
+
+    def test_land_improvement_up_and_decay(self):
+        st = state(2, [100, 0])
+        st.T[0, DI["land"]] = 1.0
+        st.improve[0] = [0.0, 0.4]
+        step_land(st, [make_race()], TP)
+        self.assertAlmostEqual(float(st.improve[0, 0]), TP["land_rate"] * 1.0 * TP["land_cap"], places=6)
+        self.assertAlmostEqual(float(st.improve[0, 1]), 0.4 * (1 - TP["land_decay"]), places=6)
diff --git a/tests/test_world.py b/tests/test_world.py
new file mode 100644
index 0000000..a9969f8
--- /dev/null
+++ b/tests/test_world.py
@@ -0,0 +1,84 @@
+import json
+import tempfile
+import unittest
+from pathlib import Path
+
+import numpy as np
+
+from tests.helpers import fields, globe_world, line_world
+from worldhistory.world import FIELDS, load_world
+
+
+class WorldTest(unittest.TestCase):
+    def test_line_neighbours_and_sum(self):
+        w = line_world(4)
+        np.testing.assert_array_equal(w.nb_sum(np.array([1.0, 2.0, 3.0, 4.0])), [2, 4, 6, 3])
+        np.testing.assert_array_equal(w.nb_any(np.array([True, False, False, False])), [False, True, False, False])
+
+    def test_nb_sum_on_stacked_arrays(self):
+        w = line_world(3)
+        x = np.array([[1.0, 0, 0], [0, 0, 5.0]])
+        np.testing.assert_array_equal(w.nb_sum(x), [[0, 1, 0], [0, 5, 0]])
+
+    def test_smooth_keeps_constant(self):
+        w = globe_world(1)
+        np.testing.assert_allclose(w.smooth(np.full(w.n, 3.0), 2), 3.0)
+
+    def test_noise_standardised_and_seeded(self):
+        w = globe_world(1)
+        a, b = w.noise(5), w.noise(5)
+        np.testing.assert_array_equal(a, b)
+        self.assertAlmostEqual(float(a.mean()), 0.0, places=6)
+        self.assertAlmostEqual(float(a.std()), 1.0, places=6)
+        self.assertFalse(np.array_equal(a, w.noise(6)))
+
+    def test_globe_neighbours_are_close(self):
+        w = globe_world(1)
+        valid = w.nb >= 0
+        self.assertTrue((valid.sum(1) >= 5).all())       # hexagons 6, pentagons 5
+        i = np.repeat(np.arange(w.n), 6)[valid.ravel()]
+        j = w.nb.ravel()[valid.ravel()]
+        self.assertLess(float(w.km(i, j).max()), 1000.0)  # res-1 cells are ~400-600 km apart on an Earth-size globe
+
+    def test_km_and_within(self):
+        w = line_world(10, spacing_km=100.0)
+        self.assertAlmostEqual(float(w.km(np.array([0]), np.array([3]))[0]), 300.0, places=3)
+        self.assertEqual(sorted(w.within(5, 150.0).tolist()), [4, 5, 6])
+
+    def test_era_switch(self):
+        w = line_world(3)
+        w.eras["later"] = fields(3, ocean=True)
+        self.assertFalse(w.fields["ocean"].any())
+        w.set_era("later")
+        self.assertTrue(w.fields["ocean"].all())
+        with self.assertRaises(KeyError):
+            w.set_era("nope")
+
+    def test_load_world_roundtrip(self):
+        g = globe_world(1)
+        import h3.api.basic_int as h3
+        ids = np.array(sorted(c for r0 in h3.get_res0_cells() for c in h3.cell_to_children(r0, 1)), np.uint64)
+        with tempfile.TemporaryDirectory() as d:
+            base = {f"g_{k}": v for k, v in dict(ids=ids, lat=g.lat, lon=g.lon, xyz=g.xyz * 2.0,
+                                                  area_km2=g.area).items()}
+            np.savez(Path(d) / "cells.npz", **base, **fields(g.n))
+            (Path(d) / "eras" / "late").mkdir(parents=True)
+            np.savez(Path(d) / "eras" / "late" / "cells.npz", **fields(g.n, ocean=True))
+            (Path(d) / "cells_meta.json").write_text(json.dumps({"radius_km": 12742.0}))
+            w = load_world(d, eras=["late"])
+        self.assertEqual(w.radius_km, 12742.0)
+        np.testing.assert_allclose(np.linalg.norm(w.xyz, axis=1), 1.0)
+        np.testing.assert_array_equal(w.nb, g.nb)
+        self.assertTrue(set(FIELDS) <= set(w.base))
+        self.assertTrue(w.eras["late"]["ocean"].all())
+
+    def test_load_world_missing_field(self):
+        g = globe_world(1)
+        import h3.api.basic_int as h3
+        ids = np.array(sorted(c for r0 in h3.get_res0_cells() for c in h3.cell_to_children(r0, 1)), np.uint64)
+        with tempfile.TemporaryDirectory() as d:
+            f = fields(g.n)
+            del f["gravity_g"]
+            np.savez(Path(d) / "cells.npz", g_ids=ids, g_lat=g.lat, g_lon=g.lon, g_xyz=g.xyz, g_area_km2=g.area, **f)
+            with self.assertRaisesRegex(ValueError, "gravity_g"):
+                load_world(d)