fix: Daisyworld and MonoLake as in TerraME; every model bit-identical…

CHANGELOG.md

@@ -2,6 +2,27 @@
 
 ## [Unreleased]
 
+- **Fix:** `Daisyworld` now follows TerraME: after updating the white area,
+  it computes the planet albedo again for the black daisies. The trajectory
+  changes (up to 0.61 in area); the equilibrium is the same. New attributes
+  `temp_near_white` and `temp_near_black`, as in TerraME.
+- **Fix:** `MonoLake` reads the lake surface and elevation with TerraME's
+  `sci.Spline` (a Catmull-Rom curve sampled at 1000 points per segment and
+  read linearly; new `dissmodel_sysdyn.models._spline.Spline`) instead of
+  scipy's `CubicSpline`, which differed by up to 26 KAF in 50 steps. As in
+  TerraME, a volume outside the curve raises an error instead of
+  extrapolating. scipy, which was used without being a declared
+  dependency, is no longer needed.
+- `PredatorPrey`, `RoomTemperature`, `SIR` and `MonoLake` do their
+  arithmetic in TerraME's order, so that the results are the same doubles.
+- `tests/test_goldens.py` (marker `golden`): every model against TerraME
+  2.0.1 runs from `terrame-goldens` v0.1.0. All 16 deterministic cases
+  (13 models, some with two parameter sets) give **the same doubles as
+  TerraME at every step**; RandomWalk is compared by distribution.
+  `tests/goldens/` holds the goldens, their SHA-256 provenance and the
+  script that copies them. The strict `xfail`s of `Daisyworld` and
+  `MonoLake` in `test_terrame_parity.py` are removed.
+
 - **Fix:** `Tub` kept its counter in `_step`, the attribute where dissmodel's
   `Model` stores the interval between ticks, so the interval grew at every
   execution (a run from t=1 to 30 executed 8 times). The counter is now

README.md

@@ -25,6 +25,16 @@ The library includes a wide range of System Dynamics models:
 *   **Environment**: Mono Lake Water Balance.
 *   **Stochasticity**: Random Walk.
 
+## Verification against TerraME
+
+Every model was run in TerraME 2.0.1 with the `sysdyn` package
+([`terrame-goldens`](https://github.com/LambdaGeo/terrame-goldens) v0.1.0).
+For the 16 deterministic cases (13 models, some with two parameter sets),
+the Python models give **the same doubles as TerraME at every step**, bit for
+bit, for every variable TerraME records; RandomWalk is compared by
+distribution. `tests/test_goldens.py` runs these comparisons
+(`pytest -m golden`, under a second).
+
 ## Examples
 
 The project includes several examples demonstrating how to run the models via CLI, Streamlit, or Jupyter Notebooks.

examples/notebooks/sysdyn_mono_lake.ipynb

@@ -88,17 +88,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import numpy as np\n",
-    "from scipy.interpolate import CubicSpline\n",
     "from dissmodel.core import Model\n",
     "from dissmodel.visualization import track_plot\n",
+    "from dissmodel_sysdyn.models._spline import Spline  # TerraME's sci.Spline\n",
     "\n",
-    "_VOLUME_KAF  = np.array([0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000], dtype=float)\n",
-    "_SURFACE_KM2 = np.array([0, 24.7, 35.3, 48.6, 54.3, 57.2, 61.6, 66.0, 69.9], dtype=float)\n",
-    "_ELEV_FT     = np.array([6224, 6335, 6369, 6392, 6412, 6430, 6447, 6463, 6477], dtype=float)\n",
+    "_VOLUME_KAF  = [0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000]\n",
+    "_SURFACE_KM2 = [0, 24.7, 35.3, 48.6, 54.3, 57.2, 61.6, 66.0, 69.9]\n",
+    "_ELEV_FT     = [6224, 6335, 6369, 6392, 6412, 6430, 6447, 6463, 6477]\n",
     "\n",
-    "_water_surface   = CubicSpline(_VOLUME_KAF, _SURFACE_KM2, extrapolate=True)\n",
-    "_water_elevation = CubicSpline(_VOLUME_KAF, _ELEV_FT,     extrapolate=True)\n",
+    "# a Catmull-Rom curve sampled at 1000 points per segment, read linearly\n",
+    "_water_surface   = Spline(_VOLUME_KAF, _SURFACE_KM2, steps=1000)\n",
+    "_water_elevation = Spline(_VOLUME_KAF, _ELEV_FT,     steps=1000)\n",
     "\n",
     "@track_plot(\"level\", \"blue\")\n",
     "@track_plot(\"water_in_lake\", \"cyan\")\n",
@@ -109,11 +109,11 @@
     "        self.evap_rate, self.other_out, self.export = evap_rate, other_out, export\n",
     "\n",
     "    def execute(self) -> None:\n",
-    "        surface = float(_water_surface(self.water_in_lake))\n",
+    "        surface = _water_surface.value(self.water_in_lake)\n",
     "        total_in = surface * self.prec_rate + self.runoff + self.other_in\n",
     "        total_out = surface * self.evap_rate + self.export + self.other_out\n",
-    "        self.water_in_lake += total_in - total_out\n",
-    "        self.level = float(_water_elevation(self.water_in_lake))"
+    "        self.water_in_lake = self.water_in_lake + total_in - total_out\n",
+    "        self.level = _water_elevation.value(self.water_in_lake)"
    ]
   },
   {

pyproject.toml

@@ -35,3 +35,6 @@ Homepage = "https://github.com/DisSModel/dissmodel-sysdyn"
 
 [tool.pytest.ini_options]
 testpaths = ["tests"]
+markers = [
+    "golden: compares the models with TerraME runs (tests/goldens)",
+]

src/dissmodel_sysdyn/models/_spline.py

@@ -0,0 +1,87 @@
+"""
+The ``Spline`` of TerraME's ``sci`` package (commit 42b0b3e, ``lua/Spline.lua``).
+
+``sci.Spline`` is not a cubic spline through the points: it samples a
+Catmull-Rom curve at ``steps`` points per segment (the first and last points
+repeated as end controls) and reads it by linear interpolation between the
+samples. :class:`Spline` reproduces both steps, with the same floating-point
+operations, so that a model that uses it matches TerraME to the last bit.
+
+One consequence of following Lua: the loop over ``t`` adds ``1 / steps``
+to ``t`` and, for some ``steps`` (1000 among them), the sum ends just below
+1, so the last sample stops short of the last point (at 7999.5 instead of
+8000 in MonoLake) and a larger ``x`` raises an error, in TerraME as here.
+"""
+from __future__ import annotations
+
+from bisect import bisect_left
+from collections.abc import Sequence
+
+
+class Spline:
+    """
+    TerraME ``sci.Spline``: a sampled Catmull-Rom curve read linearly.
+
+    Parameters
+    ----------
+    xs, ys : sequence of float
+        Control points, with increasing ``xs``.
+    steps : int, optional
+        Samples per segment, by default 10 (as in ``sci``).
+
+    Raises
+    ------
+    ValueError
+        From :meth:`value`, for an ``x`` outside the range of the points, as
+        TerraME does.
+    """
+
+    def __init__(self, xs: Sequence[float], ys: Sequence[float], steps: int = 10) -> None:
+        points = list(zip((float(x) for x in xs), (float(y) for y in ys), strict=True))
+        self.values = self._build(points, steps)
+        self.xs = [p[0] for p in self.values]
+
+    @staticmethod
+    def _build(points: list[tuple[float, float]], steps: int) -> list[tuple[float, float]]:
+        if len(points) < 3:
+            return points
+        count = len(points) - 1
+        values: list[tuple[float, float]] = []
+        for i in range(1, count + 1):  # 1-based, as in Lua
+            if i == 1:
+                p0, p1, p2, p3 = points[0], points[0], points[1], points[2]
+            elif i == count:
+                p0, p1, p2, p3 = points[-3], points[-2], points[-1], points[-1]
+            else:
+                p0, p1, p2, p3 = points[i - 2], points[i - 1], points[i], points[i + 1]
+            # Lua's ``for t = 0, 1, 1 / steps`` adds the step to t each time
+            t = 0.0
+            dt = 1 / steps
+            while t <= 1:
+                x = 0.5 * ((2 * p1[0]) + (p2[0] - p0[0]) * t
+                           + (2 * p0[0] - 5 * p1[0] + 4 * p2[0] - p3[0]) * t * t
+                           + (3 * p1[0] - p0[0] - 3 * p2[0] + p3[0]) * t * t * t)
+                y = 0.5 * ((2 * p1[1]) + (p2[1] - p0[1]) * t
+                           + (2 * p0[1] - 5 * p1[1] + 4 * p2[1] - p3[1]) * t * t
+                           + (3 * p1[1] - p0[1] - 3 * p2[1] + p3[1]) * t * t * t)
+                if not values or (values[-1][0] != x and values[-1][1] != y):
+                    values.append((x, y))
+                t += dt
+        return values
+
+    def value(self, x: float) -> float:
+        """Linear interpolation between the two samples around ``x``."""
+        xs = self.xs
+        if not xs[0] <= x <= xs[-1]:
+            raise ValueError(f"Spline: {x} is outside the range [{xs[0]}, {xs[-1]}].")
+        i = bisect_left(xs, x)
+        if xs[i] == x:
+            return self.values[i][1]
+        (x0, y0), (x1, y1) = self.values[i - 1], self.values[i]
+        w = (x - x0) / (x1 - x0)
+        return (1 - w) * y0 + w * y1
+
+    __call__ = value
+
+
+__all__ = ["Spline"]

src/dissmodel_sysdyn/models/daisyworld.py

@@ -65,6 +65,10 @@ class Daisyworld(Model):
 
     Notes
     -----
+    The step follows the TerraME ``sysdyn`` model line by line: the white
+    area is updated first, and the planet albedo used for the black daisies
+    is computed again with the new white area.
+
     The ``@track_plot`` decorators register ``white_area``,
     ``black_area``, ``empty_area``, and ``ave_temp`` for automatic
     live plotting.
@@ -89,6 +93,12 @@ class Daisyworld(Model):
     #: Current average planet temperature (K).
     ave_temp: float
 
+    #: Temperature near the white daisies in the last step (K).
+    temp_near_white: float
+
+    #: Temperature near the black daisies in the last step (K).
+    temp_near_black: float
+
     #: Solar luminosity (fractional).
     sun_luminosity: float
 
@@ -138,15 +148,16 @@ def _planet_albedo(self) -> float:
         )
 
     def execute(self) -> None:
-        pa = self._planet_albedo()
-        self.ave_temp = _planet_temp(self.sun_luminosity, pa)
+        self.ave_temp = _planet_temp(self.sun_luminosity, self._planet_albedo())
 
-        temp_white     = _local_temp(self.ave_temp, pa, self.white_albedo)
-        white_growth   = _daisy_growth_rate(temp_white) * self.empty_area
+        self.temp_near_white = _local_temp(self.ave_temp, self._planet_albedo(), self.white_albedo)
+        white_growth = _daisy_growth_rate(self.temp_near_white) * self.empty_area
         self.white_area += self.white_area * (white_growth - self.decay_rate)
 
-        temp_black     = _local_temp(self.ave_temp, pa, self.black_albedo)
-        black_growth   = _daisy_growth_rate(temp_black) * self.empty_area
+        # As in TerraME, the planet albedo is computed again here, with the
+        # new white area (and the black and empty areas of the start of the step)
+        self.temp_near_black = _local_temp(self.ave_temp, self._planet_albedo(), self.black_albedo)
+        black_growth = _daisy_growth_rate(self.temp_near_black) * self.empty_area
         self.black_area += self.black_area * (black_growth - self.decay_rate)
 
         self.empty_area = self.planet_area - (self.white_area + self.black_area)

src/dissmodel_sysdyn/models/mono_lake.py

@@ -1,17 +1,17 @@
 from __future__ import annotations
 
-from scipy.interpolate import CubicSpline
-import numpy as np
-
 from dissmodel.core import Model
 from dissmodel.visualization import track_plot
 
-_VOLUME_KAF  = np.array([0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000], dtype=float)
-_SURFACE_KM2 = np.array([0, 24.7, 35.3, 48.6, 54.3, 57.2, 61.6, 66.0, 69.9], dtype=float)
-_ELEV_FT     = np.array([6224, 6335, 6369, 6392, 6412, 6430, 6447, 6463, 6477], dtype=float)
+from ._spline import Spline
+
+_VOLUME_KAF  = [0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000]
+_SURFACE_KM2 = [0, 24.7, 35.3, 48.6, 54.3, 57.2, 61.6, 66.0, 69.9]
+_ELEV_FT     = [6224, 6335, 6369, 6392, 6412, 6430, 6447, 6463, 6477]
 
-_water_surface   = CubicSpline(_VOLUME_KAF, _SURFACE_KM2, extrapolate=True)
-_water_elevation = CubicSpline(_VOLUME_KAF, _ELEV_FT,     extrapolate=True)
+# As in TerraME: sci.Spline{points = ..., steps = 1000}
+_water_surface   = Spline(_VOLUME_KAF, _SURFACE_KM2, steps=1000)
+_water_elevation = Spline(_VOLUME_KAF, _ELEV_FT,     steps=1000)
 
 
 @track_plot("level", "blue")
@@ -43,6 +43,14 @@ class MonoLake(Model):
     export : float, optional
         Annual water exported in KAF/year, by default 100.0.
 
+    Notes
+    -----
+    The lake surface and elevation are read from the volume through
+    TerraME's ``sci.Spline`` (:class:`~dissmodel_sysdyn.models._spline.Spline`):
+    a Catmull-Rom curve sampled at 1000 points per segment and read by
+    linear interpolation, as in the TerraME model. As there, a volume
+    outside 0 to about 7999.5 KAF raises an error instead of extrapolating.
+
     Examples
     --------
     >>> from dissmodel.core import Environment
@@ -96,7 +104,7 @@ def setup(
         self.export        = export
 
     def _surface(self) -> float:
-        return float(_water_surface(self.water_in_lake))
+        return _water_surface.value(self.water_in_lake)
 
     def _total_input(self) -> float:
         return self._surface() * self.prec_rate + self.runoff + self.other_in
@@ -105,8 +113,9 @@ def _total_output(self) -> float:
         return self._surface() * self.evap_rate + self.export + self.other_out
 
     def execute(self) -> None:
-        self.water_in_lake += self._total_input() - self._total_output()
-        self.level = float(_water_elevation(self.water_in_lake))
+        # same order of operations as TerraME: (water + input) - output
+        self.water_in_lake = self.water_in_lake + self._total_input() - self._total_output()
+        self.level = _water_elevation.value(self.water_in_lake)
 
 
 __all__ = ["MonoLake"]

src/dissmodel_sysdyn/models/predatorprey.py

@@ -116,12 +116,14 @@ def execute(self) -> None:
         Updates prey and predator populations according to the
         discrete Lotka-Volterra equations.
         """
-        self.prey += (
-            self.prey_growth * self.prey
+        # same order of operations as TerraME, so that the results match it
+        # to the last bit; the predators use the prey already updated
+        self.prey = (
+            self.prey + self.prey_growth * self.prey
             - self.prey_death_pred * self.prey * self.predator
         )
-        self.predator += (
-            -self.pred_death * self.predator
+        self.predator = (
+            self.predator - self.pred_death * self.predator
             + self.pred_growth_kills * self.prey * self.predator
         )
 

src/dissmodel_sysdyn/models/room_temperature.py

@@ -92,7 +92,7 @@ def execute(self) -> None:
         self.outside = climate(t)
         inflow  = self.thermal_inertia * (self.temp_set - self.inside)
         outflow = self.loss_to_outside  * (self.inside   - self.outside)
-        self.inside += inflow - outflow
+        self.inside = self.inside + inflow - outflow  # TerraME's order of operations
 
 
 __all__ = ["RoomTemperature"]

src/dissmodel_sysdyn/models/sir.py

@@ -119,15 +119,16 @@ def execute(self) -> None:
         Computes new infections and recoveries, then updates the
         susceptible, infected, and recovered compartments.
         """
-        total = self.susceptible + self.infected + self.recovered
-        alpha = self.contacts * self.probability
+        # same order of operations as TerraME, so that the results match it
+        # to the last bit
+        proportion = self.susceptible / (self.susceptible + self.infected + self.recovered)
 
-        new_infected  = self.infected * alpha * (self.susceptible / total)
+        new_infected  = self.infected * self.contacts * self.probability * proportion
         new_recovered = self.infected / self.duration
 
-        self.susceptible -= new_infected
-        self.infected    += new_infected - new_recovered
-        self.recovered   += new_recovered
+        self.susceptible = self.susceptible - new_infected
+        self.recovered   = self.recovered + new_recovered
+        self.infected    = self.infected + new_infected - new_recovered
 
 
 __all__ = ["SIR"]

tests/goldens/README.md

@@ -0,0 +1,26 @@
+# TerraME reference values for the sysdyn models
+
+The `sysdyn` suite of
+[`terrame-goldens`](https://github.com/LambdaGeo/terrame-goldens) v0.1.0:
+what TerraME 2.0.1 computed running the models of the
+[`sysdyn`](https://github.com/TerraME/sysdyn) package (commit `9e671cb`; `sci`
+commit `42b0b3e` for MonoLake's spline).
+
+- `<case>.csv.gz`: a column `step` and one column per variable, from step 0
+  (the state after `init`) to `finalTime`, written with 17 significant digits,
+  so every double is recovered exactly. A variable that does not exist yet
+  at step 0 is empty. Lorenz keeps steps 0 to 100 and then every 100th step.
+- `random_walk.csv.gz`: `seed, step, value` for seeds 1 to 50.
+- `cases.csv`: model, arguments, final time and variables of each case.
+- `reference.json`: the source version, the package commits and the SHA-256
+  of each source file.
+
+To copy them 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 `sysdyn` and `sci` packages
+(INPE and TerraLAB/UFOP); they are offered here under CC-BY-4.0, as the
+goldens they come from.

tests/goldens/cases.csv

@@ -0,0 +1,18 @@
+case,model,arguments,seeds,final_time,variables
+chaotic_growth,ChaoticGrowth,,,300,pop
+coffee,Coffee,,,20,temperature
+daisyworld,Daisyworld,,,100,whiteArea;blackArea;emptyArea;aveTemp;tempNearWhite;tempNearBlack
+homeostasis,Homeostasis,,,30,stock
+homeostasis_100,Homeostasis,stock=100.0,,30,stock
+limited_growth,LimitedGrowth,,,100,pop
+lorenz,Lorenz,,,10000,x;y;z
+mono_lake,MonoLake,,,50,waterInLake;level
+mono_lake_export_50,MonoLake,export=50,,50,waterInLake;level
+population_growth,PopulationGrowth,,,7,population;growth
+predator_prey,PredatorPrey,,,500,prey;predator
+room_temperature,RoomTemperature,,,24,inside;outside
+sir,SIR,,,30,susceptible;infected;recovered;contacts
+tub,Tub,,,8,water
+yeast,Yeast,,,9,cells
+yeast_rate_2,Yeast,rate=2.0,,9,cells
+random_walk,RandomWalk,,50,100,value