diff --git a/CHANGELOG.md b/CHANGELOG.md index b7a2b2e..dba357b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,7 @@ ## [Unreleased] -Four models now follow the TerraME `ca` package (https://github.com/terrame/ca), +Six models now follow the TerraME `ca` package (https://github.com/terrame/ca), where they used to differ. The results of these models change; a new minor release (0.2.0) is suggested. @@ -18,6 +18,27 @@ release (0.2.0) is suggested. `initial_fire_density` now defaults to 0.0 and, above zero, starts a random fraction of the cells burning instead (with `empty=0` and a density of 0.05 the old setup is recovered, though not the same random draw). +- `GameOfLife` (vector and raster): the Moore neighborhood wraps around the + grid edges, as in TerraME's `Life` (`wrap = true`); `wrap=False` keeps the + bounded grid. A glider now crosses the edge instead of dying on it. +- `FireModel` (raster): the fire spreads to the four Von Neumann neighbors, + as in TerraME and in the vector version. The rule passed no neighborhood + to `focal_sum_mask`, which counts all eight (Moore), so the fire also + spread diagonally. +- `SolidDiffusion`: the step follows TerraME. Each vacancy draws one of all + its neighbors (a vacancy may draw another vacancy and stay put) and moves + at most once per step, reading the state at the start of the step. The + port drew only among atoms and could move a vacancy several times in one + step, so the atoms mixed faster than in TerraME (after 400 steps, 35.9 + atoms of one kind past the middle column against 26.5 in TerraME; now + 27). +- `tests/test_goldens.py` (marker `golden`): the models against TerraME 2.0.1 + runs from `terrame-goldens` v0.1.0. Deterministic models start from + TerraME's state and must match it cell by cell; stochastic ones are + compared with 50 TerraME runs by a Kolmogorov-Smirnov test and a test on + the means. About a minute by default; `DISSMODEL_GOLDEN_FULL=1` compares + every step and adds the slow ensembles. `tests/goldens/` holds the + extract, its SHA-256 provenance and the script that makes it. - New `dissmodel_ca.models._neighborhood.grid_neighbors`: Von Neumann or Moore neighborhoods built from the `"row-col"` index, with optional wrap. - `tests/test_terrame_rules.py`: Parity and Excitable are compared step by step diff --git a/README.md b/README.md index 124e991..78314cd 100644 --- a/README.md +++ b/README.md @@ -67,3 +67,21 @@ jupyter lab examples/notebooks/ca_game_of_life.ipynb | `Excitable` | Vector | Excitable medium waves (spiral/ring patterns). | | `Parasit` | Vector | Host-parasit spatial dynamics. | | `Interspecific` | Vector | Grass species competition model. | + +--- + +## ✅ Verification against TerraME + +The models were checked against TerraME 2.0.1 running the `ca` package +([`terrame-goldens`](https://github.com/LambdaGeo/terrame-goldens) v0.1.0): + +- **Deterministic rules** (`GameOfLife` vector and raster, `Parity`, + `Excitable`, `Wolfram` 30/90/110, `FireModel` vector and raster, `Anneal`, + `Oscillator`, `Parasit`): started from TerraME's state, they match it + **cell by cell at every step** (up to 500 steps). +- **Stochastic models** (`FireModel` initial state, `SolidDiffusion`, + `Growth`, `Snow`, `InterspecificCompetition`): compared with 50 TerraME + runs by distribution. + +`tests/test_goldens.py` runs these comparisons (`pytest -m golden`, about a +minute; `DISSMODEL_GOLDEN_FULL=1` for every step and the slow ensembles). diff --git a/docs/reference.md b/docs/reference.md index 0d59869..64428f5 100644 --- a/docs/reference.md +++ b/docs/reference.md @@ -196,7 +196,7 @@ initial setup (10% empty cells and a single ignition by default). ```python def rule(self, arrays: dict) -> dict: state = arrays[self.state_attr] - has_burning = self.backend.focal_sum_mask(state == int(FireState.BURNING)) > 0 + has_burning = self.backend.focal_sum_mask(state == int(FireState.BURNING), self.dirs) > 0 new_state = state.copy() new_state = np.where(state == int(FireState.BURNING), int(FireState.BURNED), new_state) new_state = np.where((state == int(FireState.FOREST)) & has_burning, int(FireState.BURNING), new_state) @@ -229,7 +229,9 @@ The classic zero-player game demonstrating mathematical emergence. `0 - Dead`, `1 - Alive`. **3. Neighborhood:** -**Moore (Queen)** neighborhood (8 neighbors). +**Moore (Queen)** neighborhood (8 neighbors), wrapped around the grid edges +(torus), as in TerraME (`createNeighborhood{wrap = true}`). +`GameOfLife(gdf=grid, wrap=False)` gives a bounded grid. **4. Transition Rules:** ```python @@ -266,13 +268,16 @@ Raster-based implementation of Conway's Game of Life using NumPy. `0 - Dead`, `1 - Alive`. **3. Neighborhood:** -**Moore (Queen)** neighborhood. +**Moore (Queen)** neighborhood, wrapped around the grid edges as in TerraME; +`wrap=False` gives a bounded grid. **4. Transition Rules:** ```python def rule(self, arrays: dict) -> dict: state = arrays[self.state_attr] - neighbors = self.backend.focal_sum_mask(state == 1) + alive = (state == 1).astype(np.int8) + neighbors = sum(np.roll(alive, (dr, dc), axis=(0, 1)) # wrap = True + for dr in (-1, 0, 1) for dc in (-1, 0, 1) if (dr, dc) != (0, 0)) survive = (state == 1) & np.isin(neighbors, [2, 3]) born = (state == 0) & (neighbors == 3) return {self.state_attr: np.where(survive | born, 1, 0).astype(np.int8)} @@ -570,7 +575,11 @@ Atomic diffusion via random vacancy swapping. **Moore (Queen)**. **4. Transition Rules:** -Uses a **sequential random sweep** in `execute()`. +As in TerraME, in `execute()`: every cell that was a vacancy at the start of +the step draws one of its neighbors at random and takes its atom, unless that +neighbor was a vacancy at the start of the step or its atom has already been +taken in this step. A vacancy moves at most once per step; cells are visited +column by column, in TerraME's order. **5. Usage Example:** ```python diff --git a/pyproject.toml b/pyproject.toml index 6222cac..c95c9dd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,3 +41,6 @@ Homepage = "https://github.com/DisSModel/dissmodel-ca" [tool.pytest.ini_options] testpaths = ["tests"] +markers = [ + "golden: compares the models with TerraME runs (tests/goldens; DISSMODEL_GOLDEN_FULL=1 for every step)", +] diff --git a/src/dissmodel_ca/models/fire_model_raster.py b/src/dissmodel_ca/models/fire_model_raster.py index a84c284..2c7cac2 100644 --- a/src/dissmodel_ca/models/fire_model_raster.py +++ b/src/dissmodel_ca/models/fire_model_raster.py @@ -183,7 +183,7 @@ def rule(self, arrays: dict) -> dict: # count BURNING neighbors (Von Neumann / Rook — 4 directions) has_burning = self.backend.focal_sum_mask( - state == int(FireState.BURNING) + state == int(FireState.BURNING), self.dirs ) > 0 new_state = state.copy() diff --git a/src/dissmodel_ca/models/game_of_life.py b/src/dissmodel_ca/models/game_of_life.py index 6c3f578..00c43c6 100644 --- a/src/dissmodel_ca/models/game_of_life.py +++ b/src/dissmodel_ca/models/game_of_life.py @@ -3,11 +3,11 @@ import random from typing import Any -from libpysal.weights import Queen - from dissmodel.geo import FillStrategy, fill from dissmodel.geo import CellularAutomaton +from ._neighborhood import grid_neighbors + # --------------------------------------------------------------------------- # Built-in patterns @@ -99,6 +99,13 @@ class GameOfLife(CellularAutomaton): Extra keyword arguments forwarded to :class:`~dissmodel.geo.CellularAutomaton`. + Notes + ----- + As in the TerraME version (``createNeighborhood{wrap = true}``), the + Moore neighborhood wraps around the grid edges (torus): a glider that + leaves on one side comes back on the other. Pass ``wrap=False`` for a + bounded grid, where cells beyond the edge count as dead. + Examples -------- >>> from dissmodel.geo import vector_grid @@ -109,9 +116,20 @@ class GameOfLife(CellularAutomaton): >>> gol.initialize() """ - def setup(self) -> None: - """Build the Queen neighborhood for the grid.""" - self.create_neighborhood(strategy=Queen, use_index=True) + def setup(self, wrap: bool = True) -> None: + """ + Build the Moore (8-direction) neighborhood. + + Parameters + ---------- + wrap : bool, optional + If ``True`` (default, as in TerraME), opposite edges of the + grid are neighbors (toroidal grid). + """ + self.wrap = wrap + self.create_neighborhood( + neighbors_dict=grid_neighbors(self.gdf.index, "moore", wrap=wrap) + ) def initialize(self) -> None: """ diff --git a/src/dissmodel_ca/models/game_of_life_raster.py b/src/dissmodel_ca/models/game_of_life_raster.py index 3fc20ea..c33e80b 100644 --- a/src/dissmodel_ca/models/game_of_life_raster.py +++ b/src/dissmodel_ca/models/game_of_life_raster.py @@ -73,6 +73,10 @@ class GameOfLife(RasterCellularAutomaton): Random seed for initialization, by default 42. state_attr : str, optional Name of the state array in the backend, by default ``"state"``. + wrap : bool, optional + If ``True`` (default, as in TerraME's ``createNeighborhood{wrap = + true}``), the neighborhood wraps around the grid edges (torus). + With ``False``, cells beyond the edge count as dead. Examples -------- @@ -91,10 +95,12 @@ def setup( density: float = 0.3, seed: int = 42, state_attr: str = "state", + wrap: bool = True, ) -> None: super().setup(backend, state_attr=state_attr) self.density = density self.seed = seed + self.wrap = wrap def initialize(self) -> None: """ @@ -112,7 +118,7 @@ def rule(self, arrays: dict) -> dict: Vectorized Conway transition rule. Applied once per step over the entire grid using Moore neighborhood - (8 directions — all 8 neighbors counted via focal_sum_mask). + (8 directions), wrapped around the edges unless ``wrap=False``. Parameters ---------- @@ -124,8 +130,15 @@ def rule(self, arrays: dict) -> dict: dict[str, np.ndarray] Updated ``"state"`` array (0 = dead, 1 = alive). """ - state = arrays[self.state_attr] - neighbors = self.backend.focal_sum_mask(state == 1) + state = arrays[self.state_attr] + alive = (state == 1).astype(np.int8) + if self.wrap: + neighbors = sum( + np.roll(alive, (dr, dc), axis=(0, 1)) + for dr in (-1, 0, 1) for dc in (-1, 0, 1) if (dr, dc) != (0, 0) + ) + else: + neighbors = self.backend.focal_sum_mask(alive == 1) survive = (state == 1) & np.isin(neighbors, [2, 3]) born = (state == 0) & (neighbors == 3) diff --git a/src/dissmodel_ca/models/solid_diffusion.py b/src/dissmodel_ca/models/solid_diffusion.py index 1bc938f..7aa225e 100644 --- a/src/dissmodel_ca/models/solid_diffusion.py +++ b/src/dissmodel_ca/models/solid_diffusion.py @@ -4,6 +4,7 @@ from enum import IntEnum from typing import Any +import numpy as np from dissmodel.geo import CellularAutomaton, parse_idx @@ -25,35 +26,12 @@ class SolidDiffusionState(IntEnum): VACANCY = 2 -def _queen_neighbors(idx: str, dim: int) -> list[str]: - """ - Return the Moore (Queen) neighbor indices for a cell on a square grid. - - Computes neighbors directly from the ``"{y}-{x}"`` index format, - clipping at grid boundaries (no wrap). - - Parameters - ---------- - idx : str - Cell index in ``"{y}-{x}"`` format. - dim : int - Grid side length (number of rows = number of columns). - - Returns - ------- - list[str] - Indices of valid neighbors (up to 8). - """ - x, y = parse_idx(idx) - neighbors = [] - for dy in (-1, 0, 1): - for dx in (-1, 0, 1): - if dy == 0 and dx == 0: - continue - nx, ny = x + dx, y + dy - if 0 <= nx < dim and 0 <= ny < dim: - neighbors.append(f"{nx}-{ny}") - return neighbors +#: (d_row, d_col) offsets of the Moore neighborhood, without wrap (TerraME default). +_MOORE: tuple[tuple[int, int], ...] = ( + (-1, -1), (-1, 0), (-1, 1), + (0, -1), (0, 1), + (1, -1), (1, 0), (1, 1), +) class SolidDiffusion(CellularAutomaton): @@ -82,16 +60,19 @@ class SolidDiffusion(CellularAutomaton): Notes ----- - Because diffusion requires atomic swaps (two cells updated in a single - operation), synchronous CA update is physically incorrect here — it - would create or destroy atoms. :meth:`execute` therefore uses a - **sequential random sweep**: vacancies are processed in a shuffled - order, each swapping with one randomly chosen non-vacant neighbor. - This matches the intended stochastic diffusion dynamics. + The step follows the TerraME version (``ca`` package, ``SolidDiffusion``): + + - every cell reads the state at the start of the step (``cell.past``); + - each cell that was a vacancy draws **one** of its Moore neighbors at + random (a vacancy may draw another vacancy, and then stays put); + - it takes the drawn atom, which becomes a vacancy, unless that atom + was a vacancy at the start of the step or has already been taken by + another vacancy in this step. - Neighbor indices are computed directly from the ``"{y}-{x}"`` grid - index format via :func:`_queen_neighbors`, avoiding reliance on any - internal weights attribute. + So a vacancy moves at most once per step, and a vacancy created in the + step does not move until the next one. Cells are visited in TerraME's + order (column by column, and row by row within a column). The swaps + keep the number of atoms of each kind. Examples -------- @@ -110,7 +91,7 @@ def setup(self, seed: int = 42) -> None: Parameters ---------- seed : int, optional - Random seed for the sequential shuffle in :meth:`execute`, + Random seed for the neighbor draws in :meth:`execute`, by default 42. """ self.seed = seed @@ -144,36 +125,38 @@ def assign(idx: str) -> int: def execute(self) -> None: """ - Perform one diffusion step via a sequential random vacancy sweep. - - Each vacancy in a randomly shuffled order picks one non-vacant - neighbor at random and swaps states with it. Vacancies that have - already been filled earlier in the same sweep are skipped. - - Neighbor indices are computed directly from grid indices using - :func:`_queen_neighbors` — no spatial weights object required. + Perform one diffusion step with TerraME's rule (see the class notes). """ assert self.dim is not None, "dim must be set — pass dim=N when instantiating" - indices = list(self.gdf.index) - self._rng.shuffle(indices) - - for idx in indices: - if self.gdf.loc[idx, "state"] != int(SolidDiffusionState.VACANCY): - continue - - neighbor_ids = _queen_neighbors(idx, self.dim) - non_vacant = [ - n for n in neighbor_ids - if self.gdf.loc[n, "state"] != int(SolidDiffusionState.VACANCY) - ] - - if not non_vacant: - continue - - chosen = self._rng.choice(non_vacant) - self.gdf.loc[idx, "state"] = self.gdf.loc[chosen, "state"] - self.gdf.loc[chosen, "state"] = int(SolidDiffusionState.VACANCY) + index = list(self.gdf.index) + pos = [parse_idx(i) for i in index] + n_rows = max(p.row for p in pos) + 1 + n_cols = max(p.col for p in pos) + 1 + rows = np.array([p.row for p in pos]) + cols = np.array([p.col for p in pos]) + + past = np.empty((n_rows, n_cols), dtype=np.int64) + past[rows, cols] = self.gdf["state"].to_numpy() + cur = past.copy() + vacancy = int(SolidDiffusionState.VACANCY) + + # TerraME visits the cells column by column (x outer, y inner) + for c in range(n_cols): + for r in range(n_rows): + if past[r, c] != vacancy: + continue + neighbors = [ + (r + dr, c + dc) + for dr, dc in _MOORE + if 0 <= r + dr < n_rows and 0 <= c + dc < n_cols + ] + nr, nc = neighbors[self._rng.randrange(len(neighbors))] + if past[nr, nc] != vacancy and cur[nr, nc] != vacancy: + cur[r, c] = past[nr, nc] + cur[nr, nc] = vacancy + + self.gdf["state"] = cur[rows, cols] def rule(self, idx: Any) -> int: # pragma: no cover """Not used. Diffusion logic is in :meth:`execute`.""" diff --git a/tests/goldens/README.md b/tests/goldens/README.md new file mode 100644 index 0000000..57f3eb9 --- /dev/null +++ b/tests/goldens/README.md @@ -0,0 +1,29 @@ +# TerraME reference values for the ca models + +An extract of the `ca` suite of +[`terrame-goldens`](https://github.com/LambdaGeo/terrame-goldens) v0.1.0: +what TerraME 2.0.1 computed running the models of the +[`ca`](https://github.com/TerraME/ca) package (commit `854a5c4`) with default +parameters. + +- `grids/.csv.gz`: the state of every cell (columns `step, y, x0, x1, + ...`, TerraME's values) for the deterministic models, at every step up to + 10 and then at steps 20, 50, 100, ... and the final one. Step 0 is the + state the test starts from. `grids/cases.csv` lists model, arguments, + seed, final time and size of each case. +- `ensembles.csv.gz`: for the stochastic models, the number of cells in each + state at a few steps, for seeds 1 to 50 (columns `case, value, step, seed, + count`). For SolidDiffusion, `atom1_right` and `atom2_left` count the atoms + past the middle column. +- `reference.json`: the source version, the package commit, the steps kept + and the SHA-256 of each source file. + +To make the extract again from a copy of `terrame-goldens`: + +```bash +python tests/goldens/make_reference.py /path/to/terrame-goldens +``` + +The values are the output of TerraME and of the `ca` package (INPE and +TerraLAB/UFOP); the extract is offered under CC-BY-4.0, as the goldens it +comes from. diff --git a/tests/goldens/ensembles.csv.gz b/tests/goldens/ensembles.csv.gz new file mode 100644 index 0000000..bd5127e Binary files /dev/null and b/tests/goldens/ensembles.csv.gz differ diff --git a/tests/goldens/grids/anneal.csv.gz b/tests/goldens/grids/anneal.csv.gz new file mode 100644 index 0000000..70bc7b1 Binary files /dev/null and b/tests/goldens/grids/anneal.csv.gz differ diff --git a/tests/goldens/grids/cases.csv b/tests/goldens/grids/cases.csv new file mode 100644 index 0000000..e9cb8a2 --- /dev/null +++ b/tests/goldens/grids/cases.csv @@ -0,0 +1,13 @@ +case,model,arguments,seed,final_time,xdim,ydim,attribute +life_glider,Life,pattern=glider,,100,30,30,state +life_rpentomino,Life,pattern=rPentomino,,100,30,30,state +life_random,Life,,1,100,30,30,state +parity,Parity,,,500,50,50,state +excitable,Excitable,,,500,50,50,state +wolfram_30,Wolfram,rule=30,,55,111,56,state +wolfram_90,Wolfram,rule=90,,55,111,56,state +wolfram_110,Wolfram,rule=110,,55,111,56,state +fire,Fire,,1,100,60,60,state +anneal,Anneal,,1,100,80,80,state +oscillator,Oscillator,,1,400,50,50,state +parasit,Parasit,,1,500,50,50,state diff --git a/tests/goldens/grids/excitable.csv.gz b/tests/goldens/grids/excitable.csv.gz new file mode 100644 index 0000000..3cdb2db Binary files /dev/null and b/tests/goldens/grids/excitable.csv.gz differ diff --git a/tests/goldens/grids/fire.csv.gz b/tests/goldens/grids/fire.csv.gz new file mode 100644 index 0000000..e58c175 Binary files /dev/null and b/tests/goldens/grids/fire.csv.gz differ diff --git a/tests/goldens/grids/life_glider.csv.gz b/tests/goldens/grids/life_glider.csv.gz new file mode 100644 index 0000000..437828a Binary files /dev/null and b/tests/goldens/grids/life_glider.csv.gz differ diff --git a/tests/goldens/grids/life_random.csv.gz b/tests/goldens/grids/life_random.csv.gz new file mode 100644 index 0000000..fbac140 Binary files /dev/null and b/tests/goldens/grids/life_random.csv.gz differ diff --git a/tests/goldens/grids/life_rpentomino.csv.gz b/tests/goldens/grids/life_rpentomino.csv.gz new file mode 100644 index 0000000..878a2cb Binary files /dev/null and b/tests/goldens/grids/life_rpentomino.csv.gz differ diff --git a/tests/goldens/grids/oscillator.csv.gz b/tests/goldens/grids/oscillator.csv.gz new file mode 100644 index 0000000..b3e7b0e Binary files /dev/null and b/tests/goldens/grids/oscillator.csv.gz differ diff --git a/tests/goldens/grids/parasit.csv.gz b/tests/goldens/grids/parasit.csv.gz new file mode 100644 index 0000000..db75be6 Binary files /dev/null and b/tests/goldens/grids/parasit.csv.gz differ diff --git a/tests/goldens/grids/parity.csv.gz b/tests/goldens/grids/parity.csv.gz new file mode 100644 index 0000000..c934c85 Binary files /dev/null and b/tests/goldens/grids/parity.csv.gz differ diff --git a/tests/goldens/grids/wolfram_110.csv.gz b/tests/goldens/grids/wolfram_110.csv.gz new file mode 100644 index 0000000..a2f5c24 Binary files /dev/null and b/tests/goldens/grids/wolfram_110.csv.gz differ diff --git a/tests/goldens/grids/wolfram_30.csv.gz b/tests/goldens/grids/wolfram_30.csv.gz new file mode 100644 index 0000000..616170e Binary files /dev/null and b/tests/goldens/grids/wolfram_30.csv.gz differ diff --git a/tests/goldens/grids/wolfram_90.csv.gz b/tests/goldens/grids/wolfram_90.csv.gz new file mode 100644 index 0000000..a4bb217 Binary files /dev/null and b/tests/goldens/grids/wolfram_90.csv.gz differ diff --git a/tests/goldens/make_reference.py b/tests/goldens/make_reference.py new file mode 100644 index 0000000..5fc20aa --- /dev/null +++ b/tests/goldens/make_reference.py @@ -0,0 +1,107 @@ +""" +Extract the reference values used by ``tests/test_goldens.py`` from the +TerraME ``ca`` goldens of ``terrame-goldens``. + +- ``grids/.csv.gz``: the grid goldens (state of every cell; columns + ``step, y, x0, x1, ...``), kept only at the steps in :func:`grid_steps`. +- ``ensembles.csv.gz``: for the stochastic models, the number of cells in + each state for seeds 1 to 50, at the steps in ``ENSEMBLE_STEPS`` + (columns ``case, value, step, seed, count``; zeros included). +- ``reference.json``: source version, package commit, steps kept and the + SHA-256 of each source file. + + python tests/goldens/make_reference.py /path/to/terrame-goldens +""" +from __future__ import annotations + +import csv +import gzip +import hashlib +import io +import json +import sys +from collections import defaultdict +from pathlib import Path + +HERE = Path(__file__).resolve().parent + +ENSEMBLE_STEPS = { + "fire": (5, 10, 20, 50, 100), + "solid_diffusion": (10, 50, 100, 200, 400), + "growth": (5, 10, 20, 50, 100), + "snow": (10, 50, 100, 300, 600), + "interspecific_competition": (5, 10, 50, 100, 200), +} + + +def grid_steps(final: int) -> list[int]: + """Every step up to 10, then a few more up to the final step.""" + steps = set(range(min(final, 10) + 1)) | {final} + steps |= {s for s in (20, 50, 100, 200, 300, 400, 500) if s <= final} + return sorted(steps) + + +def write_gz(text: str, dest: Path) -> None: + dest.parent.mkdir(parents=True, exist_ok=True) + buf = io.BytesIO() + with gzip.GzipFile(filename="", mode="wb", fileobj=buf, mtime=0, compresslevel=9) as gz: + gz.write(text.encode()) + dest.write_bytes(buf.getvalue()) + + +def sha256(path: Path) -> str: + return hashlib.sha256(path.read_bytes()).hexdigest() + + +def main(root: Path) -> None: + src = root / "goldens" / "ca" + sources: dict[str, str] = {} + grids: dict[str, list[int]] = {} + + with open(src / "grids" / "cases.csv", newline="") as f: + cases = list(csv.DictReader(f)) + for c in cases: + path = src / "grids" / f"{c['case']}.csv.gz" + sources[path.relative_to(root).as_posix()] = sha256(path) + keep = grid_steps(int(c["final_time"])) + with gzip.open(path, "rt", newline="") as f: + lines = f.read().splitlines() + out = [lines[0]] + [ln for ln in lines[1:] if int(ln.split(",", 1)[0]) in keep] + write_gz("\n".join(out) + "\n", HERE / "grids" / f"{c['case']}.csv.gz") + grids[c["case"]] = keep + (HERE / "grids" / "cases.csv").write_bytes((src / "grids" / "cases.csv").read_bytes()) + + rows = [] + for case, steps in ENSEMBLE_STEPS.items(): + path = src / "ensembles" / f"{case}.csv.gz" + sources[path.relative_to(root).as_posix()] = sha256(path) + counts: dict[tuple[str, int, int], int] = defaultdict(int) + values, seeds = set(), set() + with gzip.open(path, "rt", newline="") as f: + for r in csv.DictReader(f): + values.add(r["value"]) + seeds.add(int(r["seed"])) + if int(r["step"]) in steps: + counts[(r["value"], int(r["step"]), int(r["seed"]))] = int(r["count"]) + for v in sorted(values): + for t in steps: + for s in sorted(seeds): + rows.append(f"{case},{v},{t},{s},{counts[(v, t, s)]}") + write_gz("case,value,step,seed,count\n" + "\n".join(rows) + "\n", HERE / "ensembles.csv.gz") + + manifest = json.loads((root / "manifest.json").read_text()) + meta = { + "source": "https://github.com/LambdaGeo/terrame-goldens", + "version": manifest["version"], + "terrame": manifest["terrame"], + "packages": [p for p in manifest["packages"] if p["name"] == "ca"], + "files": sources, + "grid_steps": grids, + "ensemble_steps": {k: list(v) for k, v in ENSEMBLE_STEPS.items()}, + } + (HERE / "reference.json").write_text(json.dumps(meta, indent=2) + "\n") + print(f"{len(grids)} grids, {len(rows)} ensemble values written") + + +if __name__ == "__main__": + main(Path(sys.argv[1])) diff --git a/tests/goldens/reference.json b/tests/goldens/reference.json new file mode 100644 index 0000000..be5dcc4 --- /dev/null +++ b/tests/goldens/reference.json @@ -0,0 +1,277 @@ +{ + "source": "https://github.com/LambdaGeo/terrame-goldens", + "version": "v0.1.0", + "terrame": "2.0.1", + "packages": [ + { + "name": "ca", + "repository": "https://github.com/TerraME/ca", + "commit": "854a5c405478d712013a0baffd382f438f367642" + } + ], + "files": { + "goldens/ca/grids/life_glider.csv.gz": "4864b1ff009f11601706ea33f2d226fa87af3fab911b8a3e21595ac10bbe497e", + "goldens/ca/grids/life_rpentomino.csv.gz": "df597c76e5c271af8576577a94f69025d5cc1e2254550320982fb5333b5254e9", + "goldens/ca/grids/life_random.csv.gz": "20e28a42c86afa85b99d974cbaac3365c55ec6487e0eb9d0b9b1403f0f58660a", + "goldens/ca/grids/parity.csv.gz": "a3412ead9a48d67b2301f9e20c5d04ddb918d513a56d4b75ba1a152d75c0be52", + "goldens/ca/grids/excitable.csv.gz": "bd92ca276cf3424f1439d2b1a78533bc605cfe1e1834434649a573ed9249c477", + "goldens/ca/grids/wolfram_30.csv.gz": "12a74a836d3051da613de58e101f0bbc700f796b4b9ce4ad4f1949c77fc9d72d", + "goldens/ca/grids/wolfram_90.csv.gz": "9b15c41617f909d3355ac46e194aa314dfb8a10ea1a34eb8e8b89d9766f97303", + "goldens/ca/grids/wolfram_110.csv.gz": "5757b0593bdf25759b55ecdf6b61f8721c1538ad57ea910b3b629806e0d5a250", + "goldens/ca/grids/fire.csv.gz": "973a8b9e5024789ac944736f9c376c63acd220014f4cac328ef4e8b6ed8829ac", + "goldens/ca/grids/anneal.csv.gz": "f1d918d6ef51e0db341f3c07676ac88b491dbcf1a47148e1b01bb1ac859613c6", + "goldens/ca/grids/oscillator.csv.gz": "dd26dd943b1d5fa8bfcf82ebaebf697436ce1827f452da715404b394c57a69d7", + "goldens/ca/grids/parasit.csv.gz": "814b9c11648545d337bfe81522641abba41eaaf0a8eba27df7ee8b9508bf4328", + "goldens/ca/ensembles/fire.csv.gz": "2fcdfa432da8e250fd4d7d2a40eec186e36b82e1d5536ec2051762f8657ed24c", + "goldens/ca/ensembles/solid_diffusion.csv.gz": "dc5b71f89981e657e3eb5e4c6a021ddbafeebf7c74e05e4940e75069cb675480", + "goldens/ca/ensembles/growth.csv.gz": "249c76c7645bd269c0fa3ec6d856b9f3e7d9babb23f4452c01ce7a3f841e3f43", + "goldens/ca/ensembles/snow.csv.gz": "c57ed9260d69f7fd25b7571a4c9f66580b555c3078af63ce5367b9f5257d1875", + "goldens/ca/ensembles/interspecific_competition.csv.gz": "b5dcbfb0fe07ef3043cedcbb3594c1b54fcebb3a40371fd53743039edb992697" + }, + "grid_steps": { + "life_glider": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 100 + ], + "life_rpentomino": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 100 + ], + "life_random": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 100 + ], + "parity": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 100, + 200, + 300, + 400, + 500 + ], + "excitable": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 100, + 200, + 300, + 400, + 500 + ], + "wolfram_30": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 55 + ], + "wolfram_90": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 55 + ], + "wolfram_110": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 55 + ], + "fire": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 100 + ], + "anneal": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 100 + ], + "oscillator": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 100, + 200, + 300, + 400 + ], + "parasit": [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 20, + 50, + 100, + 200, + 300, + 400, + 500 + ] + }, + "ensemble_steps": { + "fire": [ + 5, + 10, + 20, + 50, + 100 + ], + "solid_diffusion": [ + 10, + 50, + 100, + 200, + 400 + ], + "growth": [ + 5, + 10, + 20, + 50, + 100 + ], + "snow": [ + 10, + 50, + 100, + 300, + 600 + ], + "interspecific_competition": [ + 5, + 10, + 50, + 100, + 200 + ] + } +} diff --git a/tests/test_goldens.py b/tests/test_goldens.py new file mode 100644 index 0000000..de07f8f --- /dev/null +++ b/tests/test_goldens.py @@ -0,0 +1,279 @@ +""" +The models against TerraME itself. + +``goldens/`` holds an extract of the ``ca`` suite of ``terrame-goldens``: what +TerraME 2.0.1 computed running the models of the ``ca`` package (see +``goldens/README.md``). Two kinds of tests: + +- **Grids**: models whose rule is deterministic start from TerraME's state at + step 0 and must match it **cell by cell** at every step kept in the + extract. +- **Ensembles**: models that draw random numbers while they run are compared + by distribution. At each step kept, the number of cells in each state over + ``SEEDS`` runs and over TerraME's 50 runs must pass a two-sample + Kolmogorov-Smirnov test (p >= ``P_MIN``) and their means must not differ by + more than ``Z_MAX`` standard errors. + +The vector models take 0.2 to 2 s per step, so by default the grids are +compared over the first ``QUICK_STEPS`` steps (the raster models and Wolfram +over all of them) and only the fast ensembles run, in about a minute. Set +``DISSMODEL_GOLDEN_FULL=1`` to compare every step kept and every ensemble +(Growth alone takes hours). +""" +from __future__ import annotations + +import csv +import gzip +import json +import os +import random +from collections import defaultdict +from collections.abc import Callable +from functools import cache +from pathlib import Path +from typing import Any + +import numpy as np +import pytest +from dissmodel.core import Environment, Model +from dissmodel.geo import parse_idx, vector_grid +from dissmodel.geo.raster.backend import RasterBackend + +import dissmodel_ca.models as ca + +pytestmark = pytest.mark.golden + +GOLDENS = Path(__file__).resolve().parent / "goldens" +REFERENCE = json.loads((GOLDENS / "reference.json").read_text()) +FULL = os.environ.get("DISSMODEL_GOLDEN_FULL", "") not in ("", "0") +QUICK_STEPS = 10 +SEEDS = range(1, 21) +Z_MAX = 4.0 +P_MIN = 0.001 + + +# -- grids --------------------------------------------------------------------------------- + + +def load_grid(case: str) -> dict[int, np.ndarray]: + """``{step: (y, x) array of TerraME's values as str}``.""" + with gzip.open(GOLDENS / "grids" / f"{case}.csv.gz", "rt", newline="") as f: + rows = list(csv.reader(f)) + nx = len(rows[0]) - 2 + by_step: dict[int, dict[int, list[str]]] = defaultdict(dict) + for r in rows[1:]: + by_step[int(r[0])][int(r[1])] = r[2:] + return { + t: np.array([g[y] for y in range(len(g))], dtype=object).reshape(len(g), nx) + for t, g in by_step.items() + } + + +class _Recorder(Model): + """Reads the model after each step listed (created after the model, so it runs after it).""" + + def setup(self, read: Callable[[], Any], steps: set[int]) -> None: # type: ignore[override] + self.read, self.steps, self.frames = read, steps, {} + + def execute(self) -> None: + t = int(self.env.now()) + if t in self.steps: + self.frames[t] = self.read() + + +def run_vector(factory: Callable[[Any], Any], init: np.ndarray, steps: set[int]) -> dict[int, np.ndarray]: + ny, nx = init.shape + env = Environment(start_time=1, end_time=max(steps)) + gdf = vector_grid(dimension=(nx, ny), resolution=1, attrs={"state": 0}) + model = factory(gdf) + pos = [parse_idx(i) for i in model.gdf.index] + rows = np.array([p.row for p in pos]) + cols = np.array([p.col for p in pos]) + model.gdf["state"] = init[rows, cols] + + def read() -> np.ndarray: + out = np.zeros((ny, nx), dtype=int) + out[rows, cols] = model.gdf["state"].to_numpy() + return out + + rec = _Recorder(read=read, steps=steps) + env.run() + return rec.frames + + +def run_raster(factory: Callable[[Any], Any], init: np.ndarray, steps: set[int]) -> dict[int, np.ndarray]: + env = Environment(start_time=1, end_time=max(steps)) + backend = RasterBackend(shape=init.shape) + backend.set("state", init.astype(np.int8)) + factory(backend) + rec = _Recorder(read=lambda: backend.arrays["state"].astype(int).copy(), steps=steps) + env.run() + return rec.frames + + +LIFE = {"dead": 0, "alive": 1} +FIRE = {"forest": 0, "burning": 1, "burned": 2, "empty": 3} + +#: id, golden, values -> codes (None: integers), runner, factory, quick (compare every step by default) +GRID_CASES = [ + ("life_glider", "life_glider", LIFE, run_vector, lambda g: ca.GameOfLife(gdf=g), False), + ("life_rpentomino", "life_rpentomino", LIFE, run_vector, lambda g: ca.GameOfLife(gdf=g), False), + ("life_random", "life_random", LIFE, run_vector, lambda g: ca.GameOfLife(gdf=g), False), + ("life_random-raster", "life_random", LIFE, run_raster, + lambda b: ca.GameOfLifeRaster(backend=b), True), + ("parity", "parity", {"off": 0, "on": 1}, run_vector, lambda g: ca.Parity(gdf=g, dim=50), False), + ("excitable", "excitable", None, run_vector, lambda g: ca.Excitable(gdf=g, dim=50), False), + ("wolfram_30", "wolfram_30", LIFE, run_vector, lambda g: ca.Wolfram(gdf=g, rule_number=30), True), + ("wolfram_90", "wolfram_90", LIFE, run_vector, lambda g: ca.Wolfram(gdf=g, rule_number=90), True), + ("wolfram_110", "wolfram_110", LIFE, run_vector, lambda g: ca.Wolfram(gdf=g, rule_number=110), True), + ("fire", "fire", FIRE, run_vector, lambda g: ca.FireModel(gdf=g), False), + ("fire-raster", "fire", FIRE, run_raster, lambda b: ca.FireModelRaster(backend=b), True), + ("anneal", "anneal", {"L": 0, "R": 1}, run_vector, lambda g: ca.Anneal(gdf=g), False), + ("oscillator", "oscillator", None, run_vector, lambda g: ca.Oscillator(gdf=g), False), + ("parasit", "parasit", None, run_vector, lambda g: ca.Parasit(gdf=g), False), +] + + +def codes(grid: np.ndarray, mapping: dict[str, int] | None) -> np.ndarray: + if mapping is None: + return grid.astype(int) + return np.vectorize(mapping.__getitem__, otypes=[int])(grid) + + +@pytest.mark.parametrize( + "golden,mapping,runner,factory,quick", + [c[1:] for c in GRID_CASES], + ids=[c[0] for c in GRID_CASES], +) +def test_grid_matches_terrame_cell_by_cell(golden, mapping, runner, factory, quick) -> None: + ref = {t: codes(g, mapping) for t, g in load_grid(golden).items()} + steps = {t for t in ref if t > 0 and (FULL or quick or t <= QUICK_STEPS)} + ours = runner(factory, ref[0], steps) + for t in sorted(steps): + diff = int((ours[t] != ref[t]).sum()) + assert diff == 0, f"{golden}: {diff} cells differ from TerraME at step {t}" + + +# -- ensembles ------------------------------------------------------------------------------- + + +def ks_pvalue(a: np.ndarray, b: np.ndarray) -> float: + """Two-sample Kolmogorov-Smirnov test, asymptotic p-value (Numerical Recipes ``kstwo``).""" + grid = np.union1d(a, b) + d = np.abs( + np.searchsorted(np.sort(a), grid, side="right") / len(a) + - np.searchsorted(np.sort(b), grid, side="right") / len(b) + ).max() + en = np.sqrt(len(a) * len(b) / (len(a) + len(b))) + lam = (en + 0.12 + 0.11 / en) * d + if lam < 1e-3: + return 1.0 + k = np.arange(1, 101) + return float(np.clip(2 * np.sum((-1.0) ** (k - 1) * np.exp(-2 * (k * lam) ** 2)), 0, 1)) + + +@cache +def ensemble_reference() -> dict[tuple[str, str, int], np.ndarray]: + values: dict[tuple[str, str, int], list[int]] = defaultdict(list) + with gzip.open(GOLDENS / "ensembles.csv.gz", "rt", newline="") as f: + for r in csv.DictReader(f): # sorted by seed within (case, value, step) + values[(r["case"], r["value"], int(r["step"]))].append(int(r["count"])) + return {k: np.array(v) for k, v in values.items()} + + +Counts = dict[str, int] + + +def _counts(state: np.ndarray, mapping: dict[str, int]) -> Counts: + return {value: int((state == code).sum()) for value, code in mapping.items()} + + +def _vector_ensemble(n: int, mapping: dict[str, int], factory: Callable[[Any, int], Any]): + def run(seed: int, steps: set[int]) -> dict[int, Counts]: + random.seed(seed) + np.random.seed(seed) + env = Environment(start_time=1, end_time=max(steps)) + gdf = vector_grid(dimension=(n, n), resolution=1, attrs={"state": 0}) + model = factory(gdf, seed) + model.initialize() + cols = np.array([parse_idx(i).col for i in model.gdf.index]) + + def read() -> Counts: + state = model.gdf["state"].to_numpy() + out = _counts(state, mapping) + if "atom1" in mapping: # SolidDiffusion: atoms across the middle column + middle = n // 2 + out["atom1_right"] = int(((state == mapping["atom1"]) & (cols > middle)).sum()) + out["atom2_left"] = int(((state == mapping["atom2"]) & (cols < middle)).sum()) + return out + + rec = _Recorder(read=read, steps=steps) + env.run() + return rec.frames + return run + + +def _fire_raster(seed: int, steps: set[int]) -> dict[int, Counts]: + env = Environment(start_time=1, end_time=max(steps)) + backend = RasterBackend(shape=(60, 60)) + model = ca.FireModelRaster(backend=backend, seed=seed) + model.initialize() + rec = _Recorder(read=lambda: _counts(backend.arrays["state"], FIRE), steps=steps) + env.run() + return rec.frames + + +SOLID = {"atom1": 0, "atom2": 1, "vacancy": 2} +SPECIES = {"Lolium": 0, "Agrostis": 1, "Holcus": 2, "Poa": 3, "Cynosurus": 4} + +#: case -> (runner, runs by default). SolidDiffusion only swaps atoms, so its +#: counts per state never change; ``atom1_right`` and ``atom2_left`` (atoms +#: past the middle column) measure how far they mixed. +ENSEMBLES: dict[str, tuple[Callable[[int, set[int]], dict[int, Counts]], bool]] = { + "fire": (_fire_raster, True), + "solid_diffusion": ( + _vector_ensemble(31, SOLID, lambda g, s: ca.SolidDiffusion(gdf=g, dim=31, seed=s)), True), + "growth": ( + _vector_ensemble(100, {"alive": 1, "empty": 0}, lambda g, s: ca.Growth(gdf=g, dim=100)), False), + "snow": ( + _vector_ensemble(50, {"snow": 1, "empty": 0}, + lambda g, s: ca.Snow(gdf=g, dim=50, end_time=600)), False), + "interspecific_competition": ( + _vector_ensemble(40, SPECIES, + lambda g, s: ca.InterspecificCompetition(gdf=g, dim=40, seed=s)), False), +} + + +@cache +def ensemble_runs(case: str) -> dict[tuple[str, int], np.ndarray]: + runner, _ = ENSEMBLES[case] + steps = set(REFERENCE["ensemble_steps"][case]) + out: dict[tuple[str, int], list[int]] = defaultdict(list) + for seed in SEEDS: + for t, counts in runner(seed, steps).items(): + for value, count in counts.items(): + out[(value, t)].append(count) + return {k: np.array(v) for k, v in out.items()} + + +ENSEMBLE_CASES = sorted(k for k in ensemble_reference() if k[0] in ENSEMBLES) + + +@pytest.mark.parametrize( + "case,value,step", ENSEMBLE_CASES, ids=[f"{c}-{v}-t{s}" for c, v, s in ENSEMBLE_CASES] +) +def test_ensemble_matches_terrame(case: str, value: str, step: int) -> None: + if not (FULL or ENSEMBLES[case][1]): + pytest.skip("slow ensemble; set DISSMODEL_GOLDEN_FULL=1") + tm = ensemble_reference()[(case, value, step)] + ours = ensemble_runs(case)[(value, step)] + se = np.sqrt(tm.var(ddof=1) / len(tm) + ours.var(ddof=1) / len(ours)) + diff = ours.mean() - tm.mean() + p = ks_pvalue(tm, ours) + msg = (f"{case} '{value}' at t={step}: TerraME {tm.mean():.1f} ± {tm.std():.1f}, " + f"dissmodel-ca {ours.mean():.1f} ± {ours.std():.1f}, KS p = {p:.2g}") + assert p >= P_MIN, msg + if se > 0: + assert abs(diff) <= Z_MAX * se, f"{msg}, z = {diff / se:.1f}" + else: + assert diff == 0, msg