diff options
| author | godosa <godosa@godosa.eu> | 2026-10-07 00:01:14 +0200 |
|---|---|---|
| committer | godosa <godosa@godosa.eu> | 2026-10-07 00:01:14 +0200 |
| commit | ed1dea2639b1191421de3986483aedcc14067a12 (patch) | |
| tree | 0118c6e119a84f1a9433ed0a304c8be1c9098463 | |
| download | worldhistory-ed1dea2639b1191421de3986483aedcc14067a12.tar.gz worldhistory-ed1dea2639b1191421de3986483aedcc14067a12.zip | |
worldhistory: initial public history
48 files changed, 5365 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..20d77e4 --- /dev/null +++ b/.gitignore @@ -0,0 +1,20 @@ +__pycache__/ +*.pyc +.venv +.worktrees/ +out/ +node_modules/ + +# private state (lives in the private project repo) +TASKS.md +tasks/ +workflow.toml +docs/plans/ +docs/specs/ +docs/AREA-NOTES.md +USER-NEXT.md +CATALOG.md +publish-denylist.local +archive/ +.wf/ +.superpowers/ diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..519ea94 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,26 @@ +# CLAUDE.md — worldhistory + +History simulation on a worldgen world (races spread, adapt, split, fight, invent; eras and events). Terse everywhere. +Private notes (tasks, specs, area maps) are kept outside this repo; never add them here (see `.gitignore`). + +## Layout +- `history.py` — CLI entry (→ `worldhistory/cli.py`: run, preview, compare, check-config) +- `worldhistory/` — engine package (one step = `phase_*` in engine.py, README order); `tests/` — unittest, + synthetic worlds (`tests/helpers.py`: `line_world` / `fields`, 1-D, hand-checkable) +- `docs/reference/` — research notes (canals) + +## Build / test +- Setup: `python3 -m venv .venv && .venv/bin/pip install -r requirements.txt`; in a git worktree link the main + tree's `.venv` (`ln -sfn <main>/.venv .venv`; `.venv` is ignored without trailing slash: it may be a link). +- Verify: `.venv/bin/python -m unittest discover -s tests -t .` (whole suite ~3 s). + +## Gotchas +- System python lacks scipy: always `.venv/bin/python`. +- numba optional (not in requirements): `state.py` JIT fast paths; `tests/test_fastpaths.py` skips without it. +- Input = a worldgen build (`cells.npz` + `eras/<era>/cells.npz`). A full run on a large world is big (README: + < 2 GB at r4 ≈ 290k cells). `SmokeTest` (test_engine.py) runs 100 years when `WORLDHISTORY_SMOKE_WORLD` + + `WORLDHISTORY_SMOKE_CONFIG` are set. +- Config: unknown keys are errors (config.py `_keys`); `history.py check-config <dir> [--world <build>]` validates. + +## Git +- Work in a git worktree on a branch (`.worktrees/<topic>`), not in the main tree; ff-merge when green. diff --git a/CREDITS.md b/CREDITS.md new file mode 100644 index 0000000..f0ff82d --- /dev/null +++ b/CREDITS.md @@ -0,0 +1,27 @@ +# Credits + +| Kind | Who | Note | +|---|---|---| +| Author | godosa | | +| Libraries | see [licenses.md](licenses.md) | | +| Tooling | AI assistance (Claude, Anthropic) | engineering tool; statement in README | + +## Models and sources + +| Source | Used for | Module | +|---|---|---| +| Verhulst 1838 (logistic growth) | growth toward capacity | `demography.py` | +| Allee 1931; Courchamp, Clutton-Brock & Grenfell 1999 (Allee effects) | small groups dwindle | `demography.py` | +| Fisher 1937; Kolmogorov, Petrovsky & Piskunov 1937 (Fisher–KPP wave of advance) | spread as a travelling front | `migration.py` | +| Ammerman & Cavalli-Sforza 1971, 1984 (*The Neolithic Transition and the Genetics of Populations in Europe*) | wave-of-advance speed (~1 km/yr) as a calibration reference | `migration.py` | +| Kot, Lewis & van den Driessche 1996; Clark 1998 (fat-tailed dispersal) | long-distance founder jumps | `migration.py` | +| Fretwell & Lucas 1970 (ideal free distribution); McFadden 1974 (logit choice) | drift toward headroom | `migration.py` | +| Lotka 1925; Volterra 1926 (competition) | races sharing capacity | `demography.py` | +| Kremer 1993, QJE 108(3) (population and technological change); Henrich 2004, American Antiquity 69(2) (demography and cultural evolution) | tech grows with connected population, lost in small isolated groups | `tech.py` | +| Holdridge 1947 (life zones; via worldgen) | habitat predicates | `predicates.py` | +| H3 hexagonal grid (Uber, h3geo.org) | cells and neighbours | `world.py` | +| Background: Epstein & Axtell 1996 (Sugarscape); Axtell et al. 2002, PNAS (Long House Valley) | agent-based comparison in the design spike | — | +| Real canals (Wikipedia articles on Suez, Panama, Kiel, Corinth, Göta, Caledonian, Erie, Welland, Canal du Midi, Grand Canal, Lingqu, Stecknitz, Briare, Fossa Carolina, Diolkos, Xerxes, Canal of the Pharaohs, Volga–Don, White Sea–Baltic, Rhine–Main–Danube, Bridgewater, Manchester Ship Canal) | calibration reference for great works: `docs/reference/canals.md` (planned: sub-project 2) | — | +| Lanchester 1916 | not used yet (planned: war, sub-project 2) | — | +| Turchin 2003; Turchin, Currie, Turner & Gavrilets 2013, PNAS 110(41) | not used yet (planned: polities) | — | +| Wilson 1967/1970; Christaller 1933 | not used yet (planned: trade, settlements) | — | @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2026 godosa + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..e99e612 --- /dev/null +++ b/README.md @@ -0,0 +1,116 @@ +# worldhistory + +Trials of how peoples spread, settle, adapt and change over thousands of years on a planet made with +[worldgen](../worldgen). Each people ("race") is a population field on worldgen's hexagonal cells: it grows toward +what the land can carry, competes with its neighbours, migrates (drift, founder groups, rare long jumps), adapts to +depth, gravity, air and temperature, can split into sub-species, and learns (or loses) tech — magic included. +Scheduled events (culls, era switches) and hazards (roaming monsters) shape it. Results are **trials of +algorithms**, not a history to take as given: you tune the config and compare runs. + +## Requirements + +- Python ≥ 3.11, `pip install -r requirements.txt` (numpy, scipy, pillow, h3) +- A worldgen build (`out/r<res>/cells.npz`, optionally `eras/<era>/cells.npz`) +- RAM: < 2 GB at res 4 (~290k cells) for the default races; about 10 minutes for 7,500 years. + +## Quick start + +```bash +python3 -m venv .venv && .venv/bin/pip install -r requirements.txt +.venv/bin/python history.py check-config ../myworld/config/history --world ../myworld/out/r4 +.venv/bin/python history.py run --world ../myworld/out/r4 --config ../myworld/config/history --trial t1 --seed 1 +``` + +## How it works + +One step = 10 years, in this order: + +1. **events** — scheduled: `seed` (first placement), `cull` / `die_off` (a share of people in a region), `era_switch` +2. **hazards** — static mortality layers; roaming units that raid, breed in clutches and die +3. **capacity** — habitat quality × density × tech gain (∝ q²: poor land stays in low gear) × harvest shocks × + tolerance fitness × comfort × nudges +4. **growth** — logistic with an Allee effect; races share capacity (Lotka–Volterra); optional "veins" +5. **conflict** — background war losses by temperament; yielding turns losses into emigration +6. **migration** — drift toward headroom, founder groups into empty land, fat-tailed long jumps; with + `[demography] current_bias` > 0 (default 0; e.g. 1–3 for swimmers) sea jumps and drift between sea cells favour + going with the ocean current +7. **adaptation** — each condition's optimum moves toward where people live; the widths are fixed, so the old is lost; + past the race's adaptable range [`lo`, `hi`] drift is braked (fuzzy, not a line) +8. **sub-species** — a group's progress toward a related race is read from its adapted optima (weighted over the + conditions where the two races' native curves differ); a birth (once per sub-species) resets the group to the new + race's curves; displacement squeezes the less fit of the pair in shared cells; parent groups adapted past `mix_min` + bear the dominant race; adapted-back groups give **returns** (parent-race births, any number); `mode = "ritual"` + births by chance and wins converts at a blood price +9. **tech** — domains (farming, crafts, travel, magic, war, land) grow with connected population, lost when isolated + +## Config + +A config dir holds `history.toml` (run length, regions, events, hazards, nudges, tech and +migration parameters) and `races/<race>.toml` (demography, seed, habitat terms and gates over named predicates, +travel costs, tolerance envelopes, temperament, conflict, tech priorities, sub-species `emerge` rules). Unknown keys +are errors. + +Tolerances: `optimum`, `width` (or `width_lo` / `width_hi` for a lopsided bell), `lo`, `hi` (adaptable range, soft), +`rate`, `comfort`. Conditions: `depth`, `air_pressure`, `gravity`, `o2`, `temperature`, `rain` (log10 mm/yr; no +effect at sea). Strain: an adapted curve's peak is 0.35^(x²), x = (optimum − native) / (range edge − native) — adapted +groups are livable, not good; a `comfort` list ([[value, peak], …]) replaces this for its condition. Displacement +compares fit × comfort. A sub-species inherits the parent's curves per condition and overrides the ones that define it. + +`[emerge]` (defaults): `parent`, `mode` ("adapt" | "ritual"), `birth_min` 0.4, `birth_sure` 0.8, `birth_rate` 0.2 +(chance per group per step at `birth_sure`, × P/(P+founder)), `mix_min` 0.3, `dominance` 1, `birth_shift` 1, +`displace` 0.9, `return_min` / `return_sure` (= birth values), `convert` 1, `region`, `condition`, +`settlement_density`, `isolated`, `spread` ("contact" | "region"); ritual only: `trigger` 0.01, `ritual_rate` 0.02 +(rituals per step per √person), `ritual_cost` 300 deaths, `ritual_converts` 30, `victims` ("parent" | "all": +the dead come from every nearby non-ritual people by presence), `crisis_drop` 0 (off; > 0: birth only where parents +fell to ≤ (1 − drop) of their recent peak, `crisis_min` 1000 people, memory `crisis_years` 200; `crisis_km` 0: the survivors are born +into the nearest cell within that range where the new race can live, 0 = the crisis cell itself must fit), +`backlash_victims` 0 +(off; after that many dead in all, one smackdown: the race loses `backlash_loss` 0.7 and holds rituals at +`backlash_calm` 0.2 × the rate), `raid_rate` 0 (off; before the backlash, groups of ≥ `raid_min` 1000 send +~Poisson(rate) cells of `raid_size` 100 to occupied non-ritual land within `raid_km` 300; stats `raids`), `ritual_power` 0.5 (rituals ~ Poisson(rate · P^power); 1 = every cultist sacrifices), +`backlash_years` 0 (off; the smackdown also fires this many years after the birth). The smackdown's losses fall by +exposure (own people + non-ritual people in and next to the cell): big exposed groups are wiped out, small remote +ones survive. Any mode: `after` 0 (earliest birth year), `return_rate` (= `birth_rate`). The old `class` / `years` keys and +`[exposure]` are gone (errors name the replacement). Full configs live in the world project that runs the sim (e.g. its `map/config/history/`). + +Ocean fields from newer worldgen builds (`current`, `productivity`, `upwelling`, `sst`) load when present; older +builds behave as still water. The habitat predicate `productivity` (0–1 sea life) can weight a race's terms. +`World.sail_cost(ship_ms)` gives hours per directed sea edge (still-water ship speed helped or slowed by the +current) for route work; the engine does not use it yet. + +## CLI + +``` +history.py run --world DIR --config DIR --trial NAME [--seed N] [--years N] [--out DIR] [--no-preview] +history.py preview TRIAL_DIR # PNGs into TRIAL_DIR/preview/ +history.py compare TRIAL_A TRIAL_B [--out DIR] # seed agreement: same dominant race per occupied cell +history.py check-config CONFIG_DIR [--world DIR] +``` + +Long runs (> 10 min or > 2 GB): reserve and cap via the shared ledger, e.g. +`wf res run --mem 4G --for 30m --title "history t1" -- python history.py run …` (busy → exit 3, retry later). + +## Output + +`<world>/history/<trial>/`: `run.toml` (everything resolved, seed, engine version), `snap/y<year>.npz` (sparse +population, adapted optima, tech per race and cell), `stats.json` (per-step totals, regions, progress toward each sub-species, events, emergence, returns, rituals, +hazards, timing), `preview/*.png`. + +## Layout + +``` +history.py CLI +worldhistory/ engine modules (one per phase; engine.py drives them) +tests/ python -m unittest discover -s tests -t . +``` + +## Licence and credits + +MIT, see [LICENSE](LICENSE). Third-party notices: [licenses.md](licenses.md). Credits and the models used: +[CREDITS.md](CREDITS.md). + +## From the author + +These projects are things that have been tumbling about in my head for a long time, and I now feel like trying to do something with. The bulk of my focus here is on some of the games that I love, but every time I boot them up, I just fiddle about in the menu and set up mods and fixes for a couple hours, with my drive to play them fizzling out. This is my hope to solve that, and while I am a professional software engineer this would not have happened were it not for the rise of (relatively) cheap AI that could do the bulk of the work with me. If you don't approve, that's fine. I did these things for me, sharing them is something I do in the hopes to help others in similar situations, that just want to play the games they love. Thank you to all that made these playable in the first place, with some luck this finds you, and can bring some joy. + +Written with AI assistance (Claude, by Anthropic), used as an engineering tool. diff --git a/docs/reference/canals.md b/docs/reference/canals.md new file mode 100644 index 0000000..738f9ee --- /dev/null +++ b/docs/reference/canals.md @@ -0,0 +1,102 @@ +# Real canals: a calibration reference for great works + +Purpose: ground the sim's canal model ([[t-great-works]], [[t-canal-site-scan]]) in real projects. It covers what +canals were for, what they cost, who could build them, how long they took, which failed and what made them decay. +Gathered 2026-10-03 at the author's request. + +**Figures are approximate.** They come from the standard references (mostly each canal's Wikipedia article and the +works it cites), recalled rather than re-checked line by line. Lengths are the navigable route; "dug" is the +artificial part where a canal also uses lakes or rivers. Use them for orders of magnitude, not as data. + +## 1. Pre-industrial (closest to High Fantasy technology) + +| Canal | Built | Length | Summit / locks | Joins | Purpose | Builder, labour | +|---|---|---|---|---|---|---| +| Canal of the Pharaohs (Nile–Red Sea) | begun c. 600 BC (Necho II), finished c. 500 BC (Darius I) | ≈ 150–200 km | near sea level | Nile delta ↔ Red Sea | trade, prestige | Egyptian, then Persian state; Herodotus claims 120,000 dead under Necho | +| Xerxes Canal (Athos) | 483–480 BC | ≈ 2 km, ≈ 30 m wide | sea level | across the Athos peninsula | military: the fleet avoids the cape where a Persian fleet was wrecked in 492 BC | Persian army, ≈ 3 years | +| Diolkos (Corinth) — *portage, not a canal* | c. 600 BC, used to the 1st c. AD | ≈ 6–8.5 km paved haul road | over the isthmus | Gulf of Corinth ↔ Saronic Gulf | ships or cargo dragged across instead of sailing round the Peloponnese | city of Corinth; a cheap alternative used for centuries | +| Lingqu (Qin) | 214 BC | ≈ 36 km | contour canal over a watershed; flash locks later | Xiang (Yangtze basin) ↔ Li (Pearl basin) | military supply for the conquest of the south | Qin state | +| Grand Canal (China) | sections from 486 BC; Sui trunk 605–610; Yuan route to Beijing 1280s–1290s | ≈ 1,800 km | summit ≈ 40 m (Shandong) | Yellow River ↔ Yangtze ↔ Beijing | grain tribute to the capital, army supply, unity | Sui state; millions of conscripts (figures vary) | +| Fossa Carolina | 793 | ≈ 3 km dug, never finished | over the Rhine–Danube divide | Rezat (Rhine) ↔ Altmühl (Danube) | link two great river systems | Charlemagne; abandoned (rain, soft ground) — a *failure* | +| Stecknitz Canal | 1391–1398 | ≈ 95 km | summit canal, ≈ a dozen locks | Elbe ↔ Trave → Lübeck (Baltic) | salt trade (Lüneburg) | city of Lübeck; one of Europe's first summit canals | +| Briare Canal | 1604–1642 | ≈ 55 km | summit level, ≈ 40 locks | Loire ↔ Seine | inland trade to Paris | French crown (Sully, Henri IV), then a company | +| Canal du Midi | 1666–1681 | ≈ 240 km | summit ≈ 190 m, ≈ 60–100 locks (counts vary) | Garonne (Atlantic) ↔ Mediterranean | avoid the long, hostile route round Spain via Gibraltar | Riquet with Louis XIV's backing; ≈ 12,000 workers | + +## 2. Early industrial (locks, steam; still fits a magic-assisted age) + +| Canal | Built | Length | Summit / locks | Joins | Purpose | Notes | +|---|---|---|---|---|---|---| +| Bridgewater Canal | 1759–1761 | ≈ 65 km | level | coal mines ↔ Manchester | coal | halved coal prices in Manchester | +| Caledonian Canal | 1803–1822 | ≈ 97 km (≈ ⅓ dug; rest is lochs) | summit ≈ 32 m, 29 locks | North Sea ↔ Atlantic across Scotland | avoid the stormy Pentland Firth; naval use; work for the jobless | built for frigates; ships soon outgrew it, so it was a commercial *disappointment* | +| Göta Canal | 1810–1832 | ≈ 190 km (≈ 87 km dug) | summit ≈ 92 m, 58 locks | Baltic ↔ Kattegat via lakes Vättern and Vänern | avoid the Danish Sound tolls and route | ≈ 58,000 soldiers, ≈ 7 million man-days; obsolete as a trunk route soon after (railways) | +| Erie Canal | 1817–1825 | ≈ 580 km | rise ≈ 170 m, 83 locks | Hudson ↔ Great Lakes | open the interior to the sea | freight cost Buffalo–New York fell by ≈ 90 %; made New York the leading port | +| Welland Canal | 1829 (rebuilt 1932) | ≈ 43 km | ≈ 100 m, 8 locks (today) | Lake Ontario ↔ Lake Erie | bypass Niagara Falls | a lake-to-lake access canal | + +## 3. Ship canals (sea-going ships) + +| Canal | Built | Length | Summit / locks | Saves or opens | Notes | +|---|---|---|---|---|---| +| Suez | 1859–1869 | ≈ 164 km (now ≈ 193) | sea level, no locks | London–Bombay ≈ 7,000 km shorter (≈ 40 %) | company with state backing; forced labour early on; deaths in the tens of thousands (contested) | +| Corinth | 1881–1893 | ≈ 6.4 km | sea level, cut up to ≈ 90 m into rock | ≈ 400 km round the Peloponnese | tried by Periander and Nero (begun AD 67, abandoned); too narrow for big ships, so it never paid | +| Kiel | 1887–1895 | ≈ 98 km | sea level, locks at both ends | ≈ 460 km round Denmark | military first: the German fleet between the Baltic and the North Sea | +| Manchester Ship Canal | 1887–1894 | ≈ 58 km | 5 locks, ≈ 18 m | ocean ships reach an inland city | built to bypass Liverpool's port dues | +| Panama | French 1881–1894 (failed), US 1904–1914 | ≈ 80 km | summit lake ≈ 26 m, 3 lock flights | New York–San Francisco ≈ 13,000 km shorter (vs Cape Horn) | disease (yellow fever, malaria): ≈ 22,000 dead in the French attempt; a sea-level design failed, locks worked | +| Volga–Don | 1948–1952 (attempted 1697 by Peter I) | ≈ 101 km | ≈ 88 m, 13 locks | Caspian ↔ Black Sea river systems | Peter's attempt was abandoned | +| White Sea–Baltic | 1931–1933 | ≈ 227 km (≈ 48 km dug) | 19 locks | White Sea ↔ Baltic | forced (Gulag) labour, ≈ 12,000–25,000 dead; too shallow (≈ 3.5 m) for the ships it was meant for — a *cautionary* case | +| Rhine–Main–Danube | 1960–1992 (Ludwig Canal 1836–1846 before it) | ≈ 171 km | summit ≈ 406 m, 16 locks | North Sea ↔ Black Sea by inland water | finally achieved Charlemagne's link, 1,200 years later | + +## 4. What this suggests for the sim + +**Why a polity digs** — four kinds of gain, often mixed: + +1. *Shorter route* (Suez, Panama, Kiel, Corinth): gain ∝ traffic × distance saved. +2. *Access*, where no water route existed (Erie, Welland, Manchester, Lingqu, the Grand Canal joining river systems): + the largest gains. A whole hinterland or lake joins the sea economy; freight costs fall by up to ≈ 90 %. +3. *Avoiding a hazard or toll* (Xerxes: a storm cape; Caledonian: the Pentland Firth; Göta: the Danish Sound dues; + Manchester: Liverpool's dues). +4. *Military* logistics or fleet movement (Lingqu, Sui Grand Canal, Kiel, Xerxes). States accept a poor commercial + return for this. + +**Who can dig:** centralised states that can command mass labour (Persia, Qin, Sui, Louis XIV's France, Sweden's +army). Later, chartered companies with state backing (Suez, the second Panama). Cities in leagues manage smaller works +(Lübeck's Stecknitz). A polity must hold, or control, both ends. + +**Scale of effort, as a rough guide:** + +- Big pre-industrial projects take **3–15 years** of work by **10,000 to 100,000** people, or millions for the + Grand Canal. +- Cost rises with length, with every lock (summits of 30–200 m are common; ≈ 400 m is the extreme), with rock + (Corinth), and with disease and climate (Panama). +- Deaths can be large and are a story hook (the Pharaohs' canal, Panama, White Sea). + +**Failure is common:** + +- The Fossa Carolina and Nero's Corinth were abandoned; the first Panama attempt went bankrupt. +- The Caledonian and White Sea canals were finished too small to matter. +- Peter I's Volga–Don stalled. + +So the sim should let works fail or come out undersized, not just succeed. + +**Cheaper alternatives:** a *portage road* (Diolkos, used for centuries) gives part of the gain at a fraction of the +cost. It is a natural first stage before a canal. + +**Decay:** canals need upkeep. +- The Pharaohs' canal was reopened and silted up several times (finally closed 767). +- The Grand Canal declined when the state weakened and the Yellow River moved (1855). +- A canal should decay or close when its polity collapses. That fits the decline from High into Low Fantasy. + +**Ship size matters:** a canal that only fits barges (the Midi, the Göta) serves inland trade. One that fits sea +ships (Suez, Kiel, and the Republic's neck cut onto a deep lake) joins two seas for the same fleet. The site scan +should record depth and level on both sides so the sim can tell which kind of canal a site allows. + +## 5. Mapping to the model sketch + +- **gain** = traffic × days saved, plus an *access* bonus when no water route exists, plus a hazard/toll term, plus + a military weight for strong expansionist polities. +- **cost** = length × (1 + climb / k_lock) × terrain (rock, marsh, disease) — labour-years. Labour-years / available + labour gives the build time, and the polity must stay stable that long, or the work fails. +- **options:** a portage road (cheap, partial gain) → a barge canal → a ship canal. +- **upkeep:** a yearly cost; it decays without a strong owner. + +Placeholder numbers only; tune against the cases above (e.g. the Erie access gain, the Suez distance gain, the +Göta labour-years) when sub-project 2 is built. diff --git a/history.py b/history.py new file mode 100755 index 0000000..ac2a969 --- /dev/null +++ b/history.py @@ -0,0 +1,6 @@ +#!/usr/bin/env python3 +"""worldhistory CLI: run / preview / compare / check-config (see README).""" +from worldhistory.cli import main + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/licenses.md b/licenses.md new file mode 100644 index 0000000..439f998 --- /dev/null +++ b/licenses.md @@ -0,0 +1,13 @@ +# Third-party licences + +Python packages are installed by the user from `requirements.txt` (not shipped in this repo). +Licences as published by each project, read 2026-09-30; update with any dependency change. + +| Package | Requirement | Licence | +|---|---|---| +| numpy | >=2.3 | BSD-3-Clause | +| scipy | >=1.16 | BSD-3-Clause | +| pillow | >=11 | MIT-CMU (HPND) | +| h3 (h3-py) | >=4.5,<5 | Apache-2.0 | +| numba (optional: compiled fast paths) | any | BSD-2-Clause | +| llvmlite (with numba) | any | BSD-2-Clause; bundles LLVM: Apache-2.0 with LLVM exception | diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..6f1e3e1 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,4 @@ +numpy>=2.3 +scipy>=1.16 +pillow>=11 +h3>=4.5,<5 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) diff --git a/worldhistory/__init__.py b/worldhistory/__init__.py new file mode 100644 index 0000000..6bd55b8 --- /dev/null +++ b/worldhistory/__init__.py @@ -0,0 +1 @@ +"""worldhistory: population, migration and adaptation trials over thousands of years on a worldgen world.""" diff --git a/worldhistory/adaptation.py b/worldhistory/adaptation.py new file mode 100644 index 0000000..8a2fed0 --- /dev/null +++ b/worldhistory/adaptation.py @@ -0,0 +1,53 @@ +"""Adaptation (spec §3.2): the optimum moves toward the conditions people live in; the width stays fixed, so +adapting to the new costs the old.""" +import numpy as np + +from .config import CONDITIONS + + +def adapt(st, env, races, q_eff, years, mult=None): + """The optimum drifts toward the conditions lived in (≤ rate·years × turnover). Past the race's adaptable range + [lo, hi] outward drift is braked by exp(−½(e/s)²) (e: distance past the edge, s: width on that side): a fuzzy + limit, not a line. Conditions with no meaning in a cell (rain at sea) do not move.""" + for r, race in enumerate(races): + oc = np.flatnonzero(st.P[r] > 0) # only occupied cells drift: compute there + if not len(oc): + continue + turnover = 0.3 + 0.7 * q_eff[r][oc] + m = 1.0 if mult is None else mult[r] + m = m[oc] if np.ndim(m) else m + for ci, c in enumerate(CONDITIONS): + tol = race.tolerance.get(c) + if tol is None or tol.rate <= 0: + continue + lo, hi = env.lo[c][oc], env.hi[c][oc] + live = np.isfinite(lo) & np.isfinite(hi) + target = np.where(np.abs(hi - tol.optimum) >= np.abs(lo - tol.optimum), hi, lo) + O = st.O[r, ci, oc] + lim = tol.rate * years * m + step = np.clip(np.where(live, target, O) - O, -lim, lim) * turnover + past_hi, past_lo = np.maximum(O - tol.hi, 0), np.maximum(tol.lo - O, 0) + brake = np.where(step > 0, np.exp(-0.5 * (past_hi / tol.side(True)) ** 2), + np.exp(-0.5 * (past_lo / tol.side(False)) ** 2)) + st.O[r, ci, oc] = O + step * brake + + +def prospective(st, world): + """Empty cells take the population-weighted optimum of their occupied neighbours: whether a cell is habitable is + judged for the people who would move in, not for a stale value.""" + for r in range(st.P.shape[0]): + pos = st.P[r] > 0 + if not pos.any(): # nobody anywhere: no neighbour mix to take + continue + near = world.nb_local(np.flatnonzero(pos))[1] # cells next to a group (adjacency is symmetric): only + cand = near[st.P[r][near] <= 0] # empty ones among them can have a neighbour sum > 0 + if not len(cand): + continue + sub, cols = world.nb_local(cand) + Pn = np.where(pos[cols], st.P[r][cols], 0.0) + wsum = world.nb_apply(sub, Pn) + keep = wsum > 0 # any group counts, however thin (spec 2026-09-30 §1.3) + if not keep.any(): + continue + empty = cand[keep] + st.O[r][:, empty] = world.nb_apply(sub, Pn * st.O[r][:, cols])[:, keep] / np.maximum(wsum[keep], 1e-12) diff --git a/worldhistory/cli.py b/worldhistory/cli.py new file mode 100644 index 0000000..e48ed2a --- /dev/null +++ b/worldhistory/cli.py @@ -0,0 +1,69 @@ +"""history.py run | preview | compare | check-config.""" +import argparse +import json +import sys +from pathlib import Path + +from .config import ConfigError, load_config + + +def main(argv=None): + ap = argparse.ArgumentParser(prog="history.py", description="Population and migration trials on a worldgen world.") + sub = ap.add_subparsers(dest="cmd", required=True) + run = sub.add_parser("run", help="run a trial") + run.add_argument("--world", required=True, help="worldgen build dir, e.g. out/r4") + run.add_argument("--config", required=True, help="dir with history.toml and races/") + run.add_argument("--trial", required=True, help="trial name (output: <world>/history/<trial>)") + run.add_argument("--seed", type=int) + run.add_argument("--years", type=int) + run.add_argument("--out", help="output dir instead of <world>/history/<trial>") + run.add_argument("--no-preview", action="store_true") + pv = sub.add_parser("preview", help="render PNGs for a finished trial") + pv.add_argument("trial_dir") + cp = sub.add_parser("compare", help="compare two trials (seed agreement)") + cp.add_argument("a") + cp.add_argument("b") + cp.add_argument("--out") + ck = sub.add_parser("check-config", help="validate a config dir without running") + ck.add_argument("config") + ck.add_argument("--world", help="also check fields and eras against this build dir") + args = ap.parse_args(argv) + try: + if args.cmd == "check-config": + fields = None + if args.world: + import numpy as np + fields = set(np.load(Path(args.world) / "cells.npz").files) + h, races = load_config(args.config, fields) + if args.world: + for e in h.eras: + if not (Path(args.world) / "eras" / e / "cells.npz").exists(): + raise ConfigError(f"era {e!r}: {args.world}/eras/{e}/cells.npz missing") + print(f"ok: {len(races)} races ({', '.join(r.id for r in races)}), {len(h.events)} events, " + f"{len(h.hazards)} hazards") + return 0 + if args.cmd == "run": + from .engine import Engine + from .preview import preview + from .world import load_world + h, races = load_config(args.config) + w = load_world(args.world, eras=h.eras, extra_fields=h.world_fields) + out = Path(args.out) if args.out else Path(args.world) / "history" / args.trial + stats = Engine(w, h, races, seed=args.seed).run(out, years=args.years, progress=print) + print(f"done in {stats['timing']['total_s']:.0f} s, peak {stats['timing']['peak_rss_mb']:.0f} MB -> {out}") + if not args.no_preview: + preview(out, world=w) + return 0 + if args.cmd == "preview": + from .preview import preview + for p in preview(args.trial_dir): + print(p) + return 0 + if args.cmd == "compare": + from .preview import compare + print(json.dumps(compare(args.a, args.b, out=args.out), indent=1)) + return 0 + except ConfigError as e: + print(f"config error: {e}", file=sys.stderr) + return 1 + return 2 diff --git a/worldhistory/config.py b/worldhistory/config.py new file mode 100644 index 0000000..0782c56 --- /dev/null +++ b/worldhistory/config.py @@ -0,0 +1,465 @@ +"""Config: history.toml + races/<race>.toml -> History and Race objects, with defaults and validation. + +Unknown keys, predicates, conditions, regions, races and eras are errors (ConfigError names file and key).""" +from __future__ import annotations + +import copy +import json +import re +import tomllib +from dataclasses import asdict, dataclass, field +from pathlib import Path + +from . import predicates +from .regions import KINDS as REGION_KINDS + +CONDITIONS = ("depth", "air_pressure", "gravity", "o2", "temperature", "rain") +DOMAINS = ("farming", "crafts", "travel", "magic", "war", "land") +EVENT_KINDS = ("seed", "cull", "die_off", "era_switch") +HAZARD_KINDS = ("static", "roaming") + + +class ConfigError(ValueError): + pass + + +@dataclass +class Tolerance: + optimum: float + lo: float + hi: float + width: float | None = None # sets both sides + width_lo: float | None = None # how fast fitness falls below the optimum + width_hi: float | None = None # … and above it + rate: float = 0.0 # condition units per year the optimum can move + comfort: list = field(default_factory=list) # [[optimum, multiplier], ...] + + def __post_init__(self): + self.width_lo = self.width if self.width_lo is None else self.width_lo + self.width_hi = self.width if self.width_hi is None else self.width_hi + if self.width_lo is None or self.width_hi is None: + raise TypeError("give width, or both width_lo and width_hi") + + def side(self, up): + """Width on the side above (up=True) or below the optimum.""" + return self.width_hi if up else self.width_lo + + +@dataclass +class Race: + id: str + name: str + family: str + colour: list + density: float # people/km² at habitat quality 1, before tech + growth: float # per step + allee: float + founder: float + mobility: float + roam: float # share of every group that moves on each step (nomads) + jump: float # long-jump rate multiplier (e.g. swimmers crossing chasms) + range: int # own crowding over this many neighbour rings; -1 = the whole world + group_max: float # groups above this split: half leaves as a new group (0 = never) + jump_path: str # any (through the air) | sea (swimmers: the path stays in water) + current_bias: float # sea jumps/drift favour going with the current (0 = ignore) + clusters: int + heads: float + seed_realm: str # "" = anywhere habitable | land | sea (gifting placement) + realm: str # land | sea | both + terms: list + gates: list + hook: str + qmax: float + travel: dict + tolerance: dict # condition -> Tolerance + temperament: dict + overlap: dict # other race id -> alpha override + conflict: dict + tech: dict # domain -> priority + reach_cap: float + veins: float + emerge: dict | None + + +@dataclass +class History: + years: int + step: int + snapshot_every: int + seed: int + fringe: float + f_min: float + eras: list + world_fields: list + regions: dict + events: list + hazards: list + nudges: list + tech: dict + migration: dict + competition: dict + stats_regions: list + + +RACE_SECTIONS = { + "demography": {"density": None, "growth": None, "allee": 0.0, "founder": None, "mobility": 0.05, "roam": 0.0, + "jump": 1.0, "jump_path": "any", "current_bias": 0.0, "range": 0, "group_max": 0.0}, + "seed": {"clusters": 0, "heads": 0.0, "realm": ""}, + "habitat": {"realm": "land", "terms": None, "gates": [], "hook": "", "qmax": 1.0}, + "travel": {"river": -1.0, "coast": -1.0, "mountain": 0.0, "desert": 0.0, "sea": 0.0}, + "temperament": {"pace": 0.5, "structure": 0.5, "warlike": 0.5, "xenophobia": 0.5, "tolerated_share": 0.2}, + "conflict": {"internal": 0.0, "border": 0.0, "aggression": 0.5, "defend": 1.0, "yield": 0.0, "power": 1.0, + "dread": 0.0, "curse": 0.0}, + "tech": {**{d: 1.0 for d in DOMAINS}, "reach_cap": 0.2}, + "special": {"veins": 0.0}, +} +RACE_TOP = {"id", "name", "family", "colour", "tolerance", "overlap", "emerge", *RACE_SECTIONS} +INHERIT = ("tolerance", "temperament", "conflict", "tech", "travel") # sub-species copy these from the parent +TOL_KEYS = {"optimum", "width", "width_lo", "width_hi", "lo", "hi", "rate", "comfort"} +EMERGE_KEYS = {"parent", "mode", "birth_min", "birth_sure", "birth_rate", "mix_min", "dominance", "birth_shift", + "displace", "return_min", "return_sure", "convert", "region", "condition", "settlement_density", + "isolated", "spread", "trigger", "ritual_rate", "ritual_cost", "ritual_converts", "after", + "return_rate", "crisis_drop", "crisis_min", "crisis_years", "victims", "backlash_victims", + "backlash_loss", "backlash_calm", "raid_rate", "raid_size", "raid_min", "raid_km", "ritual_power", + "backlash_years", "crisis_km"} +RITUAL_ONLY = ("trigger", "ritual_rate", "ritual_cost", "ritual_converts", "crisis_drop", "crisis_min", "crisis_years", + "victims", "backlash_victims", "backlash_loss", "backlash_calm", "raid_rate", "raid_size", "raid_min", + "raid_km", "ritual_power", "backlash_years", "crisis_km") +EMERGE_REMOVED = {"class", "years"} +EMERGE_DEFAULTS = {"mode": "adapt", "birth_min": 0.4, "birth_sure": 0.8, "birth_rate": 0.2, "mix_min": 0.3, + "dominance": 1.0, "birth_shift": 1.0, "displace": 0.9, "return_min": None, "return_sure": None, + "convert": 1.0, "region": "", "condition": {}, "settlement_density": 0.0, "isolated": False, + "spread": "contact", "trigger": 0.01, "ritual_rate": 0.02, "ritual_cost": 300.0, + "ritual_converts": 30.0, "after": 0, + "return_rate": None, "crisis_drop": 0.0, "crisis_min": 1000.0, "crisis_years": 200.0, + "victims": "parent", "backlash_victims": 0.0, "backlash_loss": 0.7, "backlash_calm": 0.2, + "raid_rate": 0.0, "raid_size": 100.0, "raid_min": 1000.0, "raid_km": 300.0, "ritual_power": 0.5, + "backlash_years": 0, "crisis_km": 0.0} + +INF = float("inf") +EMERGE_RANGES = [ # key, lo, hi, lo excluded + ("displace", 0, 1, False), ("dominance", 0, 1, False), ("birth_shift", 0, 1, False), ("convert", 0, 1, True), + ("mix_min", 0, 1, False), ("birth_rate", 0, 1, False), ("return_rate", 0, 1, False), ("after", 0, INF, False), + ("trigger", 0, 1, False), ("ritual_rate", 0, INF, False), ("ritual_cost", 0, INF, False), + ("ritual_converts", 0, INF, False), ("ritual_power", 0, INF, True), ("crisis_min", 0, INF, False), + ("crisis_years", 0, INF, True), ("backlash_victims", 0, INF, False), ("backlash_years", 0, INF, False), + ("raid_rate", 0, INF, False), ("raid_size", 0, INF, True), ("raid_min", 0, INF, False), ("raid_km", 0, INF, False), ("crisis_km", 0, INF, False), + ("settlement_density", 0, INF, False)] + +HISTORY_SECTIONS = { + "run": {"years": 7500, "step": 10, "snapshot_every": 50, "seed": 1, "fringe": 0.2, "f_min": 0.1}, + "world": {"eras": [], "fields": []}, + "tech": {"rate": 0.004, "n_half": 5000.0, "loss_below": 200.0, "loss_rate": 0.02, "diffuse": 0.1, "passes": 1, + "farming_gain": 11.0, "crafts_gain": 1.0, "crafts_shock": 0.5, "crafts_reach": 0.3, "travel_cost": 1.0, + "travel_jump": 1.0, "magic_share": 0.5, "magic_reach": 0.5, "land_rate": 0.05, "land_cap": 0.5, + "land_decay": 0.02}, + "migration": {"pf": 0.25, "temp": 4.0, "hostile": 1.5, "jump_rate": 0.015, "jump_xm_km": 80.0, + "jump_alpha": 1.5, "jump_cap_km": 1500.0}, + "competition": {"intolerance": 2.0, "curse_decay": 0.8}, + "stats": {"regions": []}, +} +HISTORY_TOP = {"regions", "event", "hazard", "nudge", *HISTORY_SECTIONS} +EVENT_KEYS = {"year", "kind", "region", "share", "by_race", "noise", "era", "race", "clusters", "heads"} +HAZARD_KEYS = { + "static": {"kind", "name", "region", "mortality", "years", "decay", "push"}, + "roaming": {"kind", "name", "start", "home_region", "home_near", "local", "wanderers", "radius_km", + "local_mortality", "raid_chance", "raid_burst", "raid_death", "density_ref", "wander_step_km", + "wander_stop", "wander_burst", "growth", "cap", "clutch", "infight"}, +} +ROAMING_DEFAULTS = {"home_near": "", "local": 20, "wanderers": 2, "radius_km": 150.0, "local_mortality": 0.01, + "raid_chance": 0.02, "raid_burst": 0.3, "raid_death": 0.05, "density_ref": 5.0, + "wander_step_km": 800.0, "wander_stop": 0.3, "wander_burst": 0.2, "growth": 0.01, "cap": 0, + "clutch": [2, 4], "infight": 0.002} +NUDGE_KEYS = {"race", "region", "years", "growth", "capacity", "expansion", "adapt"} + + +def _keys(d, allowed, where): + if not isinstance(d, dict): + raise ConfigError(f"{where}: expected a table") + bad = sorted(set(d) - set(allowed)) + if bad: + raise ConfigError(f"{where}: unknown key(s) {bad}") + + +def _section(d, name, defaults, where): + sec = d.get(name, {}) + _keys(sec, defaults, f"{where} [{name}]") + return {**defaults, **copy.deepcopy(sec)} + + +def _pred(expr, where, field_names=None): + try: + predicates.check(expr, field_names) + except ValueError as e: + raise ConfigError(f"{where}: {e}") from None + + +def parse_race(d, where): + _keys(d, RACE_TOP, where) + if "id" not in d: + raise ConfigError(f"{where}: missing id") + s = {k: _section(d, k, v, where) for k, v in RACE_SECTIONS.items()} + for sec, key in (("demography", "density"), ("demography", "growth"), ("habitat", "terms")): + if s[sec][key] is None: + raise ConfigError(f"{where} [{sec}]: missing {key}") + hab = s["habitat"] + if hab["realm"] not in ("land", "sea", "both"): + raise ConfigError(f"{where} [habitat]: realm must be land, sea or both") + if s["demography"]["jump_path"] not in ("any", "sea"): + raise ConfigError(f"{where} [demography]: jump_path must be any or sea") + if s["seed"]["realm"] not in ("", "land", "sea"): + raise ConfigError(f"{where} [seed]: realm must be land or sea") + for t in hab["terms"]: + _keys(t, {"p", "w"}, f"{where} [habitat] term") + _pred(t["p"], f"{where} [habitat] term") + for g in hab["gates"]: + _keys(g, {"require", "penalty", "factor", "realm"}, f"{where} [habitat] gate") + if ("require" in g) == ("penalty" in g): + raise ConfigError(f"{where} [habitat] gate: give exactly one of require / penalty") + _pred(g.get("require", g.get("penalty")), f"{where} [habitat] gate") + if "penalty" in g and "factor" not in g: + raise ConfigError(f"{where} [habitat] gate: penalty needs factor") + tol = {} + for c, t in d.get("tolerance", {}).items(): + if c not in CONDITIONS: + raise ConfigError(f"{where} [tolerance.{c}]: unknown condition; known: {list(CONDITIONS)}") + _keys(t, TOL_KEYS, f"{where} [tolerance.{c}]") + try: + tol[c] = Tolerance(**t) + except TypeError as e: + raise ConfigError(f"{where} [tolerance.{c}]: {e}") from None + em = d.get("emerge") + if em is not None: + old = sorted(EMERGE_REMOVED & set(em)) + if old: + raise ConfigError(f"{where} [emerge]: {old} removed (gradual adaptation, 2026-09-30); use birth_min / " + f"birth_sure / birth_rate and the sub-species' own [tolerance.*] curves") + _keys(em, EMERGE_KEYS, f"{where} [emerge]") + if "parent" not in em: + raise ConfigError(f"{where} [emerge]: missing parent") + mode = em.get("mode", "adapt") + if mode not in ("adapt", "ritual"): + raise ConfigError(f"{where} [emerge]: mode must be adapt or ritual") + if mode == "adapt" and any(k in em for k in RITUAL_ONLY): + raise ConfigError(f"{where} [emerge]: {', '.join(RITUAL_ONLY)} only with mode = \"ritual\"") + em = {**copy.deepcopy(EMERGE_DEFAULTS), **em} + em["return_min"] = em["birth_min"] if em["return_min"] is None else em["return_min"] + em["return_sure"] = em["birth_sure"] if em["return_sure"] is None else em["return_sure"] + em["return_rate"] = em["birth_rate"] if em["return_rate"] is None else em["return_rate"] + if em["victims"] not in ("parent", "all"): + raise ConfigError(f"{where} [emerge]: victims must be parent or all") + for key, lo, hi, lo_open in EMERGE_RANGES: + x = em[key] + if not isinstance(x, (int, float)) or x < lo or x > hi or (lo_open and x == lo): + raise ConfigError(f"{where} [emerge]: {key} = {x!r} out of range " + f"{'(' if lo_open else '['}{lo}, {hi}]") + if not 0 <= em["crisis_drop"] < 1 or not 0 <= em["backlash_loss"] <= 1 or em["backlash_calm"] < 0: + raise ConfigError(f"{where} [emerge]: need 0 ≤ crisis_drop < 1, 0 ≤ backlash_loss ≤ 1, backlash_calm ≥ 0") + if not 0 <= em["birth_min"] <= em["birth_sure"] <= 1 or not 0 <= em["return_min"] <= em["return_sure"] <= 1: + raise ConfigError(f"{where} [emerge]: need 0 ≤ birth_min ≤ birth_sure ≤ 1 (and the same for return_*)") + if em["spread"] not in ("contact", "region"): + raise ConfigError(f"{where} [emerge]: spread must be contact or region") + demo = s["demography"] + race = Race( + id=d["id"], name=d.get("name", d["id"]), family=d.get("family", d["id"]), colour=list(d.get("colour", [200, 200, 200])), + density=float(demo["density"]), growth=float(demo["growth"]), allee=float(demo["allee"]), + founder=float(demo["founder"] if demo["founder"] is not None else max(demo["allee"], 1.0)), + mobility=float(demo["mobility"]), roam=float(demo["roam"]), jump=float(demo["jump"]), jump_path=demo["jump_path"], current_bias=float(demo["current_bias"]), + range=int(demo["range"]), group_max=float(demo["group_max"]), + clusters=0 if em else int(s["seed"]["clusters"]), heads=float(s["seed"]["heads"]), + seed_realm=s["seed"]["realm"], + realm=hab["realm"], terms=hab["terms"], gates=hab["gates"], hook=hab["hook"], qmax=float(hab["qmax"]), + travel=s["travel"], tolerance=tol, temperament=s["temperament"], overlap=dict(d.get("overlap", {})), + conflict=s["conflict"], tech={k: v for k, v in s["tech"].items() if k != "reach_cap"}, + reach_cap=float(s["tech"]["reach_cap"]), veins=float(s["special"]["veins"]), emerge=em) + race._given = {k for k in INHERIT if k in d} # which sections the file set (for inheritance) + return race + + +def parse_history(d, where): + if "exposure" in d: + raise ConfigError(f"{where}: [exposure] removed (gradual adaptation, 2026-09-30); sub-species now come from " + f"tolerance curves — see [emerge] birth_min / birth_sure") + _keys(d, HISTORY_TOP, where) + s = {k: _section(d, k, v, where) for k, v in HISTORY_SECTIONS.items()} + regions = d.get("regions", {}) + for name, spec in regions.items(): + _keys(spec, {"kind", "lat", "lon", "era", "fields", "eps", "ids", "p", "land", "exclude"}, + f"{where} [regions.{name}]") + if spec.get("kind", "all") not in REGION_KINDS: + raise ConfigError(f"{where} [regions.{name}]: unknown kind {spec.get('kind')!r}") + events = [] + for i, ev in enumerate(d.get("event", [])): + _keys(ev, EVENT_KEYS, f"{where} [[event]] #{i}") + if ev.get("kind") not in EVENT_KINDS: + raise ConfigError(f"{where} [[event]] #{i}: kind must be one of {list(EVENT_KINDS)}") + events.append({"region": "all", "share": 0.0, "by_race": {}, "noise": 0.0, **ev}) + hazards = [] + for i, hz in enumerate(d.get("hazard", [])): + kind = hz.get("kind") + if kind not in HAZARD_KINDS: + raise ConfigError(f"{where} [[hazard]] #{i}: kind must be one of {list(HAZARD_KINDS)}") + _keys(hz, HAZARD_KEYS[kind], f"{where} [[hazard]] #{i}") + base = ({"name": f"hazard{i}", "region": "all", "mortality": 0.0, "years": [0, 10 ** 9], "decay": 1.0, + "push": 0.0} + if kind == "static" else {"name": f"hazard{i}", "start": 0, "home_region": "all", **ROAMING_DEFAULTS}) + hazards.append({**base, **hz}) + nudges = [] + for i, nd in enumerate(d.get("nudge", [])): + _keys(nd, NUDGE_KEYS, f"{where} [[nudge]] #{i}") + nudges.append({"region": "all", "years": [0, 10 ** 9], "growth": 1.0, "capacity": 1.0, "expansion": 1.0, + "adapt": 1.0, **nd}) + run = s["run"] + return History(years=int(run["years"]), step=int(run["step"]), snapshot_every=int(run["snapshot_every"]), + seed=int(run["seed"]), fringe=float(run["fringe"]), f_min=float(run["f_min"]), + eras=list(s["world"]["eras"]), world_fields=list(s["world"]["fields"]), regions=regions, + events=events, hazards=hazards, nudges=nudges, tech=s["tech"], + migration=s["migration"], competition=s["competition"], stats_regions=list(s["stats"]["regions"])) + + +def distinct_axes(r_from, r_to): + """Conditions whose optima differ by more than ¼ of r_from's width on that side (the axes progress counts).""" + out = [] + for c in CONDITIONS: + a, b = r_from.tolerance.get(c), r_to.tolerance.get(c) + if a is not None and b is not None and abs(b.optimum - a.optimum) > 0.25 * a.side(b.optimum > a.optimum): + out.append(c) + return out + + +def link(history, races, field_names=None): + """Cross-checks and sub-species inheritance. Returns races with every parent before its sub-species.""" + ids = [r.id for r in races] + if len(set(ids)) != len(ids): + raise ConfigError(f"duplicate race ids in {ids}") + by_id = {r.id: r for r in races} + region_names = set(history.regions) | {"all"} + + def region(name, where): + if name and name not in region_names: + raise ConfigError(f"{where}: unknown region {name!r}; known: {sorted(region_names)}") + + for name, spec in history.regions.items(): + if spec.get("exclude"): + region(spec["exclude"], f"[regions.{name}] exclude") + if spec.get("kind") == "changed" and spec.get("era") not in history.eras: + raise ConfigError(f"[regions.{name}]: era {spec.get('era')!r} not in [world] eras") + for i, ev in enumerate(history.events): + w = f"[[event]] #{i} ({ev['kind']}, year {ev['year']})" + region(ev["region"], w) + for rid in ev["by_race"]: + if rid not in by_id: + raise ConfigError(f"{w}: unknown race {rid!r} in by_race") + if ev.get("race") and ev["race"] not in by_id: + raise ConfigError(f"{w}: unknown race {ev['race']!r}") + if ev["kind"] == "era_switch" and ev.get("era") not in history.eras: + raise ConfigError(f"{w}: era {ev.get('era')!r} not in [world] eras") + for i, hz in enumerate(history.hazards): + w = f"[[hazard]] #{i} ({hz['name']})" + region(hz.get("region", ""), w) + region(hz.get("home_region", ""), w) + if hz.get("home_near") and hz["home_near"] not in by_id: + raise ConfigError(f"{w}: unknown race {hz['home_near']!r} in home_near") + for i, nd in enumerate(history.nudges): + region(nd["region"], f"[[nudge]] #{i}") + if nd.get("race") not in by_id: + raise ConfigError(f"[[nudge]] #{i}: unknown race {nd.get('race')!r}") + for name in history.stats_regions: + region(name, "[stats] regions") + if field_names is not None: + for r in races: + for t in r.terms: + _pred(t["p"], f"race {r.id}: term", field_names) + for r in races: + for other in r.overlap: + if other not in by_id: + raise ConfigError(f"race {r.id} [overlap]: unknown race {other!r}") + if r.emerge is None: + continue + p = by_id.get(r.emerge["parent"]) + if p is None: + raise ConfigError(f"race {r.id} [emerge]: unknown parent {r.emerge['parent']!r}") + region(r.emerge["region"], f"race {r.id} [emerge]") + r.family = p.family + for sec in INHERIT: + if sec == "tolerance": + merged = copy.deepcopy(p.tolerance) + merged.update(r.tolerance) # own curves win, per condition + r.tolerance = merged + elif sec not in getattr(r, "_given", set()): + setattr(r, sec, copy.deepcopy(getattr(p, sec))) + em = r.emerge + if em["mode"] == "adapt" and em["birth_min"] > 0 and not distinct_axes(p, r): + raise ConfigError(f"race {r.id}: no tolerance differs enough from {p.id}'s (> ¼ width), so progress " + f"stays 0 < birth_min and it would never be born") + parents = [r for r in races if r.emerge is None] + subs = [r for r in races if r.emerge is not None] + return parents + subs + + +def load_config(directory, field_names=None): + d = Path(directory) + try: + history = parse_history(tomllib.loads((d / "history.toml").read_text()), str(d / "history.toml")) + races = [parse_race(tomllib.loads(p.read_text()), str(p)) for p in sorted((d / "races").glob("*.toml"))] + except tomllib.TOMLDecodeError as e: + raise ConfigError(f"{d}: {e}") from None + if not races: + raise ConfigError(f"{d / 'races'}: no race files") + return history, link(history, races, field_names) + + +def _plain(x): + if isinstance(x, dict): + return {str(k): _plain(v) for k, v in x.items() if v is not None} + if isinstance(x, (list, tuple)): + return [_plain(v) for v in x] + if hasattr(x, "__dataclass_fields__"): + return _plain(asdict(x)) + return x + + +def resolved(history, races, **extra): + return _plain({**extra, "history": asdict(history), "race": [asdict(r) for r in races]}) + + +_BARE = re.compile(r"^[A-Za-z0-9_-]+$") + + +def _key(k): + return k if _BARE.match(k) else json.dumps(k) + + +def _val(v): + if isinstance(v, bool): + return "true" if v else "false" + if isinstance(v, (int, float)): + return repr(v) + if isinstance(v, str): + return json.dumps(v) + if isinstance(v, list): + return "[" + ", ".join(_val(x) for x in v) + "]" + if isinstance(v, dict): + return "{ " + ", ".join(f"{_key(k)} = {_val(x)}" for k, x in v.items()) + " }" + raise TypeError(f"cannot write {type(v).__name__} to TOML") + + +def dumps_toml(d, prefix=""): + """Minimal TOML writer for the run record: scalars, lists, tables, arrays of tables.""" + lines, tables = [], [] + for k, v in d.items(): + if isinstance(v, dict) and v: + tables.append((k, v)) + elif isinstance(v, list) and v and all(isinstance(x, dict) for x in v): + tables.append((k, v)) + else: + lines.append(f"{_key(k)} = {_val(v)}") + out = "\n".join(lines) + for k, v in tables: + name = f"{prefix}.{_key(k)}" if prefix else _key(k) + if isinstance(v, dict): + out += f"\n\n[{name}]\n" + dumps_toml(v, name) + else: + for item in v: + out += f"\n\n[[{name}]]\n" + dumps_toml(item, name) + return out.strip() + "\n" diff --git a/worldhistory/conflict.py b/worldhistory/conflict.py new file mode 100644 index 0000000..fe7fa04 --- /dev/null +++ b/worldhistory/conflict.py @@ -0,0 +1,56 @@ +"""Background war losses until polities exist (spec §5.2): internal wars, border fights, defended homelands, +yielding (moving on instead of fighting), dread (feared peoples are left alone and avoided) and curses (those who +kill a cursing people's folk die of it for years after).""" +import numpy as np + + +def conflict(P, q_eff, crowd, races, core_q=0.6, blame=None): + R = len(races) + u = np.flatnonzero((P != 0).any(0)) # cells with anyone: elsewhere every term is 0 (loss and press 0) + Pu = P[:, u] + tot = Pu[0].copy() # race by race, as P.sum(0) adds the rows (Pu is column-major: + for i in range(1, R): # its own sum(0) would add pairwise, other roundings) + tot += Pu[i] + share = np.where(tot > 0, Pu / np.maximum(tot, 1e-12), 0.0) + fam = [r.family for r in races] + C = [r.conflict for r in races] + loss = np.zeros_like(Pu) + press = np.zeros_like(Pu) + for i in range(R): + loss[i] += C[i]["internal"] * np.minimum(crowd[i, u], 2.0) * Pu[i] + for j in range(R): # per attacker j, all its targets at once: each (i, cell) still + x = np.flatnonzero(share[j] != 0) # gets the same terms, in j order (where j has no share they + I = [i for i in range(R) if fam[i] != fam[j]] # all add 0) + if not len(x) or not I: + continue + cj, ci = C[j], [C[i] for i in I] + sj = share[j, x] + k_att = np.array([cj["aggression"] * (1 - c["dread"]) * cj["power"] / c["power"] for c in ci]) + k_atk = np.array([c["aggression"] * cj["defend"] for c in ci]) + p_i = np.array([c["power"] for c in ci]) + attacked = k_att[:, None] * sj + attacking = k_atk[:, None] * (q_eff[j, u[x]] >= core_q) * cj["power"] / p_i[:, None] * sj + ix = np.ix_(I, x) + press[ix] += attacked + cj["dread"] * sj + loss[ix] += np.array([c["border"] for c in ci])[:, None] * (attacked + attacking) * Pu[ix] + if blame is not None: + for n, c in enumerate(ci): + if c["curse"] > 0: + blame[j, u[x]] += c["curse"] * c["border"] * attacked[n] + out, press_all = np.zeros_like(P), np.zeros_like(P) + out[:, u] = np.minimum(loss, 0.9 * Pu) + press_all[:, u] = press + return out, press_all + + +def apply_conflict(st, loss, races): + for r, race in enumerate(races): + y = race.conflict["yield"] + st.P[r] = np.maximum(st.P[r] - (1 - y) * loss[r], 0.0) + st.push[r] += y * loss[r] + + +def apply_curse(st, blame, decay): + """Curses on perpetrators: a per-step death share that travels with the cursed and fades by decay per step.""" + st.C[:, 0] = st.C[:, 0] * decay + blame + st.P *= 1 - np.clip(st.C[:, 0], 0, 0.5) diff --git a/worldhistory/demography.py b/worldhistory/demography.py new file mode 100644 index 0000000..8eda86f --- /dev/null +++ b/worldhistory/demography.py @@ -0,0 +1,106 @@ +"""Capacity, growth and competition. Logistic growth (Verhulst 1838) with an Allee effect (Allee 1931; +Courchamp et al. 1999); races share capacity by Lotka–Volterra competition (Lotka 1925; Volterra 1926).""" +import numpy as np + + +def overlap_matrix(suits, races): + """alpha[i, j]: how much race j's people count against race i's capacity. Habitat similarity (cosine over + cells) × temperament friction; asymmetric (i's xenophobia, j's warlikeness). TOML overrides win.""" + S = np.asarray(suits, float) + norm = np.linalg.norm(S, axis=1) + 1e-12 + sim = (S @ S.T) / np.outer(norm, norm) + R = len(races) + a = np.zeros((R, R)) + for i, ri in enumerate(races): + for j, rj in enumerate(races): + if i == j: + continue + if rj.id in ri.overlap: + a[i, j] = ri.overlap[rj.id] + elif ri.family == rj.family: + a[i, j] = sim[i, j] + else: + ti, tj = ri.temperament, rj.temperament + friction = (0.5 + 0.5 * ti["xenophobia"] + 0.5 * abs(ti["pace"] - tj["pace"]) + + 0.5 * tj["warlike"] * abs(ti["structure"] - tj["structure"])) + a[i, j] = np.clip(sim[i, j] * friction, 0, 2) + return a + + +def capacity(world, race, q_eff, fit, comf, kgain, shock, mult): + """Low gear: the tech gain scales with q², so poor land gains little.""" + return race.density * world.area * q_eff * (1 + kgain * q_eff ** 2) * shock * fit * comf * mult + + +def weather(world, rng): + """One regionally correlated standard-normal draw per cell for this step (shared by all races).""" + e = rng.normal(size=world.n) + return (e + world.nb_sum(e)) / np.sqrt(1 + (world.nb >= 0).sum(1)) + + +def harvest_shocks(world, q, scale, e): + sig = 0.3 * q * scale + return np.exp(sig * e - 0.5 * sig ** 2) + + +def _spread(world, x, passes): + for _ in range(passes): + x = x + world.nb_sum(x) + return x + + +def crowding(P, K, alpha, races, intolerance, world=None): + R = len(races) + tot = P.sum(0) + fam = [r.family for r in races] + out = np.zeros_like(P) + for i, ri in enumerate(races): + if ri.range == 0: # local crowding: only where the race can live (K > 0) + x = np.flatnonzero(K[i] > 0) + own = sum(P[j, x] for j in range(R) if fam[j] == fam[i]) + t = tot[x] + other_share = np.where(t > 0, (t - own) / np.maximum(t, 1e-12), 0.0) + m = np.where(other_share > ri.temperament["tolerated_share"], intolerance, 1.0) + comp = sum(alpha[i, j] * P[j, x] for j in range(R) if j != i) if R > 1 else 0.0 + out[i, x] = (P[i, x] + m * comp) / np.maximum(K[i, x], 1e-12) + continue + own = sum(P[j] for j in range(R) if fam[j] == fam[i]) + other_share = np.where(tot > 0, (tot - own) / np.maximum(tot, 1e-12), 0.0) + m = np.where(other_share > ri.temperament["tolerated_share"], intolerance, 1.0) + comp = sum(alpha[i, j] * P[j] for j in range(R) if j != i) if R > 1 else 0.0 + num, cap = P[i] + m * comp, K[i] + if ri.range < 0: # whole people vs whole capacity (few, scattered clans) + num, cap = np.full_like(num, num.sum()), np.full_like(cap, cap.sum()) + elif ri.range > 0 and world is not None: # wide-ranging: crowding over the neighbourhood + num, cap = _spread(world, num, ri.range), _spread(world, cap, ri.range) + out[i] = np.where(K[i] > 0, num / np.maximum(cap, 1e-12), 0.0) + return out + + +def allee(P, a): + return np.clip((P - a) / (P + a + 1e-12), -0.3, 1.0) + + +def grow(P, K, crowd, q_eff, races, gmult): + """Returns this step's gross births per race and cell (r·P where the group can live): the generational + turnover through which a dominant trait spreads, even in a group at capacity.""" + births = np.zeros_like(P) + for i, race in enumerate(races): + r = race.growth * (0.3 + 0.7 * q_eff[i]) * gmult[i] + g = r * P[i] * np.maximum(1 - crowd[i], -1.0) * allee(P[i], race.allee) + g = np.where(K[i] > 0, g, -0.5 * P[i]) + births[i] = np.where(K[i] > 0, r * P[i], 0.0) + P[i] = np.maximum(P[i] + g, 0.0) + return births + + +def veins(P_r, K_r, rate, a): + """New people wake from the land (not births), only beside a viable community (P >= a).""" + return P_r + rate * np.maximum(K_r - P_r, 0) * (P_r >= a) + + +def stochastic_round(P, rng, below=5.0): + small = (P > 0) & (P < below) + v = P[small] + fl = np.floor(v) + P[small] = fl + (rng.random(v.shape) < (v - fl)) diff --git a/worldhistory/engine.py b/worldhistory/engine.py new file mode 100644 index 0000000..4519d72 --- /dev/null +++ b/worldhistory/engine.py @@ -0,0 +1,223 @@ +"""The engine: one step = events → hazards → capacity → growth → conflict → migration → adaptation → +sub-species → tech. Later sub-projects add phases to Engine.phases.""" +import resource +import time +from pathlib import Path + +import numpy as np + +from . import tech as techm +from .adaptation import adapt, prospective +from .config import dumps_toml, resolved +from .conflict import apply_conflict, apply_curse, conflict +from .demography import capacity, crowding, grow, harvest_shocks, overlap_matrix, stochastic_round, veins, weather +from .events import apply_event +from .habitat import comfort, environment, fitness, native_optima, quality, suitability +from .hazards import make_hazard +from .migration import drift_and_bud, long_jumps, travel_cost +from .nudges import nudge_multipliers +from .regions import Regions +from .snapshot import engine_version, write_json, write_snapshot +from .state import new_state +from .lineage import (birth_step, displacement, inherit, init_rules, native, progress, raids, return_step, + rituals, summarise) + + +class Engine: + def __init__(self, world, history, races, seed=None): + self.w, self.h, self.races = world, history, races + self.R, self.n = len(races), world.n + self.seed = history.seed if seed is None else seed + self.regions = Regions(world, history.regions) + self.micro = np.exp(0.3 * world.noise(self.seed + 1000)) + self.state = new_state(self.R, self.n, native_optima(races), self.seed) + self.rules = init_rules(races) + self.returns, self.ritual_log, self.raid_log = [], [], [] + self.births = np.zeros((self.R, self.n)) + self.hazards = [make_hazard(s) for s in history.hazards] + self.events_log, self.emergence, self.stats_steps = [], [], [] + self.hazard_stats = [] + self.years_per_step = history.step + self._fc = {} # per race: inputs and values of _fit_comf + self._nudge_memo = None + self._world_changed() + self._capacity() + self.hostile = np.zeros((self.R, self.n)) + self.crowd = np.zeros((self.R, self.n)) + self.phases = [("events", self.phase_events), ("hazards", self.phase_hazards), + ("capacity", self._capacity), ("growth", self.phase_growth), ("conflict", self.phase_conflict), + ("migration", self.phase_migration), ("adaptation", self.phase_adaptation), + ("subspecies", self.phase_subspecies), ("tech", self.phase_tech)] + + # ---- derived state ---- + def _world_changed(self): + self.env = environment(self.w) + self.suit = np.stack([suitability(self.w, r) for r in self.races]) + self.q = np.stack([quality(self.suit[i], r, self.micro) for i, r in enumerate(self.races)]) + self.alpha = overlap_matrix(self.suit, self.races) + self.native_ok = np.stack([(self.suit[i] > 0) & (fitness(self.env, np.repeat(native(r)[:, None], self.n, 1), r) + >= self.h.f_min) for i, r in enumerate(self.races)]) + self.tcost0 = np.stack([travel_cost(self.w, r, 0.0) for r in self.races]) + + def _capacity(self): + st, tp = self.state, self.h.tech + prospective(st, self.w) + self.nudge_g, nudge_c, self.nudge_e, self.nudge_a = self._nudges(st.t) + e = weather(self.w, st.rng) + K = np.zeros((self.R, self.n)) + self.q_eff = np.zeros((self.R, self.n)) + self.eff = [] + self.OK = np.zeros((self.R, self.n), bool) + self.fit = np.zeros((self.R, self.n)) + self.nudge_share = {} + for r, race in enumerate(self.races): + eff = techm.effective(st.T[r], tp) + self.eff.append(eff) + q = self.q[r] + q_eff = np.where(q > 0, q + st.improve[r] * (race.qmax - q), 0.0) + fit, comf = self._fit_comf(r, race, st.O[r], techm.reach(eff, race, tp)) + self.fit[r] = fit * comf # displacement compares the whole position on the curve + shock = harvest_shocks(self.w, q_eff, techm.shock_scale(eff, tp), e) + k = capacity(self.w, race, q_eff, fit, comf, techm.k_gain(eff, tp), shock, nudge_c[r]) + ok = (self.suit[r] > 0) & (fit >= self.h.f_min) + K[r] = np.where(ok, k, 0.0) + self.q_eff[r], self.OK[r] = q_eff, ok + tot = K[r].sum() + self.nudge_share[race.id] = float((K[r] * (1 - 1 / nudge_c[r])).sum() / tot) if tot > 0 else 0.0 + displacement(K, st.P, self.fit, self.rules, self.races) + self.K = K + + def _nudges(self, t): + """nudge_multipliers, rebuilt only when the set of active nudges or the era changes (read-only arrays).""" + key = (tuple(i for i, nd in enumerate(self.h.nudges) if nd["years"][0] <= t < nd["years"][1]), self.w.era) + if self._nudge_memo is None or self._nudge_memo[0] != key: + arrs = nudge_multipliers(self.h.nudges, self.races, self.regions, t, self.n) + for a in arrs: + a.flags.writeable = False + self._nudge_memo = (key, arrs) + return self._nudge_memo[1] + + def _fit_comf(self, r, race, O, rch): + """fitness and comfort of race r, recomputed only in cells whose optima or reach changed since the last + step (both are per-cell functions of those and the environment): the same values, much less work.""" + rch = np.broadcast_to(np.asarray(rch), (self.n,)) # keeps its dtype (float32 from tech): same rounding + c = self._fc.get(r) + if c is None or c[0] is not self.env: + fit, comf = fitness(self.env, O, race, rch), comfort(O, race) + else: + _, O0, r0, fit, comf = c + ch = np.flatnonzero((O != O0).any(0) | (rch != r0)) + if len(ch): + fit[ch] = fitness(self.env, O[:, ch], race, rch[ch], cells=ch) + comf[ch] = comfort(O[:, ch], race) + self._fc[r] = (self.env, O.copy(), rch.copy(), fit, comf) + return fit, comf + + # ---- phases ---- + def phase_events(self): + for i, ev in enumerate(self.h.events): + if ev["year"] != self.state.t: + continue + log = apply_event(self.state, self.w, self.races, ev, self.regions, self.q, self.OK, self.h, i) + self.events_log.append(log) + if ev["kind"] == "era_switch": + self._world_changed() + + def phase_hazards(self): + self.hazard_stats = [hz.step(self.state, self.w, self.races, self.regions, self.state.t, self.state.rng) + for hz in self.hazards] + + def phase_growth(self): + st = self.state + self.crowd = crowding(st.P, self.K, self.alpha, self.races, self.h.competition["intolerance"], world=self.w) + self.births = grow(st.P, self.K, self.crowd, self.q_eff, self.races, self.nudge_g) + for r, race in enumerate(self.races): + if race.veins > 0: + st.P[r] = veins(st.P[r], self.K[r], race.veins, race.allee) + + def phase_conflict(self): + blame = np.zeros_like(self.state.P) + loss, self.hostile = conflict(self.state.P, self.q_eff, self.crowd, self.races, blame=blame) + apply_conflict(self.state, loss, self.races) + apply_curse(self.state, blame, self.h.competition["curse_decay"]) + + def phase_migration(self): + st, tp, mp = self.state, self.h.tech, self.h.migration + tcost = self.tcost0 - np.stack([techm.travel_bonus(self.eff[i], tp) for i in range(self.R)]) + danger = sum((hz.danger for hz in self.hazards if getattr(hz, "danger", None) is not None), np.zeros(self.n)) + drift_and_bud(st, self.w, self.races, self.K, self.OK, tcost, self.hostile + danger[None], self.nudge_e, st.rng, + mp) + boost = np.stack([techm.jump_boost(self.eff[i], tp) for i in range(self.R)]) * self.nudge_e + long_jumps(st, self.w, self.races, self.K, self.OK, st.rng, mp, boost, env=self.env, suit=self.suit, + f_min=self.h.f_min) + + def phase_adaptation(self): + adapt(self.state, self.env, self.races, self.q_eff, self.years_per_step, self.nudge_a) + + def phase_subspecies(self): + st, rng = self.state, self.state.rng + for e in birth_step(st, self.w, self.races, self.rules, self.regions, rng, self.native_ok): + e["lat"], e["lon"] = float(self.w.lat[e["cell"]]), float(self.w.lon[e["cell"]]) + self.emergence.append(e) + self.returns += summarise(return_step(st, self.w, self.races, self.rules, rng, self.native_ok), ("people",)) + inherit(st, self.w, self.races, self.rules, self.births, self.regions, self.native_ok) + self.raid_log += summarise(raids(st, self.w, self.races, self.rules, rng, self.native_ok), ("raids",)) + self.ritual_log += summarise(rituals(st, self.w, self.races, self.rules, rng), ("victims", "converts", "backlash", "dead")) + + def phase_tech(self): + steps = self.years_per_step / 10.0 + techm.step_tech(self.state, self.w, self.races, self.h.tech, steps) + techm.step_land(self.state, self.races, self.h.tech, steps) + stochastic_round(self.state.P, self.state.rng) + + # ---- driving ---- + def step(self): + """Process year state.t, then advance.""" + for _, fn in self.phases: + fn() + self.state.t += self.h.step + + def _stats(self, t): + P = self.state.P + ids = [r.id for r in self.races] + fam = {} + for i, r in enumerate(self.races): + fam[r.family] = fam.get(r.family, 0.0) + float(P[i].sum()) + return {"year": t, "pop": {ids[i]: float(P[i].sum()) for i in range(self.R)}, + "family": fam, + "occupied": {ids[i]: int((P[i] >= 1).sum()) for i in range(self.R)}, + "regions": {name: {ids[i]: float(P[i, self.regions(name)].sum()) for i in range(self.R)} + for name in self.h.stats_regions}, + "hazards": self.hazard_stats, "nudge_capacity": self.nudge_share, + "progress": {f"{self.races[r['p']].id}>{self.races[r['s']].id}": self._mean_progress(r) + for r in self.rules}} + + def _mean_progress(self, rule): + p, s = rule["p"], rule["s"] + P = self.state.P[p] + tot = P.sum() + return float((progress(self.state.O[p], self.races[p], self.races[s]) * P).sum() / tot) if tot > 0 else 0.0 + + def run(self, out_dir, years=None, progress=None): + years = self.h.years if years is None else years + out = Path(out_dir) + (out / "snap").mkdir(parents=True, exist_ok=True) + meta = {"seed": self.seed, "engine": engine_version(), "world": str(self.w.meta.get("dir", "")), + "years": years} + (out / "run.toml").write_text(dumps_toml(resolved(self.h, self.races, **meta))) + t0 = time.time() + while self.state.t <= years: + t = self.state.t + self.step() + self.stats_steps.append(self._stats(t)) + if t % self.h.snapshot_every == 0 or t == years: + write_snapshot(out / "snap" / f"y{t:04d}.npz", self.state, t, self.races) + if progress and t % 500 == 0: + progress(f"year {t}: " + ", ".join(f"{k} {v:,.0f}" for k, v in self.stats_steps[-1]["pop"].items())) + stats = {"steps": self.stats_steps, "events": self.events_log, "emergence": self.emergence, + "returns": self.returns, "rituals": self.ritual_log, "raids": self.raid_log, + "timing": {"total_s": time.time() - t0, + "peak_rss_mb": resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024}, + "meta": meta} + write_json(out / "stats.json", stats) + return stats diff --git a/worldhistory/events.py b/worldhistory/events.py new file mode 100644 index 0000000..5bd6294 --- /dev/null +++ b/worldhistory/events.py @@ -0,0 +1,68 @@ +"""Scheduled events: seed (the gifting), cull / die_off (a share of people in a region), era_switch.""" +import numpy as np + + +def seed(st, world, races, q, OK, fringe, rng, only=None, mask=None, clusters=None, heads=None): + """Clusters per race: most by q², a share on the fringe (0.05 < q < 0.3) to push them. Targeted (only = a race + id): that race's clusters within mask, by q², no fringe.""" + placed = {} + for r, race in enumerate(races): + if only is not None and race.id != only: + continue + n_cl = race.clusters if clusters is None else int(clusters) + size = race.heads if heads is None else float(heads) + if n_cl <= 0: + continue + qq = q[r] * OK[r] * (1.0 if mask is None else mask) + if race.seed_realm: + ocean = np.asarray(world.fields["ocean"], bool) + qq = qq * (ocean if race.seed_realm == "sea" else ~ocean) + n_fringe = 0 if only is not None else int(round(n_cl * fringe)) + n_main = n_cl - n_fringe + picks = [] + p = qq ** 2 + if n_main and p.sum() > 0: + k = min(n_main, int((p > 0).sum())) + picks += list(rng.choice(world.n, size=k, replace=False, p=p / p.sum())) + fr = np.flatnonzero((qq > 0.05) & (qq < 0.3)) + if n_fringe and len(fr): + picks += list(rng.choice(fr, size=min(n_fringe, len(fr)), replace=False)) + np.add.at(st.P[r], np.array(picks, np.int64), size) + placed[race.id] = [int(c) for c in picks] + return placed + + +def cull(st, races, mask, share, by_race, noise=None, amp=0.0): + extra = amp * noise if noise is not None else 0.0 + for r, race in enumerate(races): + s = np.clip(by_race.get(race.id, share) + extra, 0, 1) + st.P[r] = np.where(mask, st.P[r] * (1 - s), st.P[r]) + + +def era_switch(st, world, races, era): + world.set_era(era) + sea = np.asarray(world.fields["ocean"], bool) + for r, race in enumerate(races): + if race.realm == "land": + st.P[r, sea] = 0.0 + elif race.realm == "sea": + st.P[r, ~sea] = 0.0 + + +def apply_event(st, world, races, ev, regions, q, OK, history, index): + before = {r.id: float(st.P[i].sum()) for i, r in enumerate(races)} + kind = ev["kind"] + info = {} + if kind == "seed": + if ev.get("race"): + info["placed"] = seed(st, world, races, q, OK, history.fringe, st.rng, only=ev["race"], + mask=regions(ev["region"]), clusters=ev.get("clusters"), heads=ev.get("heads")) + else: + info["placed"] = seed(st, world, races, q, OK, history.fringe, st.rng) + elif kind in ("cull", "die_off"): + noise = world.noise(history.seed * 7919 + index) if ev["noise"] else None + cull(st, races, regions(ev["region"]), ev["share"], ev["by_race"], noise, ev["noise"]) + elif kind == "era_switch": + era_switch(st, world, races, ev["era"]) + after = {r.id: float(st.P[i].sum()) for i, r in enumerate(races)} + return {"year": int(ev["year"]), "kind": kind, "region": ev["region"], "before": before, "after": after, **info} diff --git a/worldhistory/habitat.py b/worldhistory/habitat.py new file mode 100644 index 0000000..13c5ea5 --- /dev/null +++ b/worldhistory/habitat.py @@ -0,0 +1,111 @@ +"""Habitat: environment conditions per cell, suitability, quality, tolerance fitness and comfort.""" +import importlib +from dataclasses import dataclass + +import numpy as np + +from .config import CONDITIONS +from .predicates import depth, evaluate + + +@dataclass +class Env: + x: dict # condition -> value per cell + lo: dict # condition -> low end of the cell's span (toward neighbour midpoints) + hi: dict + + +def intervals(world, x): + """A cell spans its own value out to the midpoints toward its neighbours (cells are big: ~85 km at r4).""" + mids = np.where(world._nbv, (x[:, None] + x[world._nbj]) / 2, x[:, None]) + return np.minimum(x, mids.min(1)), np.maximum(x, mids.max(1)) + + +def environment(world, F=None): + F = world.fields if F is None else F + sea = np.asarray(F["ocean"], bool) + x = {"depth": depth(F), + "air_pressure": np.asarray(F["pressure_bar"], float), + "gravity": np.asarray(F["gravity_g"], float), + "o2": np.asarray(F["o2_fraction"], float), + "temperature": np.where(sea, F["bottom_temp_c"], F["T_mean"]).astype(float), + "rain": np.log10(np.maximum(np.asarray(F["P_ann"], float), 1.0))} + lo, hi = {}, {} + for c in CONDITIONS: + lo[c], hi[c] = intervals(world, x[c]) + lo["rain"] = np.where(sea, -np.inf, lo["rain"]) # rain means nothing to sea dwellers + hi["rain"] = np.where(sea, np.inf, hi["rain"]) + return Env(x, lo, hi) + + +def suitability(world, race, F=None): + F = world.fields if F is None else F + s = np.zeros(world.n) + for t in race.terms: + s += t["w"] * evaluate(t["p"], world, F) + if race.hook: + mod, fn = race.hook.split(":") + s += np.asarray(getattr(importlib.import_module(mod), fn)(F), float) + for g in race.gates: + realm = g.get("realm", "both") + where = np.ones(world.n, bool) if realm == "both" else evaluate(realm, world, F) > 0.5 + m = (evaluate(g["require"], world, F) if "require" in g + else 1 - (1 - g["factor"]) * evaluate(g["penalty"], world, F)) + s = s * np.where(where, m, 1.0) + sea = np.asarray(F["ocean"], bool) + if race.realm == "land": + s = s * ~sea + elif race.realm == "sea": + s = s * sea + return np.clip(s, 0, None) + + +def quality(s, race, micro): + """0..qmax: suitability (times fixed micro-habitat noise) over the race's own 99th percentile.""" + v = s * micro + pos = v[v > 0] + ref = np.percentile(pos, 99) if len(pos) else 1.0 + return race.qmax * np.clip(v / ref, 0, 1) + + +def native_optima(races): + return np.array([[r.tolerance[c].optimum if c in r.tolerance else 0.0 for c in CONDITIONS] for r in races], float) + + +def fitness(env, O_r, race, reach=0.0, cells=None): + """cells: judge O_r (C, len(cells)) against those cells' conditions instead of the whole world. Lopsided bell: + width_hi for cells above the optimum, width_lo below.""" + f = np.ones(O_r.shape[1]) + for ci, c in enumerate(CONDITIONS): + tol = race.tolerance.get(c) + if tol is None: + continue + o = O_r[ci] + lo, hi = (env.lo[c], env.hi[c]) if cells is None else (env.lo[c][cells], env.hi[c][cells]) + above = np.maximum(lo - o, 0) / (tol.width_hi * (1 + reach)) + below = np.maximum(o - hi, 0) / (tol.width_lo * (1 + reach)) + f = f * np.exp(-0.5 * (above ** 2 + below ** 2)) + return f + + +STRAIN_AT_EDGE = 0.35 # core spec §3.4: fully adapted to the range edge = low gear + + +def comfort(O_r, race): + """Peak height of each adapted curve: 1 at the race's native optimum, falling smoothly to STRAIN_AT_EDGE at the + edge of its adaptable range (lower beyond) — adapted groups are livable, not good (author 2026-09-30). A + `comfort` list in config overrides this for its condition.""" + m = np.ones(O_r.shape[1]) + for ci, c in enumerate(CONDITIONS): + tol = race.tolerance.get(c) + if tol is None: + continue + if tol.comfort: + xs, ys = zip(*sorted(tuple(p) for p in tol.comfort)) + m = m * np.interp(O_r[ci], xs, ys) + continue + o = O_r[ci] + span = np.where(o >= tol.optimum, tol.hi - tol.optimum, tol.optimum - tol.lo) + x = np.where(span > 0, (o - tol.optimum) / np.maximum(span, 1e-12), 0.0) + m = m * STRAIN_AT_EDGE ** (x * x) + return m diff --git a/worldhistory/hazards.py b/worldhistory/hazards.py new file mode 100644 index 0000000..78c8592 --- /dev/null +++ b/worldhistory/hazards.py @@ -0,0 +1,108 @@ +"""Hazards: static mortality layers and roaming units (few, very powerful creatures that raid).""" +import numpy as np + + +class Static: + def __init__(self, spec): + self.spec = spec + self.danger = None # cells drift avoids (push hazards) + + def step(self, st, world, races, regions, t, rng): + s = self.spec + y0, y1 = s["years"] + if not (y0 <= t < y1): + self.danger = np.zeros(world.n) if s["push"] > 0 else None + return {"name": s["name"]} + f = s["decay"] ** ((t - y0) / 10.0) + m = s["mortality"] * f + mask = regions(s["region"]) + st.P[:, mask] *= 1 - m + if s["push"] > 0: # on land, a share flees each step (e.g. megafauna driving people into the sea) + land = mask & ~np.asarray(world.fields["ocean"], bool) + st.push[:, land] += s["push"] * f * st.P[:, land] + self.danger = land.astype(float) + return {"name": s["name"], "mortality": m} + + +class Roaming: + def __init__(self, spec): + self.spec = spec + self.home = self.cur = self.kind = None + + def _start(self, st, world, races, regions, rng): + s = self.spec + region = regions(s["home_region"]) + ids = [r.id for r in races] + w = np.zeros(world.n) + if s["home_near"] and s["home_near"] in ids: + w = st.P[ids.index(s["home_near"])] * region + if w.sum() <= 0: + w = (region & ~np.asarray(world.fields["ocean"], bool)).astype(float) + if w.sum() <= 0: + w = region.astype(float) + n = s["local"] + s["wanderers"] + self.home = rng.choice(world.n, size=n, p=w / w.sum()).astype(np.int64) + self.cur = self.home.copy() + self.kind = np.r_[np.zeros(s["local"], np.int8), np.ones(s["wanderers"], np.int8)] + + def _burst(self, st, world, cell, share): + cells = world.within(int(cell), self.spec["radius_km"]) + st.P[:, cells] *= 1 - share + + def _walk(self, world, cell, km, rng): + v = world.xyz[cell] + rnd = rng.normal(size=3) + t = rnd - rnd.dot(v) * v + t /= np.linalg.norm(t) + d = km / world.radius_km + return int(world.nearest(v * np.cos(d) + t * np.sin(d))) + + def step(self, st, world, races, regions, t, rng): + s = self.spec + if t < s["start"]: + return {"name": s["name"], "units": 0, "raids": 0, "births": 0, "deaths": 0} + if self.home is None: + self._start(st, world, races, regions, rng) + tot = st.P.sum(0) + occ = np.flatnonzero(regions(s["home_region"]) & (tot >= 1)) + raids, dead = 0, [] + for i in range(len(self.home)): + if self.kind[i] == 0: + self.cur[i] = self.home[i] + if len(occ) and rng.random() < s["raid_chance"]: + tgt = int(rng.choice(occ, p=tot[occ] / tot[occ].sum())) + self.cur[i] = tgt + self._burst(st, world, tgt, s["raid_burst"]) + raids += 1 + dens = tot[tgt] / world.area[tgt] + if rng.random() < s["raid_death"] * min(1.0, dens / s["density_ref"]): + dead.append(i) + else: + self._burst(st, world, self.home[i], s["local_mortality"]) + else: + self.cur[i] = self._walk(world, self.cur[i], s["wander_step_km"], rng) + if rng.random() < s["wander_stop"]: + self._burst(st, world, self.cur[i], s["wander_burst"]) + n = len(self.home) + dead = set(dead) | {i for i in range(n) if rng.random() < s["infight"]} + keep = np.array([i not in dead for i in range(n)], bool) + self.home, self.cur, self.kind = self.home[keep], self.cur[keep], self.kind[keep] + births = 0 + lo, hi = s["clutch"] + local = np.flatnonzero(self.kind == 0) + if len(self.home) and s["growth"] > 0: + n = len(self.home) + room = max(0.0, 1.0 - n / s["cap"]) if s["cap"] > 0 else 1.0 # logistic when capped + for _ in range(rng.poisson(s["growth"] * n * room / ((lo + hi) / 2))): + size = int(rng.integers(lo, hi + 1)) + at = self.home[rng.choice(local)] if len(local) else self.home[rng.integers(len(self.home))] + self.home = np.r_[self.home, np.full(size, at)] + self.cur = np.r_[self.cur, np.full(size, at)] + self.kind = np.r_[self.kind, np.zeros(size, np.int8)] + births += size + return {"name": s["name"], "units": int(len(self.home)), "raids": raids, "births": births, + "deaths": int(len(dead))} + + +def make_hazard(spec): + return Static(spec) if spec["kind"] == "static" else Roaming(spec) diff --git a/worldhistory/lineage.py b/worldhistory/lineage.py new file mode 100644 index 0000000..6fc6a81 --- /dev/null +++ b/worldhistory/lineage.py @@ -0,0 +1,316 @@ +"""Sub-species as birth events (spec 2026-09-30 gradual adaptation §2). A group's progress toward a related race is +read from its adapted optima; births (once per sub-species) and returns (any number) happen with a chance that +rises with progress; the born group takes the new race's native curves.""" +import numpy as np + +from .config import CONDITIONS +from .state import convert, move + + +def native(race): + return np.array([race.tolerance[c].optimum if c in race.tolerance else 0.0 for c in CONDITIONS], float) + + +def progress(O, r_from, r_to): + """0..1 per group: how far its optima have moved from r_from's native curves toward r_to's (spec §2.2).""" + num, den = np.zeros(O.shape[1]), 0.0 + for ci, c in enumerate(CONDITIONS): + a, b = r_from.tolerance.get(c), r_to.tolerance.get(c) + if a is None or b is None: + continue + gap = b.optimum - a.optimum + w = a.side(gap > 0) + if abs(gap) <= 0.25 * w: + continue + weight = abs(gap) / w + num += weight * np.clip((O[ci] - a.optimum) / gap, 0.0, 1.0) + den += weight + return num / den if den > 0 else num + + +def smoothstep(x): + x = np.clip(x, 0.0, 1.0) + return x * x * (3 - 2 * x) + + +def chance(a, lo, hi, rate): + a = np.asarray(a, float) + if hi <= lo: + return np.where(a >= hi, rate, 0.0) + return np.where(a < lo, 0.0, rate * smoothstep((a - lo) / (hi - lo))) + + +def init_rules(races): + idx = {r.id: i for i, r in enumerate(races)} + return [{"s": s, "p": idx[r.emerge["parent"]], "origin": None, "year": None, "peak": None, "peak_t": None, + "victims": 0.0, "backlash": False, "rate": r.emerge["ritual_rate"]} + for s, r in enumerate(races) if r.emerge is not None] + + +def _crisis(rule, spec, P, t): + """Parent groups almost eradicated: P ≤ (1 − crisis_drop) of their recent peak (memory crisis_years, decaying), + with a peak of at least crisis_min people. Updates the rule's peak record.""" + if rule["peak"] is None: + rule["peak"] = P.copy() + else: + rule["peak"] = np.maximum(P, rule["peak"] * np.exp(-(t - rule["peak_t"]) / spec["crisis_years"])) + rule["peak_t"] = t + return (P > 0) & (P <= (1 - spec["crisis_drop"]) * rule["peak"]) & (rule["peak"] >= spec["crisis_min"]) + + +def _where(spec, world, regions, F): + m = regions(spec["region"]) if spec["region"] else np.ones(world.n, bool) + cond = spec["condition"] + if cond: + x = np.asarray(F[cond["field"]], float) + if "below" in cond: + m = m & (x < cond["below"]) + if "above" in cond: + m = m & (x > cond["above"]) + return m + + +def _shelter(world, where, ok, km): + """Per cell: for crisis cells (where), the nearest cell within km the new race can live in (ok), else −1; km 0 → the + crisis cell itself must fit.""" + dest = np.full(world.n, -1, np.int64) + for c in np.flatnonzero(where): + near = world.within(c, km) + near = near[ok[near]] + if len(near): + dest[c] = near[np.argmin(world.km(np.full(len(near), c), near))] + return dest + + +def _shifted(st, p, s, cells, races, shift): + """Optima the born people arrive with: `shift` of the way from their own to the new race's native curves.""" + return (1 - shift) * st.O[p][:, cells] + shift * native(races[s])[:, None] + + +def birth_step(st, world, races, rules, regions, rng, ok, F=None): + F = world.fields if F is None else F + new = [] + for rule in rules: + if rule["origin"] is not None: + continue + s, p = rule["s"], rule["p"] + spec, P = races[s].emerge, st.P[p] + crisis = _crisis(rule, spec, P, st.t) if spec["crisis_drop"] > 0 else True + if st.t < spec["after"]: + continue + ritual = spec["mode"] == "ritual" + where = (P > 0) & _where(spec, world, regions, F) & crisis + if not ritual: + where &= ok[s] + if not where.any(): + continue + a = progress(st.O[p], races[p], races[s]) + if ritual: + dest = _shelter(world, where, ok[s], spec["crisis_km"]) + where &= dest >= 0 + if not where.any() or rng.random() >= spec["trigger"]: + continue + h = np.where(where, P, 0.0) + else: + h = np.where(where, chance(a, spec["birth_min"], spec["birth_sure"], spec["birth_rate"]) + * P / (P + races[p].founder), 0.0) + if spec["settlement_density"] > 0 and h.any(): + dens = np.where(h > 0, P / world.area, -1.0) + best = int(np.argmax(dens)) + h = np.where(np.arange(world.n) == best, h, 0.0) * (dens[best] >= spec["settlement_density"]) + if not h.any() or rng.random() >= 1 - np.prod(1 - np.clip(h, 0, 1)): + continue + c = np.flatnonzero(h > 0) + wgt = h[c] + if spec["isolated"]: + around = world.nb_sum(st.P.sum(0))[c] + wgt = wgt * (P[c] / (P[c] + around)) ** 4 + cell = int(rng.choice(c, p=wgt / wgt.sum())) + convert(st, p, s, [cell], [spec["convert"] * P[cell]], O_new=_shifted(st, p, s, [cell], races, + spec["birth_shift"])) + ev = {"race": races[s].id, "cell": cell, "year": int(st.t), "a": float(a[cell])} + if ritual: # survivors of the crisis flee to the nearest land they can live in + ev["cell"], ev["crisis"] = int(dest[cell]), cell + if ev["cell"] != cell: + move(st, s, [cell], [ev["cell"]], [st.P[s, cell]]) + rule["origin"], rule["year"] = ev["cell"], int(st.t) + new.append(ev) + return new + + +def return_step(st, world, races, rules, rng, ok): + """Sub-species groups adapted back toward the parent give birth to parent-race people (spec §2.6).""" + out = [] + for rule in rules: + s, p = rule["s"], rule["p"] + spec = races[s].emerge + if rule["origin"] is None or spec["mode"] == "ritual": + continue + P = st.P[s] + a = progress(st.O[s], races[s], races[p]) + h = np.where((P > 0) & ok[p], chance(a, spec["return_min"], spec["return_sure"], spec["return_rate"]) + * P / (P + races[s].founder), 0.0) + cells = np.flatnonzero(rng.random(world.n) < h) + if not len(cells): + continue + amt = spec["convert"] * P[cells] + convert(st, s, p, cells, amt, O_new=_shifted(st, s, p, cells, races, spec["birth_shift"])) + out += [{"race": races[p].id, "from": races[s].id, "cell": int(c), "year": int(st.t), "people": float(x)} + for c, x in zip(cells, amt)] + return out + + +def _born(rules): + return [r for r in rules if r["origin"] is not None] + + +def displacement(K, P, fit, rules, races): + """In cells a race shares with its parent / sub-species, the less fit one keeps (1 − displace·other's share) + of its capacity (spec §2.4).""" + for rule in _born(rules): + s, p = rule["s"], rule["p"] + d = races[s].emerge["displace"] + tot = P[s] + P[p] + share_s = np.where(tot > 0, P[s] / np.maximum(tot, 1e-12), 0.0) + both = (P[s] > 0) & (P[p] > 0) + K[p] *= np.where(both & (fit[p] < fit[s]), 1 - d * share_s, 1.0) + K[s] *= np.where(both & (fit[s] < fit[p]), 1 - d * (1 - share_s), 1.0) + + +def inherit(st, world, races, rules, births, regions, ok): + """Dominant trait: a share of the parent groups' new births in or next to the sub-species is born as it, gated + by the parents' own adaptation (spec §2.5; ritual races: no gate, §2.7).""" + for rule in _born(rules): + s, p = rule["s"], rule["p"] + spec = races[s].emerge + Ps, Pp = st.P[s], st.P[p] + c = np.flatnonzero((Pp > 0) & (births[p] > 0) & ok[s]) # only where parents are born and the new + if not len(c): # curves fit: compute there + continue + sub, cols = world.nb_local(c) + near_s = Ps[c] + world.nb_apply(sub, Ps[cols]) + near = near_s + Pp[c] + world.nb_apply(sub, Pp[cols]) + share = np.where(near > 0, near_s / np.maximum(near, 1e-12), 0.0) + where = near_s > 0 + if spec["spread"] == "region" and spec["region"]: + where &= regions(spec["region"])[c] + if spec["mode"] == "ritual": + gate = np.ones(len(c)) + else: + a = progress(st.O[p][:, c], 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, c], 0.0) + keep = amt > 0 + cells, amt = c[keep], amt[keep] + if len(cells): + convert(st, p, s, cells, amt, O_new=np.repeat(native(races[s])[:, None], len(cells), 1)) + + +def raids(st, world, races, rules, rng, ok): + """Before the smackdown, ritual races send out cells: groups of at least raid_min send ~Poisson(raid_rate) raids of + raid_size people each to occupied non-ritual land within raid_km (weighted by its people) where the race can + live; the rituals there do the killing. Each group keeps at least raid_min / 2.""" + out = [] + for rule in _born(rules): + s = rule["s"] + spec = races[s].emerge + if spec["mode"] != "ritual" or spec["raid_rate"] <= 0 or rule["backlash"]: + continue + others = np.sum(np.delete(st.P, s, axis=0), axis=0) + for c in np.flatnonzero(st.P[s] >= spec["raid_min"]): + k = min(rng.poisson(spec["raid_rate"]), int((st.P[s, c] - spec["raid_min"] / 2) // spec["raid_size"])) + if k <= 0: + continue + near = world.within(c, spec["raid_km"]) + near = near[(near != c) & (others[near] >= 1) & ok[s][near]] + if not len(near): + continue + to = rng.choice(near, size=k, p=others[near] / others[near].sum()) + move(st, s, np.full(k, c), to, np.full(k, float(spec["raid_size"]))) + out.append({"race": races[s].id, "cell": int(c), "year": int(st.t), "raids": k}) + return out + + +def rituals(st, world, races, rules, rng): + """Ritual races win converts at a blood price: rituals per cell ~ Poisson(rate·√people); each kills ritual_cost + and converts ritual_converts, taken from the neighbour cell (or own) with most parents. victims = "all": the dead + come from every nearby non-ritual people by presence. After backlash_victims dead in all, the neighbours strike + back once: the race loses backlash_loss of its people and holds rituals at backlash_calm × the rate.""" + out = [] + for rule in _born(rules): + s, p = rule["s"], rule["p"] + spec = races[s].emerge + if spec["mode"] != "ritual": + continue + others = np.array([r for r in range(len(races)) if r != s]) if spec["victims"] == "all" else np.array([p]) + per = spec["ritual_cost"] + spec["ritual_converts"] + for c in np.flatnonzero(st.P[s] >= 1): + k = rng.poisson(rule["rate"] * st.P[s, c] ** spec["ritual_power"]) + if k == 0: + continue + nb = world.nb[c] + opts = np.r_[c, nb[nb >= 0]] + v = int(opts[np.argmax(st.P[p, opts])]) + pool = st.P[np.ix_(others, opts)] + k = min(k, int(st.P[p, v] // per), int(pool.sum() // per)) + if k == 0: + continue + dead = k * spec["ritual_cost"] + taken = dead * pool / pool.sum() + st.P[np.ix_(others, opts)] = pool - taken + conv = k * spec["ritual_converts"] + st.P[p, v] -= conv + st.P[s, c] += conv # the converts join the ritual group, taking its curves (unchanged O) + rule["victims"] += dead + out.append({"race": races[s].id, "cell": int(c), "year": int(st.t), "victims": float(dead), + "converts": float(conv), + "dead": {races[o].id: float(x) for o, x in zip(others, taken.sum(1)) if x > 0}}) + by_count = spec["backlash_victims"] > 0 and rule["victims"] >= spec["backlash_victims"] + by_time = spec["backlash_years"] > 0 and st.t >= rule["year"] + spec["backlash_years"] + if not rule["backlash"] and (by_count or by_time): + rule["backlash"] = True + lost = _smackdown(st, world, s, others, spec["backlash_loss"]) + rule["rate"] = spec["ritual_rate"] * spec["backlash_calm"] + out.append({"race": races[s].id, "cell": -1, "year": int(st.t), "victims": 0.0, "converts": 0.0, + "backlash": lost}) + return out + + +def _smackdown(st, world, s, others, loss): + """Kill `loss` of race s; each group's share ∝ its exposure (own people + non-ritual people in and next to its + cell), capped at all of it: big, exposed groups are wiped out, small remote ones survive. Returns people lost.""" + P = st.P[s] + enemies = st.P[others].sum(0) + h = np.where(P > 0, P + enemies + world.nb_sum(enemies), 0.0) + target = loss * P.sum() + if target <= 0 or not h.any(): + return 0.0 + lo, hi = 0.0, 1.0 / h[h > 0].min() # at hi every group is wiped out + for _ in range(200): + k = (lo + hi) / 2 + lo, hi = (k, hi) if (P * np.minimum(1.0, k * h)).sum() < target else (lo, k) + frac = np.minimum(1.0, hi * h) + st.P[s] = P * (1 - frac) + return float(target) + + +def summarise(events, sums): + """Per-cell events of one step → one record per race (and source race, if given): groups and summed `sums` + fields; dict fields are summed per key (keeps stats small).""" + out = {} + for e in events: + k = (e["race"], e.get("from"), e["year"]) + head = {"race": e["race"], **({"from": e["from"]} if "from" in e else {}), "year": e["year"], "groups": 0} + rec = out.setdefault(k, {**head, **{f: 0.0 for f in sums}}) + rec["groups"] += e["cell"] >= 0 # cell −1: a race-wide record (the backlash), not a group + for f in sums: + if f not in e: + continue + x = e[f] + if isinstance(x, dict): + rec[f] = rec[f] or {} + for key, v in x.items(): + rec[f][key] = rec[f].get(key, 0.0) + v + else: + rec[f] += x + return list(out.values()) diff --git a/worldhistory/migration.py b/worldhistory/migration.py new file mode 100644 index 0000000..15fbb4f --- /dev/null +++ b/worldhistory/migration.py @@ -0,0 +1,149 @@ +"""Migration: drift toward headroom between occupied cells (ideal free distribution, Fretwell & Lucas 1970, by +logit choice, McFadden 1974), budding of founder groups into empty cells (a travelling front: Fisher 1937; +Ammerman & Cavalli-Sforza 1984), and fat-tailed long jumps (Kot, Lewis & van den Driessche 1996; Clark 1998).""" +import numpy as np + +from .habitat import fitness +from .predicates import evaluate +from .state import move + + +V_REF = 0.5 # m/s: a current this strong counts as "1" for current_bias +PATH_FRAC = np.array([0.0, 0.25, 0.5, 0.75, 1.0]) # points along a jump where the current is read + + +def _downstream(world, v, d, rng, bias, k=4): + """Per jumper, one tangent direction of k random candidates, weighted exp(bias·a), a = mean current along the + great-circle path (5 points) projected on the heading, / V_REF.""" + cur = world.current + m = len(v) + rnd = rng.normal(size=(k, m, 3)) + t = rnd - np.sum(rnd * v[None], axis=-1, keepdims=True) * v[None] + t /= np.linalg.norm(t, axis=-1, keepdims=True) + s = d[None, :, None] * PATH_FRAC[None, None, :] # (1, m, 5) + pts = v[None, :, None] * np.cos(s)[..., None] + t[:, :, None] * np.sin(s)[..., None] + head = -v[None, :, None] * np.sin(s)[..., None] + t[:, :, None] * np.cos(s)[..., None] + c = cur[world.nearest(pts.reshape(-1, 3))].reshape(k, m, len(PATH_FRAC), 3) + a = np.sum(c * head, axis=-1).mean(axis=-1) / V_REF # (k, m) + wgt = np.exp(bias * (a - a.max(axis=0, keepdims=True))) + pick = (rng.random(m)[None] * wgt.sum(0) < np.cumsum(wgt, axis=0)).argmax(axis=0) + return t[pick, np.arange(m)] + + +def travel_cost(world, race, bonus, F=None): + F = world.fields if F is None else F + tr = race.travel + river = np.asarray(F["river"], bool) | np.asarray(F["lake"], bool) + coast = evaluate("coast", world, F) + c = np.minimum(tr["river"] * river, tr["coast"] * coast) + c = c + tr["mountain"] * evaluate("mountain", world, F) + tr["desert"] * evaluate("desert", world, F) + c = c + tr["sea"] * evaluate("sea", world, F) + return c - bonus + + +def drift_and_bud(st, world, races, K, OK, tcost, hostile, expansion, rng, mp): + for r, race in enumerate(races): + P = st.P[r] + push = st.push[r].copy() + if race.group_max > 0: # big groups split: half moves off + push = push + np.where(P > race.group_max, P / 2, 0.0) + st.push[r] = 0.0 + act = np.flatnonzero(P >= 1) + if not len(act): + continue + Pa, Ka = P[act], K[r, act] + crowd = np.clip(Pa / np.maximum(Ka, 1e-9), 0, 2) + drift = np.minimum(race.mobility * expansion[r, act] * (1 + crowd) * Pa, np.maximum(Pa - race.allee, 0.0)) + out = np.minimum(drift + push[act], 0.8 * Pa) # drift never takes a viable group below the Allee size + j = world._nbj[act] + valid = world._nbv[act] & OK[r][j] + Kj, Hj = K[r][j], hostile[r][j] # attraction only where it is read + att_j = (Kj - P[j]) / (Kj + 1) - mp["hostile"] * Hj + a = mp["temp"] * (att_j - tcost[r][j]) + if race.current_bias and world.current.any(): # sea edges: go with the current + sea = np.asarray(world.fields["ocean"], bool) + along = np.einsum("mk,mjk->mj", world.current[act], world._nbt[act]) / V_REF + a = a + race.current_bias * np.where(sea[act][:, None] & sea[j], along, 0.0) + amax = np.max(np.where(valid, a, -1e300), axis=1, keepdims=True) + with np.errstate(over="ignore", invalid="ignore"): + wgt = np.where(valid, np.exp(np.where(valid, a - amax, 0.0)), 0.0) + s = wgt.sum(1, keepdims=True) + wgt = np.where(s > 0, wgt / np.maximum(s, 1e-300), 0.0) + moved = out[:, None] * wgt + empty = P[j] < 1 + pushf = np.clip(crowd / 0.3, 0, 1)[:, None] + with np.errstate(over="ignore"): + pf = np.clip(mp["pf"] * pushf * (6 * wgt) ** 1.5 * np.exp(-tcost[r][j]) * expansion[r, act][:, None], 0, 1) + found = (rng.random(moved.shape) < pf) & (Pa[:, None] >= 2 * race.founder) + flee = np.minimum(push[act], 0.8 * Pa)[:, None] * wgt # driven-out people settle empty land too + moved = np.where(empty, np.maximum(np.where(found, race.founder, 0.0), flee), moved) + moved = np.where(valid, np.minimum(moved, Pa[:, None] / 6), 0.0) + move(st, r, np.repeat(act, j.shape[1]), j.ravel(), moved.ravel()) + gather(st, world, r, race, OK) + if race.roam > 0: + roam_groups(st, world, r, race, OK, rng) + + +def roam_groups(st, world, r, race, OK, rng): + """Caravans: with chance roam per step, a whole group moves on together to one random livable neighbour.""" + P = st.P[r] + act = np.flatnonzero(P >= 1) + go = act[rng.random(len(act)) < race.roam] + if not len(go): + return + j = world._nbj[go] + valid = world._nbv[go] & OK[r][j] + pick = np.where(valid, rng.random(valid.shape), -1.0).argmax(1) + ok = valid[np.arange(len(go)), pick] + move(st, r, go[ok], j[np.arange(len(go)), pick][ok], P[go[ok]]) + + +def gather(st, world, r, race, OK): + """Groups below the Allee size join their largest viable neighbour instead of dwindling alone.""" + P = st.P[r] + small = np.flatnonzero((P >= 1) & (P < race.allee)) + if not len(small): + return + j = world._nbj[small] + Pj = np.where(world._nbv[small] & OK[r][j] & (P[j] >= race.allee), P[j], -1.0) + best = Pj.argmax(1) + ok = Pj[np.arange(len(small)), best] > 0 + move(st, r, small[ok], j[np.arange(len(small)), best][ok], P[small[ok]]) + + +def long_jumps(st, world, races, K, OK, rng, mp, boost, env=None, suit=None, f_min=0.1): + """env/suit given: a landing cell is judged for the migrants' own adapted optima (they can skip ground they + could not live on, e.g. swim an abyss), not by the cell's stored optimum.""" + for r, race in enumerate(races): + P = st.P[r] + crowd = P / np.maximum(K[r], 1e-9) + cand = np.flatnonzero((P >= 2 * max(race.allee, race.founder)) & (crowd > 0.5)) + if not len(cand): + continue + go = cand[rng.random(len(cand)) < mp["jump_rate"] * race.jump * boost[r, cand]] + if not len(go): + continue + u = 1.0 - rng.random(len(go)) + d = np.minimum(mp["jump_xm_km"] * u ** (-1.0 / mp["jump_alpha"]), mp["jump_cap_km"]) / world.radius_km + v = world.xyz[go] + if race.jump_path == "sea" and race.current_bias and world.current.any(): + t = _downstream(world, v, d, rng, race.current_bias) + else: + rnd = rng.normal(size=v.shape) + t = rnd - (rnd * v).sum(1, keepdims=True) * v + t /= np.linalg.norm(t, axis=1, keepdims=True) + tgt = world.nearest(v * np.cos(d)[:, None] + t * np.sin(d)[:, None]) + wet = np.ones(len(go), bool) + if race.jump_path == "sea": # swimmers: every point along the way must be water + ocean = np.asarray(world.fields["ocean"], bool) + k = int(np.ceil(d.max() * world.radius_km / 25.0)) + 1 + s = d[:, None] * np.linspace(0, 1, k)[None] + pts = v[:, None] * np.cos(s)[..., None] + t[:, None] * np.sin(s)[..., None] + wet = ocean[world.nearest(pts.reshape(-1, 3))].reshape(len(go), k).all(1) + if env is None: + ok = OK[r, tgt] & (K[r, tgt] > 0) + else: + ok = (suit[r, tgt] > 0) & (fitness(env, st.O[r][:, go], race, cells=tgt) >= f_min) + ok = ok & wet & (P[tgt] < 1) & (tgt != go) + size = np.maximum(race.founder, 0.02 * P[go]) + move(st, r, go[ok], tgt[ok], size[ok]) diff --git a/worldhistory/nudges.py b/worldhistory/nudges.py new file mode 100644 index 0000000..4e7e352 --- /dev/null +++ b/worldhistory/nudges.py @@ -0,0 +1,19 @@ +"""Light-touch nudges: weak multipliers on growth, capacity, expansion and adaptation by race, region and years.""" +import numpy as np + + +def nudge_multipliers(nudges, races, regions, t, n): + R = len(races) + g, c, e, a = np.ones((R, n)), np.ones((R, n)), np.ones((R, n)), np.ones((R, n)) + ids = {r.id: i for i, r in enumerate(races)} + for nd in nudges: + y0, y1 = nd["years"] + if not (y0 <= t < y1): + continue + m = regions(nd["region"]) + i = ids[nd["race"]] + g[i, m] *= nd["growth"] + c[i, m] *= nd["capacity"] + e[i, m] *= nd["expansion"] + a[i, m] *= nd["adapt"] + return g, c, e, a diff --git a/worldhistory/predicates.py b/worldhistory/predicates.py new file mode 100644 index 0000000..89686f4 --- /dev/null +++ b/worldhistory/predicates.py @@ -0,0 +1,97 @@ +"""Named habitat predicates (0..1 per cell) that race configs combine. + +Syntax: "name" or "name:arg:arg"; an arg "3,4" is a list of ints. A list of expressions is their product. +Holdridge groupings follow worldgen's 38 life zones (Holdridge 1947).""" +import inspect + +import numpy as np + +HOLDRIDGE = { + "forest": [7, 8, 9, 13, 14, 15, 19, 20, 21, 22, 26, 27, 28, 29, 33, 34, 35, 36, 37], + "steppe": [6, 11, 12, 17, 18, 24, 25, 31, 32], + "desert": [0, 5, 10, 16, 23, 30], + "tundra": [1, 2, 3, 4], +} + + +def depth(F): + """Water depth in m (0 on land).""" + return np.where(F["ocean"], np.maximum(-np.asarray(F["elevation_m"], float), 0.0), 0.0) + + +def _land(w, F): return ~np.asarray(F["ocean"], bool) +def _sea(w, F): return np.asarray(F["ocean"], bool) +def _ids(a): return a if isinstance(a, list) else [int(a)] + + +def p_land(w, F): return _land(w, F) +def p_sea(w, F): return _sea(w, F) +def p_forest(w, F): return np.isin(F["holdridge"], HOLDRIDGE["forest"]) & _land(w, F) +def p_steppe(w, F): return np.isin(F["holdridge"], HOLDRIDGE["steppe"]) & _land(w, F) +def p_desert(w, F): return np.isin(F["holdridge"], HOLDRIDGE["desert"]) & _land(w, F) +def p_tundra(w, F): return np.isin(F["holdridge"], HOLDRIDGE["tundra"]) & _land(w, F) +def p_water(w, F): return np.clip(F["river"] * (0.4 + 0.15 * F["strahler"]) + F["lake"] * 0.6, 0, 1) +def p_coast(w, F): return _land(w, F) & w.nb_any(_sea(w, F)) +def p_mountain(w, F): return np.isin(F["landform"], [3, 8]) & _land(w, F) +def p_hills(w, F): return (F["landform"] == 2) & _land(w, F) +def p_karst(w, F): return F["ground"] == 8 +def p_wetland(w, F): return np.isin(F["ground"], [1, 2]) +def p_ice(w, F): return F["ice"] > 0 +def p_ice_free(w, F): return F["ice"] == 0 +def p_underground(w, F): return (p_mountain(w, F) | p_hills(w, F) | p_karst(w, F)) & _land(w, F) +def p_coal(w, F): return np.asarray(F["coal_potential"], float) > 0 +def p_iron(w, F): return np.asarray(F["iron_potential"], float) > 0 +def p_vent(w, F): return np.clip(F["vent_potential"], 0, 1) * _sea(w, F) +def p_productivity(w, F): return np.clip(np.asarray(F.get("productivity", np.zeros(len(F["ocean"]))), float), 0, 1) +def p_shelf(w, F): return _sea(w, F) & (depth(F) < 200) +def p_deep_sea(w, F): return _sea(w, F) & (depth(F) >= 1000) +def p_lithology(w, F, ids): return np.isin(F["lithology"], _ids(ids)) +def p_holdridge(w, F, ids): return np.isin(F["holdridge"], _ids(ids)) +def p_landform(w, F, ids): return np.isin(F["landform"], _ids(ids)) +def p_seabed(w, F, ids): return np.isin(F["seabed_type"], _ids(ids)) +def p_ground(w, F, ids): return np.isin(F["ground"], _ids(ids)) +def p_warm(w, F, lo, hi): return np.clip((F["T_mean"] - lo) / (hi - lo), 0, 1) +def p_wet(w, F, lo, hi): return np.clip((F["P_ann"] - lo) / (hi - lo), 0, 1) +def p_above(w, F, f, v): return np.asarray(F[f], float) > v +def p_below(w, F, f, v): return np.asarray(F[f], float) < v +def p_field(w, F, f): return np.clip(np.asarray(F[f], float), 0, 1) + + +PRED = {k[2:]: v for k, v in globals().items() if k.startswith("p_") and callable(v)} +ARITY = {k: len(inspect.signature(v).parameters) - 2 for k, v in PRED.items()} +FIELD_ARG = {"above", "below", "field"} + + +def _arg(s): + if "," in s: + return [int(v) for v in s.split(",") if v] + try: + return float(s) + except ValueError: + return s + + +def parse(expr): + name, *args = expr.split(":") + return name, [_arg(a) for a in args] + + +def check(exprs, field_names=None): + """ValueError for an unknown predicate, a wrong number of args, or a field the world lacks.""" + for e in [exprs] if isinstance(exprs, str) else exprs: + name, args = parse(e) + if name not in PRED: + raise ValueError(f"unknown predicate {name!r} in {e!r}; known: {sorted(PRED)}") + if len(args) != ARITY[name]: + raise ValueError(f"predicate {e!r} takes {ARITY[name]} argument(s)") + if name in FIELD_ARG and field_names is not None and args[0] not in field_names: + raise ValueError(f"predicate {e!r}: the world has no field {args[0]!r}") + + +def evaluate(exprs, world, F=None): + F = world.fields if F is None else F + out = np.ones(world.n) + for e in [exprs] if isinstance(exprs, str) else exprs: + name, args = parse(e) + out = out * np.asarray(PRED[name](world, F, *args), float) + return out diff --git a/worldhistory/preview.py b/worldhistory/preview.py new file mode 100644 index 0000000..5bc015b --- /dev/null +++ b/worldhistory/preview.py @@ -0,0 +1,172 @@ +"""Preview PNGs (pillow only): race map, per-race maps, population chart, contact sheet; compare two trials.""" +import json +import tomllib +from pathlib import Path + +import numpy as np +from PIL import Image, ImageDraw + +from .snapshot import load_snapshot + +SEA, LAND, INK = (12, 22, 40), (60, 55, 45), (235, 235, 235) + + +def _xy(lat, lon, W, H): + return (((lon + 180) / 360 * (W - 1)).astype(int), ((90 - lat) / 180 * (H - 1)).astype(int)) + + +STAMP = [(dx, dy) for dx in (-1, 0, 1, 2) for dy in (-1, 0, 1)] # ~ one r4 cell at 1440 px (as the spike) + + +def _base(lat, lon, ocean, W, H): + x, y = _xy(lat, lon, W, H) + arr = np.zeros((H, W, 3), np.uint8) + arr[:] = SEA + col = np.where(np.asarray(ocean, bool)[:, None], [SEA], [LAND]).astype(np.uint8) + for dx, dy in STAMP: + arr[np.clip(y + dy, 0, H - 1), np.clip(x + dx, 0, W - 1)] = col + return arr, x, y + + +def _paint(arr, x, y, cells, cols, alpha, W, H): + for dx, dy in STAMP: + yy, xx = np.clip(y[cells] + dy, 0, H - 1), np.clip(x[cells] + dx, 0, W - 1) + arr[yy, xx] = (arr[yy, xx] * (1 - alpha[:, None]) + cols * alpha[:, None]).astype(np.uint8) + + +def race_map(lat, lon, ocean, area, P, colours, path, title, W=1440, H=720): + arr, x, y = _base(lat, lon, ocean, W, H) + tot = P.sum(0) + cells = np.flatnonzero(tot >= 1) + dens = tot[cells] / area[cells] + alpha = np.clip(np.log1p(dens) / np.log1p(40), 0.15, 1) + cols = np.asarray(colours, float)[P[:, cells].argmax(0)] + _paint(arr, x, y, cells, cols, alpha, W, H) + img = Image.fromarray(arr) + d = ImageDraw.Draw(img) + d.text((8, 6), title, fill=INK) + img.save(path) + + +def small_multiples(lat, lon, ocean, area, P, meta, path, cols=3, W=480, H=240): + R = len(meta) + rows = (R + cols - 1) // cols + sheet = Image.new("RGB", (cols * W, rows * H), SEA) + for i, m in enumerate(meta): + arr, x, y = _base(lat, lon, ocean, W, H) + cells = np.flatnonzero(P[i] >= 1) + alpha = np.clip(np.log1p(P[i, cells] / area[cells]) / np.log1p(40), 0.15, 1) + _paint(arr, x, y, cells, np.asarray(m["colour"], float), alpha, W, H) + img = Image.fromarray(arr) + ImageDraw.Draw(img).text((6, 4), f"{m['id']} {P[i].sum():,.0f}", fill=INK) + sheet.paste(img, ((i % cols) * W, (i // cols) * H)) + sheet.save(path) + + +def chart(stats, meta, path, W=1000, H=500): + img = Image.new("RGB", (W, H), (250, 250, 248)) + d = ImageDraw.Draw(img) + steps = stats["steps"] + years = np.array([s["year"] for s in steps], float) + L, Rm, T, B = 70, 150, 20, 40 + top = max(max(s["pop"].values()) for s in steps) or 1.0 + ymax = np.ceil(np.log10(top + 1)) + X = lambda yr: L + (yr - years.min()) / max(years.max() - years.min(), 1) * (W - L - Rm) + Y = lambda p: H - B - np.log10(p + 1) / ymax * (H - T - B) + d.line([(L, T), (L, H - B), (W - Rm, H - B)], fill=(60, 60, 60)) + for k in range(int(ymax) + 1): + d.text((8, Y(10 ** k) - 6), f"1e{k}", fill=(90, 90, 90)) + for tick in np.linspace(years.min(), years.max(), 6): + d.text((X(tick) - 14, H - B + 6), f"{tick:.0f}", fill=(90, 90, 90)) + for i, m in enumerate(meta): + pts = [(X(s["year"]), Y(s["pop"][m["id"]])) for s in steps] + d.line(pts, fill=tuple(m["colour"]), width=2) + d.text((W - Rm + 8, T + 14 * i), m["id"], fill=tuple(m["colour"])) + img.save(path) + + +def contact_sheet(paths, path, cols=3): + ims = [Image.open(p) for p in paths] + w, h = ims[0].size + w2, h2 = w // 2, h // 2 + rows = (len(ims) + cols - 1) // cols + sheet = Image.new("RGB", (cols * w2, rows * h2), SEA) + for i, im in enumerate(ims): + sheet.paste(im.resize((w2, h2)), ((i % cols) * w2, (i // cols) * h2)) + sheet.save(path) + + +def _run(trial_dir): + run = tomllib.loads((Path(trial_dir) / "run.toml").read_text()) + meta = [{"id": r["id"], "colour": r["colour"], "family": r["family"]} for r in run["race"]] + return run, meta + + +def _world(run, world): + if world is None: + from .world import load_world + world = load_world(run["world"], eras=run["history"].get("eras", [])) + switches = sorted((e["year"], e["era"]) for e in run["history"].get("events", []) if e["kind"] == "era_switch") + + def ocean(year): + era = None + for y, e in switches: + if y <= year: + era = e + return world.base["ocean"] if era is None else world.eras[era]["ocean"] + return world, ocean + + +CONTACT_YEARS = (0, 20, 30, 200, 500, 1000, 2500, 5000, 7500) + + +def preview(trial_dir, world=None): + trial_dir = Path(trial_dir) + run, meta = _run(trial_dir) + world, ocean = _world(run, world) + out = trial_dir / "preview" + out.mkdir(exist_ok=True) + files, maps = [], {} + for snap in sorted((trial_dir / "snap").glob("y*.npz")): + s = load_snapshot(snap, world.n) + p = out / f"race_{snap.stem}.png" + race_map(world.lat, world.lon, ocean(s["year"]), world.area, s["P"], [m["colour"] for m in meta], p, + f"year {s['year']}") + files.append(p) + maps[s["year"]] = p + last = s + p = out / "races_final.png" + small_multiples(world.lat, world.lon, ocean(last["year"]), world.area, last["P"], meta, p) + files.append(p) + stats = json.loads((trial_dir / "stats.json").read_text()) + p = out / "chart.png" + chart(stats, meta, p) + files.append(p) + years = sorted(maps) + pick = [maps[min(years, key=lambda y: abs(y - c))] for c in CONTACT_YEARS if c <= years[-1]] + pick = list(dict.fromkeys(pick)) + p = out / "contact.png" + contact_sheet(pick, p) + files.append(p) + return files + + +def compare(dir_a, dir_b, world=None, out=None): + """Seed agreement: share of occupied land cells with the same dominant race in both final snapshots.""" + run_a, meta = _run(dir_a) + world, ocean = _world(run_a, world) + last = lambda d: sorted((Path(d) / "snap").glob("y*.npz"))[-1] + a, b = load_snapshot(last(dir_a), world.n), load_snapshot(last(dir_b), world.n) + occ = (a["P"].sum(0) >= 1) | (b["P"].sum(0) >= 1) + agree = float((a["P"].argmax(0)[occ] == b["P"].argmax(0)[occ]).mean()) if occ.any() else 1.0 + res = {"agreement": agree, "year": a["year"], + "pop": {m["id"]: [float(a["P"][i].sum()), float(b["P"][i].sum())] for i, m in enumerate(meta)}} + if out is not None: + out = Path(out) + out.mkdir(parents=True, exist_ok=True) + cols = [m["colour"] for m in meta] + race_map(world.lat, world.lon, ocean(a["year"]), world.area, a["P"], cols, out / "a.png", f"A {dir_a}") + race_map(world.lat, world.lon, ocean(b["year"]), world.area, b["P"], cols, out / "b.png", f"B {dir_b}") + contact_sheet([out / "a.png", out / "b.png"], out / "compare.png", cols=2) + (out / "compare.json").write_text(json.dumps(res, indent=1)) + return res diff --git a/worldhistory/regions.py b/worldhistory/regions.py new file mode 100644 index 0000000..e46da55 --- /dev/null +++ b/worldhistory/regions.py @@ -0,0 +1,54 @@ +"""Named regions: masks over cells used by events, hazards, nudges, sub-species and stats.""" +import numpy as np + +from .predicates import evaluate + +KINDS = ("all", "box", "changed", "plate", "predicate") + + +def region_mask(world, spec): + kind = spec.get("kind", "all") + n = world.n + if kind == "all": + m = np.ones(n, bool) + elif kind == "box": + (a, b), (c, d) = spec["lat"], spec["lon"] + m = (world.lat >= a) & (world.lat <= b) & (world.lon >= c) & (world.lon <= d) + elif kind == "changed": # what an era changed relative to the base world + era = world.eras[spec["era"]] + m = np.zeros(n, bool) + for f in spec.get("fields", ["gravity_g"]): + m |= np.abs(np.asarray(era[f], float) - np.asarray(world.base[f], float)) > spec.get("eps", 1e-3) + elif kind == "plate": + m = np.isin(world.base["plate"], spec["ids"]) + elif kind == "predicate": + m = evaluate(spec["p"], world, world.base) > 0.5 + else: + raise ValueError(f"unknown region kind {kind!r}") + if kind != "predicate" and spec.get("p"): # any kind can be narrowed by a predicate + m = m & (evaluate(spec["p"], world, world.base) > 0.5) + land = spec.get("land") + if land == "base": + m = m & ~np.asarray(world.base["ocean"], bool) + elif land == "current": + m = m & ~np.asarray(world.fields["ocean"], bool) + return m + + +class Regions: + def __init__(self, world, specs): + self.world, self.specs, self._cache = world, dict(specs), {} + + def __call__(self, name): + if name in (None, "all"): + return np.ones(self.world.n, bool) + spec = self.specs[name] + if spec.get("land") == "current": + m = region_mask(self.world, spec) + return m & ~self(spec["exclude"]) if spec.get("exclude") else m + if name not in self._cache: + m = region_mask(self.world, spec) + if spec.get("exclude"): + m = m & ~self(spec["exclude"]) + self._cache[name] = m + return self._cache[name] diff --git a/worldhistory/snapshot.py b/worldhistory/snapshot.py new file mode 100644 index 0000000..19d83de --- /dev/null +++ b/worldhistory/snapshot.py @@ -0,0 +1,34 @@ +"""Output: sparse snapshots (npz), stats (json), engine version.""" +import json +import subprocess +from pathlib import Path + +import numpy as np + + +def write_snapshot(path, st, t, races): + r, c = np.nonzero(st.P >= 0.5) + np.savez_compressed(path, year=t, races=np.array([x.id for x in races]), race=r.astype(np.int16), + cell=c.astype(np.int32), pop=st.P[r, c].astype(np.float32), + opt=st.O[r, :, c].astype(np.float32), tech=st.T[r, :, c].astype(np.float32)) + + +def load_snapshot(path, n): + z = np.load(path) + races = [str(x) for x in z["races"]] + P = np.zeros((len(races), n), np.float32) + P[z["race"], z["cell"]] = z["pop"] + return {"year": int(z["year"]), "races": races, "P": P, "race": z["race"], "cell": z["cell"], + "opt": z["opt"], "tech": z["tech"]} + + +def write_json(path, obj): + Path(path).write_text(json.dumps(obj, indent=1)) + + +def engine_version(): + try: + return subprocess.run(["git", "-C", str(Path(__file__).parent), "rev-parse", "--short", "HEAD"], + capture_output=True, text=True, check=True).stdout.strip() + except (OSError, subprocess.CalledProcessError): + return "unknown" diff --git a/worldhistory/state.py b/worldhistory/state.py new file mode 100644 index 0000000..d6211f2 --- /dev/null +++ b/worldhistory/state.py @@ -0,0 +1,183 @@ +"""Simulation state and the two ways people change place: move (cell to cell) and convert (race slot to slot). +Per-cell attributes (adapted optima O, tech T, curse C) travel with people as population-weighted mixes.""" +import math +from dataclasses import dataclass + +import numpy as np + +from .config import DOMAINS + +ATTRS = ("O", "T", "C") + + +@dataclass +class State: + P: np.ndarray # (R, n) people + O: np.ndarray # (R, C, n) adapted optimum per condition + T: np.ndarray # (R, D, n) tech level per domain, 0..1 + C: np.ndarray # (R, 1, n) curse: per-step death share carried by perpetrators + improve: np.ndarray # (R, n) land improvement, 0..land_cap + push: np.ndarray # (R, n) extra emigrants this step (yielding from conflict) + rng: np.random.Generator + t: int = 0 + + +def new_state(R, n, native, seed): + from .config import CONDITIONS + native = np.asarray(native, np.float32) + if native.shape != (R, len(CONDITIONS)): + raise ValueError(f"native optima {native.shape}: need ({R}, {len(CONDITIONS)}) — one per race × conditions") + O = np.repeat(native[:, :, None], n, axis=2) + return State(P=np.zeros((R, n)), O=O, + T=np.zeros((R, len(DOMAINS), n), np.float32), C=np.zeros((R, 1, n)), improve=np.zeros((R, n), np.float32), + push=np.zeros((R, n)), rng=np.random.default_rng(seed)) + + +def _change_cells_np(P, wide, src, dst): + mask = np.zeros(len(P), bool) + mask[src] = True + mask[dst] = True + mask |= np.signbit(P) | ((P > 0) & (P < 1e-12)) + for A in wide: + mask |= (A != 0).any(0) & (P > 0) + return np.flatnonzero(mask) + + +def _make_change_cells_jit(): + """_change_cells_np in one compiled pass (numba optional; WORLDHISTORY_NO_JIT=1 turns it off).""" + import os + if os.environ.get("WORLDHISTORY_NO_JIT"): + return None + try: + import numba + except ImportError: + return None + + @numba.njit(cache=True) + def cells(P, A, src, dst): + n = len(P) + mask = np.zeros(n, np.bool_) + for s in src: + mask[s] = True + for d in dst: + mask[d] = True + k = 0 + for i in range(n): + p = P[i] + if not mask[i]: + if math.copysign(1.0, p) < 0 or (p > 0 and p < 1e-12): + mask[i] = True + elif p > 0: + for a in range(A.shape[0]): + if A[a, i] != 0: + mask[i] = True + break + if mask[i]: + k += 1 + out = np.empty(k, np.int64) + k = 0 + for i in range(n): + if mask[i]: + out[k] = i + k += 1 + return out + return cells + + +_change_cells_jit = _make_change_cells_jit() + + +def _make_mix_jit(): + """move's attribute mix, compiled (numba optional, as change_cells): the inflow sums in move order from 0 (as + bincount adds), then (stay·A + S) / tot where tot > 0, cast back to A's dtype — the same floats.""" + import os + if os.environ.get("WORLDHISTORY_NO_JIT"): + return None + try: + import numba + except ImportError: + return None + + @numba.njit(cache=True) + def mix(A, cs, sl, dl, amt, stay, tot): + m = len(cs) + S = np.empty(m) + for x in range(A.shape[0]): + S[:] = 0.0 + for i in range(len(dl)): + S[dl[i]] += amt[i] * np.float64(A[x, cs[sl[i]]]) + for c in range(m): + if tot[c] > 0: + A[x, cs[c]] = (stay[c] * np.float64(A[x, cs[c]]) + S[c]) / max(tot[c], 1e-12) + return mix + + +_mix_jit = _make_mix_jit() + + +def change_cells(P, wide, src, dst): + """Cells a move can change: its ends, plus every cell where the full-array update would not be a no-op + (P < 0 or −0, 0 < P < 1e-12, nonzero float64 attributes where P > 0; float32 ones round back exactly).""" + if _change_cells_jit is not None and len(wide) <= 1 and P.dtype == np.float64 and P.flags.c_contiguous: + A = wide[0] if wide else np.zeros((0, len(P))) + if A.dtype == np.float64: + return _change_cells_jit(P, np.ascontiguousarray(A), src, dst) + return _change_cells_np(P, wide, src, dst) + + +def move(st, r, src, dst, amt): + """Move amt people of race r from src to dst cells (parallel arrays); sources give at most what they have. + Works on the cells it can change only: the moves' ends, plus every cell where the full-array update would not + be a no-op (P < 0 or −0, 0 < P < 1e-12, nonzero float64 attributes; float32 ones round back exactly). Same + results, bit for bit, as updating every cell.""" + src = np.asarray(src, np.int64) + dst = np.asarray(dst, np.int64) + amt = 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] + wide = [getattr(st, name)[r] for name in ATTRS if getattr(st, name)[r].shape[0] and getattr(st, name).dtype != np.float32] + cs = change_cells(P, wide, src, dst) + m = len(cs) + loc = np.empty(n, np.int64) + loc[cs] = np.arange(m) + sl, dl = loc[src], loc[dst] + Pc = P[cs] + out = np.bincount(sl, amt, m) + scale = np.where(out > Pc, Pc / np.maximum(out, 1e-12), 1.0) + amt = amt * scale[sl] + stay = np.maximum(Pc - np.bincount(sl, amt, m), 0.0) + inflow = np.bincount(dl, amt, m) + tot = stay + inflow + for name in ATTRS: + A = getattr(st, name)[r] + if A.shape[0] == 0: + continue + if _mix_jit is not None and A.flags.c_contiguous: + _mix_jit(A, cs, sl, dl, amt, stay, tot) + continue + Ac = A[:, cs] + S = np.stack([np.bincount(dl, amt * Ac[x, sl], m) for x in range(A.shape[0])]) + A[:, cs] = np.where(tot > 0, (stay * Ac + S) / np.maximum(tot, 1e-12), Ac) + P[cs] = tot + + +def convert(st, r_from, r_to, cells, amt, O_new=None): + """Turn amt people of slot r_from into slot r_to in the same (unique) cells, carrying their attributes; with + O_new (C, len(cells)) they arrive with those optima instead (a birth resets the curves).""" + cells = np.asarray(cells, np.int64) + amt = np.minimum(np.asarray(amt, float), st.P[r_from, cells]) + p_to = st.P[r_to, cells] + tot = p_to + amt + for name in ATTRS: + A = getattr(st, name) + if A.shape[1] == 0: + continue + a_to = A[r_to][:, cells] + a_from = O_new if (name == "O" and O_new is not None) else A[r_from][:, cells] + A[r_to][:, cells] = np.where(tot > 0, (p_to * a_to + amt * a_from) / np.maximum(tot, 1e-12), a_to) + st.P[r_from, cells] -= amt + st.P[r_to, cells] += amt diff --git a/worldhistory/tech.py b/worldhistory/tech.py new file mode 100644 index 0000000..13ec20a --- /dev/null +++ b/worldhistory/tech.py @@ -0,0 +1,78 @@ +"""Tech domains (incl. magic): levels 0..1 per race per cell. Progress grows with the connected population +(Kremer 1993) and small isolated groups lose tech (Henrich 2004). Effects: capacity, shocks, travel, reach, land.""" +import numpy as np + +from .config import DOMAINS + +DI = {d: i for i, d in enumerate(DOMAINS)} + + +def effective(T_r, tp): + m = T_r[DI["magic"]] + fold = lambda d: np.clip(T_r[DI[d]] + tp["magic_share"] * m, 0, 1) + return {"farming": fold("farming"), "crafts": fold("crafts"), "travel": fold("travel"), "land": fold("land"), + "magic": m, "war": T_r[DI["war"]]} + + +def k_gain(eff, tp): + return tp["farming_gain"] * eff["farming"] + tp["crafts_gain"] * eff["crafts"] + + +def shock_scale(eff, tp): + return 1 - tp["crafts_shock"] * eff["crafts"] + + +def travel_bonus(eff, tp): + return tp["travel_cost"] * eff["travel"] + + +def jump_boost(eff, tp): + return 1 + tp["travel_jump"] * eff["travel"] + + +def reach(eff, race, tp): + return np.minimum(race.reach_cap, tp["crafts_reach"] * eff["crafts"] + tp["magic_reach"] * eff["magic"]) + + +def neigh_pop(world, P_r, passes): + x = P_r + for _ in range(passes): + x = x + world.nb_sum(x) + return x + + +def step_tech(st, world, races, tp, steps=1.0): + for r, race in enumerate(races): + P = st.P[r] + occ = P >= 1 + if not occ.any(): + continue + oc = np.flatnonzero(occ) # only occupied cells change: compute there + loc = world.nb_local(oc) # their neighbours: all the diffusion below reads + if tp["passes"] == 1: + N = P[oc] + world.nb_apply(loc[0], P[loc[1]]) + else: + 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: + sub, cols = loc + Pc, Tc = P[oc], T[:, oc] + Pn, Tn = P[cols], T[:, cols] + M = (Pc * Tc + world.nb_apply(sub, Pn * Tn)) / np.maximum(Pc + world.nb_apply(sub, Pn), 1e-12) + T[:, oc] = Tc + tp["diffuse"] * np.clip(M - Tc, 0, None) + + +def step_land(st, races, tp, steps=1.0): + for r in range(len(races)): + occ = st.P[r] >= 1 + L = effective(st.T[r], tp)["land"] + imp = st.improve[r] + up = imp + tp["land_rate"] * steps * L * (tp["land_cap"] - imp) + down = imp * (1 - tp["land_decay"] * steps) + st.improve[r] = np.clip(np.where(occ, up, down), 0, tp["land_cap"]) diff --git a/worldhistory/world.py b/worldhistory/world.py new file mode 100644 index 0000000..f21bcbe --- /dev/null +++ b/worldhistory/world.py @@ -0,0 +1,190 @@ +"""The world a history runs on: cells, neighbours and fields (base world + eras). Game-agnostic.""" +from __future__ import annotations + +import json +from dataclasses import dataclass, field +from functools import cached_property +from pathlib import Path + +import numpy as np +from scipy.spatial import cKDTree + +# worldgen cells.npz fields the engine reads (more via extra_fields) +FIELDS = ("ocean", "elevation_m", "T_mean", "P_ann", "holdridge", "landform", "river", "strahler", "lake", "ice", + "coal_potential", "iron_potential", "vent_potential", "seabed_type", "lithology", "ground", "gravity_g", + "o2_fraction", "pressure_bar", "bottom_temp_c", "plate") + +OPTIONAL = ("current", "productivity", "upwelling", "sst") # ocean fields of newer worldgen builds; absent → still water + + +@dataclass +class World: + nb: np.ndarray # (n, 6) neighbour indices, -1 = none (pentagons, line worlds) + lat: np.ndarray # degrees + lon: np.ndarray + xyz: np.ndarray # (n, 3) unit vectors + area: np.ndarray # km² + radius_km: float + base: dict # field name -> (n,) array, the world before any era + eras: dict = field(default_factory=dict) + era: str | None = None + meta: dict = field(default_factory=dict) + + @property + def n(self) -> int: + return len(self.lat) + + @property + def fields(self) -> dict: + return self.base if self.era is None else self.eras[self.era] + + def set_era(self, name): + if name is not None and name not in self.eras: + raise KeyError(f"era {name!r} not loaded") + self.era = name + self.__dict__.pop("_sail", None) + + @cached_property + def tree(self): + return cKDTree(self.xyz) + + @cached_property + def _nbj(self): + return np.clip(self.nb, 0, None) + + @cached_property + def _nbv(self): + return self.nb >= 0 + + @cached_property + def _adj(self): + from scipy.sparse import csr_matrix + i = np.repeat(np.arange(self.n), self.nb.shape[1])[self._nbv.ravel()] + j = self.nb.ravel()[self._nbv.ravel()] + return csr_matrix((np.ones(len(i)), (i, j)), shape=(self.n, self.n)) + + @property + def current(self): + c = self.fields.get("current") + return np.zeros((self.n, 3)) if c is None else np.asarray(c, dtype=np.float64) + + @cached_property + def _nbt(self): + """(n, 6, 3) unit tangents from each cell toward each neighbour (0 where there is none).""" + p, q = self.xyz[:, None, :], self.xyz[self._nbj] + t = q - np.sum(p * q, axis=-1, keepdims=True) * p + t = t / np.maximum(np.linalg.norm(t, axis=-1, keepdims=True), 1e-15) + return np.where(self._nbv[..., None], t, 0.0) + + def sail_cost(self, ship_ms): + """Hours to sail each directed edge (n, 6) at ship_ms through still water, helped or slowed by the mean + current of its two cells; inf where there is no edge or either end is land.""" + cache = self.__dict__.setdefault("_sail", {}) + if ship_ms not in cache: + sea = np.asarray(self.fields["ocean"], bool) + cur = self.current + ubar = 0.5 * (cur[:, None, :] + cur[self._nbj]) + along = np.sum(ubar * self._nbt, axis=-1) + d_m = self.radius_km * 1e3 * np.arccos(np.clip(np.sum(self.xyz[:, None, :] * self.xyz[self._nbj], -1), -1, 1)) + v = np.maximum(ship_ms + along, 1e-3) + ok = self._nbv & sea[:, None] & sea[self._nbj] + cache[ship_ms] = np.where(ok, d_m / v / 3600.0, np.inf) + return cache[ship_ms] + + def nb_sum(self, x): + """Sum over each cell's neighbours; x has cells on its last axis.""" + x = np.asarray(x) + if x.ndim == 1: + return self._adj @ x + flat = x.reshape(-1, self.n) + return (self._adj @ flat.T).T.reshape(x.shape) + + def nb_local(self, rows): + """(sub, cols): the adjacency rows `rows` over only the columns they touch (sorted `cols`), so + sub @ x[..., cols] == nb_sum(x)[..., rows] with the same terms in the same order (the same floats), + and the caller need only compute x at cols.""" + from scipy.sparse import csr_matrix + adj = self._adj[rows] + mark = np.zeros(self.n, bool) + mark[adj.indices] = True + cols = np.flatnonzero(mark) + inv = np.empty(self.n, adj.indices.dtype) + inv[cols] = np.arange(len(cols), dtype=adj.indices.dtype) + sub = csr_matrix((adj.data, inv[adj.indices], adj.indptr), shape=(len(rows), len(cols))) + return sub, cols + + @staticmethod + def nb_apply(sub, xc): + """sub @ xc over the last axis (xc: values at nb_local's cols).""" + xc = np.asarray(xc) + if xc.ndim == 1: + return sub @ xc + flat = xc.reshape(-1, xc.shape[-1]) + return (sub @ flat.T).T.reshape(xc.shape[:-1] + (sub.shape[0],)) + + def nb_sum_at(self, x, rows, local=None): + """nb_sum(x)[..., rows], summing only those rows (same per-row order: the same floats).""" + sub, cols = self.nb_local(rows) if local is None else local + return self.nb_apply(sub, np.asarray(x)[..., cols]) + + def nb_any(self, mask): + return (self._nbv & mask[self._nbj]).any(1) + + def smooth(self, x, passes=1): + cnt = 1 + self._nbv.sum(1) + for _ in range(passes): + x = (x + self.nb_sum(x)) / cnt + return x + + def noise(self, seed, passes=3): + """Smooth random field, mean 0, std 1, fixed by the seed.""" + z = self.smooth(np.random.default_rng(seed).normal(size=self.n), passes) + return (z - z.mean()) / (z.std() + 1e-12) + + def km(self, a, b): + d = np.clip((self.xyz[a] * self.xyz[b]).sum(-1), -1, 1) + return np.arccos(d) * self.radius_km + + def within(self, cell, km): + chord = 2 * np.sin(km / self.radius_km / 2) + return np.asarray(self.tree.query_ball_point(self.xyz[cell], chord + 1e-12), dtype=np.int64) + + def nearest(self, xyz): + return self.tree.query(xyz)[1] + + +def neighbours_from_ids(ids): + import h3.api.basic_int as h3 + idx = {int(c): i for i, c in enumerate(ids)} + nb = np.full((len(ids), 6), -1, np.int64) + for i, c in enumerate(ids): + k = 0 + for d in h3.grid_disk(int(c), 1): + j = idx.get(int(d)) + if j is not None and j != i: + nb[i, k] = j + k += 1 + return nb + + +def _fields(z, extra): + return {f: np.asarray(z[f]) for f in (*FIELDS, *OPTIONAL, *extra) if f in z.files} + + +def load_world(build_dir, eras=(), extra_fields=()): + """A worldgen build dir (out/r<res>): cells.npz (base world) + eras/<era>/cells.npz.""" + d = Path(build_dir) + z = np.load(d / "cells.npz") + meta_p = d / "cells_meta.json" + meta = json.loads(meta_p.read_text()) if meta_p.exists() else {} + xyz = np.asarray(z["g_xyz"], float) + xyz = xyz / np.linalg.norm(xyz, axis=1, keepdims=True) + w = World(nb=neighbours_from_ids(z["g_ids"]), lat=np.asarray(z["g_lat"], float), lon=np.asarray(z["g_lon"], float), + xyz=xyz, area=np.asarray(z["g_area_km2"], float), radius_km=float(meta.get("radius_km", 6371.0)), + base=_fields(z, extra_fields), meta={**meta, "dir": str(d)}) + missing = [f for f in FIELDS if f not in w.base] + if missing: + raise ValueError(f"{d / 'cells.npz'} lacks fields {missing}") + for e in eras: + w.eras[e] = _fields(np.load(d / "eras" / e / "cells.npz"), extra_fields) + return w |
