fix: Life wrap, raster Fire neighborhood and SolidDiffusion

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

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).

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

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)",
+]

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()

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:
         """

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)

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`."""