From ec19439ab670ce463d75868fa0896aae14774a8f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9rgio=20Souza=20Costa?= Date: Tue, 20 Oct 2026 20:10:46 -0300 Subject: [PATCH 1/3] fix(tub): stop overwriting Model._step with the inflow counter Tub kept its execution counter in ``self._step``, the attribute where dissmodel's Model stores the interval between ticks. Each increment widened the interval, so a run from t=1 to 30 executed only at ticks 1, 2, 4, 7, 11, 16, 22 and 29 (8 executions instead of 30), and the periodic inflow never matched the TerraME Tub example. The counter is now ``_tick``. A new test checks, for every exported model, that it executes exactly once per tick; it fails on the old Tub and passes on all models now. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_011GTEeQcZQHHFJCk7Jdqe61 --- CHANGELOG.md | 7 +++++++ src/dissmodel_sysdyn/models/tub.py | 11 +++++++---- tests/test_models_smoke.py | 23 +++++++++++++++++++++++ 3 files changed, 37 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ee7b0f5..809109f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,12 @@ # Changelog +## [Unreleased] + +- **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 + `_tick`; a test checks that every model executes once per tick. + ## [0.1.1] — 2026-10-08 First citable release. The `v0.1.0` tag was created before the release diff --git a/src/dissmodel_sysdyn/models/tub.py b/src/dissmodel_sysdyn/models/tub.py index b2e4007..53c6273 100644 --- a/src/dissmodel_sysdyn/models/tub.py +++ b/src/dissmodel_sysdyn/models/tub.py @@ -39,7 +39,7 @@ class Tub(Model): Examples -------- >>> from dissmodel.core import Environment - >>> env = Environment(end_time=8) + >>> env = Environment(start_time=1, end_time=30) >>> tub = Tub(water=40, out_flow=5, in_flow=20, in_period=10) >>> env.run() """ @@ -81,7 +81,10 @@ def setup( self.out_flow = out_flow self.in_flow = in_flow self.in_period = in_period - self._step = 0 + # Executions counter. Must not be named ``_step``: dissmodel's Model + # keeps the scheduling interval in ``_step``, so incrementing it here + # made the interval grow (ticks 1, 2, 4, 7, 11, ...). + self._tick = 0 def execute(self) -> None: """ @@ -90,13 +93,13 @@ def execute(self) -> None: Applies the outflow every step and the inflow every ``in_period`` steps. Water is clamped to zero from below. """ - self._step += 1 + self._tick += 1 self.water -= self.out_flow if self.water < 0.0: self.water = 0.0 - if self._step % self.in_period == 0: + if self._tick % self.in_period == 0: self.water += self.in_flow diff --git a/tests/test_models_smoke.py b/tests/test_models_smoke.py index c06657b..09ace54 100644 --- a/tests/test_models_smoke.py +++ b/tests/test_models_smoke.py @@ -44,3 +44,26 @@ def test_model_runs(name, cls): assert numeric, f"{name} exposes no numeric state" for key, value in numeric.items(): assert math.isfinite(value), f"{name}.{key} = {value}" + + +@pytest.mark.parametrize("name, cls", MODELS, ids=[n for n, _ in MODELS]) +def test_model_executes_once_per_tick(name, cls, monkeypatch): + """A model must not change its own scheduling. + + ``dissmodel.core.Model`` keeps the interval between ticks in ``_step``; + a model that reuses that name (as ``Tub`` did up to 0.1.1) silently + skips ticks. + """ + ticks: list[float] = [] + original = cls.execute + + def counting(self): + ticks.append(self.env.now()) + original(self) + + monkeypatch.setattr(cls, "execute", counting) + np.random.seed(0) + env = Environment(start_time=1, end_time=STEPS) + cls() + env.run() + assert ticks == list(range(1, STEPS + 1)) From deac8dae4a41cc04be9447dae5bff87a33568dfe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9rgio=20Souza=20Costa?= Date: Tue, 20 Oct 2026 20:22:24 -0300 Subject: [PATCH 2/3] docs(examples): run from t=1 to finalTime, as TerraME does The examples used Environment() with env.run(N), or start_time=0, which runs N + 1 steps (t = 0..N, end_time is inclusive). A TerraME Timer with finalTime = N runs N steps at t = 1..N. The difference is one extra step in every model and, in RoomTemperature, a one-hour shift of the outside temperature curve (up to 2.3 degrees at the same step). Notebooks, CLI and Streamlit examples, the demo app and the docs now use Environment(start_time=1, end_time=N) and env.run(). docs/models.md explains the convention. All notebooks execute (no stored outputs). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_011GTEeQcZQHHFJCk7Jdqe61 --- CHANGELOG.md | 4 +++ demo/app.py | 2 +- docs/architecture.md | 2 +- docs/models.md | 11 ++++++-- docs/reference.md | 28 +++++++++---------- examples/cli/sysdyn_sir.py | 4 +-- .../notebooks/sysdyn_chaotic_growth.ipynb | 4 +-- examples/notebooks/sysdyn_coffee.ipynb | 4 +-- examples/notebooks/sysdyn_daisyworld.ipynb | 4 +-- examples/notebooks/sysdyn_homeostasis.ipynb | 4 +-- .../notebooks/sysdyn_limited_growth.ipynb | 4 +-- examples/notebooks/sysdyn_lorenz.ipynb | 4 +-- examples/notebooks/sysdyn_mono_lake.ipynb | 4 +-- .../notebooks/sysdyn_population_growth.ipynb | 4 +-- examples/notebooks/sysdyn_predator_prey.ipynb | 4 +-- examples/notebooks/sysdyn_random_walk.ipynb | 4 +-- .../notebooks/sysdyn_room_temperature.ipynb | 4 +-- examples/notebooks/sysdyn_sir.ipynb | 4 +-- examples/notebooks/sysdyn_tub.ipynb | 4 +-- examples/notebooks/sysdyn_yeast.ipynb | 4 +-- examples/streamlit/sysdyn_all.py | 2 +- examples/streamlit/sysdyn_sir.py | 2 +- 22 files changed, 61 insertions(+), 50 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 809109f..f69872b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,10 @@ `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 `_tick`; a test checks that every model executes once per tick. +- Examples, notebooks, demo and docs use `Environment(start_time=1, end_time=N)`, + which runs N steps at t = 1..N like a TerraME `Timer` with `finalTime = N`. + The default `start_time=0` runs one step more and shifts the clock by one hour + in `RoomTemperature`. A note in the docs explains the convention. ## [0.1.1] — 2026-10-08 diff --git a/demo/app.py b/demo/app.py index 17b6e81..98fdf51 100644 --- a/demo/app.py +++ b/demo/app.py @@ -64,7 +64,7 @@ # --------------------------------------------------------------------------- # 1. Environment -env = Environment(start_time=0, end_time=steps) +env = Environment(start_time=1, end_time=steps) # 2. Model model = ModelClass() diff --git a/docs/architecture.md b/docs/architecture.md index 84094ff..3ae0317 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -76,7 +76,7 @@ class Decay(Model): self.amount -= decay # Running the simulation -env = Environment(end_time=50) +env = Environment(start_time=1, end_time=50) model = Decay() env.run() ``` diff --git a/docs/models.md b/docs/models.md index b58e75a..b8fe103 100644 --- a/docs/models.md +++ b/docs/models.md @@ -17,9 +17,16 @@ from dissmodel_sysdyn.models import SIR import matplotlib.pyplot as plt # Create the simulation environment -env = Environment(end_time=60) +env = Environment(start_time=1, end_time=60) ``` +!!! note "Counting steps as in TerraME" + `end_time` is inclusive. With `start_time=1, end_time=N` the model runs + N steps at t = 1, ..., N, like a TerraME `Timer` with `finalTime = N`. + The `Environment` default, `start_time=0`, runs one step more (t = 0, ..., N) + and shifts the clock by one, which matters for models that read `env.now()`, + such as `RoomTemperature`. + ### Step 2: Model Configuration We instantiate the SIR model. We start with a nearly susceptible population and a small number of "Patient Zeros". @@ -76,7 +83,7 @@ from dissmodel.core import Environment from dissmodel_sysdyn.models import PredatorPrey import matplotlib.pyplot as plt -env = Environment(end_time=200) +env = Environment(start_time=1, end_time=200) ``` ### Step 2: Define Biological Parameters diff --git a/docs/reference.md b/docs/reference.md index 50c308c..333e4ad 100644 --- a/docs/reference.md +++ b/docs/reference.md @@ -29,7 +29,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import ChaoticGrowth -env = Environment(end_time=100) +env = Environment(start_time=1, end_time=100) model = ChaoticGrowth(pop=0.1, rate=3.9) env.run() ``` @@ -62,7 +62,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import Coffee -env = Environment(end_time=30) +env = Environment(start_time=1, end_time=30) coffee = Coffee(temperature=90, room_temperature=22, cooling_rate=0.15) env.run() ``` @@ -105,7 +105,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import Daisyworld -env = Environment(end_time=200) +env = Environment(start_time=1, end_time=200) model = Daisyworld(sun_luminosity=1.0) env.run() ``` @@ -138,7 +138,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import Homeostasis -env = Environment(end_time=50) +env = Environment(start_time=1, end_time=50) model = Homeostasis(stock=10, gain=5, rate=-0.2) env.run() ``` @@ -171,7 +171,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import LimitedGrowth -env = Environment(end_time=100) +env = Environment(start_time=1, end_time=100) model = LimitedGrowth(pop=10, rate=0.2, capacity=1000) env.run() ``` @@ -211,7 +211,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import Lorenz -env = Environment(end_time=5000) +env = Environment(start_time=1, end_time=5000) model = Lorenz() env.run() ``` @@ -245,7 +245,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import MonoLake -env = Environment(end_time=50) +env = Environment(start_time=1, end_time=50) model = MonoLake(export=150) env.run() ``` @@ -279,7 +279,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import PopulationGrowth -env = Environment(end_time=20) +env = Environment(start_time=1, end_time=20) model = PopulationGrowth(population=100, growth=0.05, growth_change=1.1) env.run() ``` @@ -311,7 +311,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import PredatorPrey -env = Environment(end_time=200) +env = Environment(start_time=1, end_time=200) model = PredatorPrey() env.run() ``` @@ -347,7 +347,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import RandomWalk -env = Environment(end_time=100) +env = Environment(start_time=1, end_time=100) model = RandomWalk(prob=0.55, seed=42) env.run() ``` @@ -383,7 +383,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import RoomTemperature -env = Environment(end_time=24) +env = Environment(start_time=1, end_time=24) model = RoomTemperature() env.run() ``` @@ -420,7 +420,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import SIR -env = Environment(end_time=60) +env = Environment(start_time=1, end_time=60) sir = SIR(infected=5, probability=0.3) env.run() ``` @@ -454,7 +454,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import Tub -env = Environment(end_time=50) +env = Environment(start_time=1, end_time=50) tub = Tub(water=100, out_flow=2, in_flow=50, in_period=20) env.run() ``` @@ -486,7 +486,7 @@ def execute(self) -> None: from dissmodel.core import Environment from dissmodel_sysdyn.models import Yeast -env = Environment(end_time=15) +env = Environment(start_time=1, end_time=15) yeast = Yeast() env.run() ``` diff --git a/examples/cli/sysdyn_sir.py b/examples/cli/sysdyn_sir.py index 4ff160d..3845df4 100644 --- a/examples/cli/sysdyn_sir.py +++ b/examples/cli/sysdyn_sir.py @@ -26,7 +26,7 @@ # --------------------------------------------------------------------------- # 1. Environment -env = Environment() +env = Environment(start_time=1, end_time=30) # 2. Model SIR( @@ -44,4 +44,4 @@ # --------------------------------------------------------------------------- # Run # --------------------------------------------------------------------------- -env.run(30) +env.run() diff --git a/examples/notebooks/sysdyn_chaotic_growth.ipynb b/examples/notebooks/sysdyn_chaotic_growth.ipynb index fc0e85b..05d9cb8 100644 --- a/examples/notebooks/sysdyn_chaotic_growth.ipynb +++ b/examples/notebooks/sysdyn_chaotic_growth.ipynb @@ -65,7 +65,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=100)" ] }, { @@ -147,7 +147,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(100)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_coffee.ipynb b/examples/notebooks/sysdyn_coffee.ipynb index 1d04a2c..2fba009 100644 --- a/examples/notebooks/sysdyn_coffee.ipynb +++ b/examples/notebooks/sysdyn_coffee.ipynb @@ -62,7 +62,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=50)" ] }, { @@ -147,7 +147,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(50)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_daisyworld.ipynb b/examples/notebooks/sysdyn_daisyworld.ipynb index e2e8ac4..5a08911 100644 --- a/examples/notebooks/sysdyn_daisyworld.ipynb +++ b/examples/notebooks/sysdyn_daisyworld.ipynb @@ -65,7 +65,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=100)" ] }, { @@ -178,7 +178,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(100)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_homeostasis.ipynb b/examples/notebooks/sysdyn_homeostasis.ipynb index a557647..a3d26a3 100644 --- a/examples/notebooks/sysdyn_homeostasis.ipynb +++ b/examples/notebooks/sysdyn_homeostasis.ipynb @@ -62,7 +62,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=30)" ] }, { @@ -147,7 +147,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(30)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_limited_growth.ipynb b/examples/notebooks/sysdyn_limited_growth.ipynb index 588e348..3c28b73 100644 --- a/examples/notebooks/sysdyn_limited_growth.ipynb +++ b/examples/notebooks/sysdyn_limited_growth.ipynb @@ -62,7 +62,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=100)" ] }, { @@ -147,7 +147,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(100)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_lorenz.ipynb b/examples/notebooks/sysdyn_lorenz.ipynb index 4e4b173..9c87f2e 100644 --- a/examples/notebooks/sysdyn_lorenz.ipynb +++ b/examples/notebooks/sysdyn_lorenz.ipynb @@ -64,7 +64,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=2000)" ] }, { @@ -150,7 +150,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(2000)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_mono_lake.ipynb b/examples/notebooks/sysdyn_mono_lake.ipynb index d372dbe..3d66b34 100644 --- a/examples/notebooks/sysdyn_mono_lake.ipynb +++ b/examples/notebooks/sysdyn_mono_lake.ipynb @@ -63,7 +63,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=50)" ] }, { @@ -156,7 +156,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(50)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_population_growth.ipynb b/examples/notebooks/sysdyn_population_growth.ipynb index e2d8bcb..b156e25 100644 --- a/examples/notebooks/sysdyn_population_growth.ipynb +++ b/examples/notebooks/sysdyn_population_growth.ipynb @@ -63,7 +63,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=20)" ] }, { @@ -142,7 +142,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(20)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_predator_prey.ipynb b/examples/notebooks/sysdyn_predator_prey.ipynb index faa61d6..9c9c487 100644 --- a/examples/notebooks/sysdyn_predator_prey.ipynb +++ b/examples/notebooks/sysdyn_predator_prey.ipynb @@ -63,7 +63,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=500)" ] }, { @@ -146,7 +146,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(500)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_random_walk.ipynb b/examples/notebooks/sysdyn_random_walk.ipynb index c795c63..cfd230b 100644 --- a/examples/notebooks/sysdyn_random_walk.ipynb +++ b/examples/notebooks/sysdyn_random_walk.ipynb @@ -63,7 +63,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=100)" ] }, { @@ -145,7 +145,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(100)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_room_temperature.ipynb b/examples/notebooks/sysdyn_room_temperature.ipynb index 084736c..97eef2e 100644 --- a/examples/notebooks/sysdyn_room_temperature.ipynb +++ b/examples/notebooks/sysdyn_room_temperature.ipynb @@ -63,7 +63,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=24)" ] }, { @@ -149,7 +149,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(24)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_sir.ipynb b/examples/notebooks/sysdyn_sir.ipynb index e7c8a1c..87d96cc 100644 --- a/examples/notebooks/sysdyn_sir.ipynb +++ b/examples/notebooks/sysdyn_sir.ipynb @@ -68,7 +68,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=30)" ] }, { @@ -156,7 +156,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(30)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_tub.ipynb b/examples/notebooks/sysdyn_tub.ipynb index cced38f..4f0c424 100644 --- a/examples/notebooks/sysdyn_tub.ipynb +++ b/examples/notebooks/sysdyn_tub.ipynb @@ -62,7 +62,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=100)" ] }, { @@ -144,7 +144,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(100)" + "env.run()" ] }, { diff --git a/examples/notebooks/sysdyn_yeast.ipynb b/examples/notebooks/sysdyn_yeast.ipynb index 8c95efb..ba938d7 100644 --- a/examples/notebooks/sysdyn_yeast.ipynb +++ b/examples/notebooks/sysdyn_yeast.ipynb @@ -63,7 +63,7 @@ "metadata": {}, "outputs": [], "source": [ - "env = Environment()" + "env = Environment(start_time=1, end_time=15)" ] }, { @@ -150,7 +150,7 @@ "metadata": {}, "outputs": [], "source": [ - "env.run(15)" + "env.run()" ] }, { diff --git a/examples/streamlit/sysdyn_all.py b/examples/streamlit/sysdyn_all.py index 6f53708..2bd4cad 100644 --- a/examples/streamlit/sysdyn_all.py +++ b/examples/streamlit/sysdyn_all.py @@ -61,7 +61,7 @@ # --------------------------------------------------------------------------- # 1. Environment -env = Environment(start_time=0, end_time=steps) +env = Environment(start_time=1, end_time=steps) # 2. Model — instantiated dynamically from sidebar selection ModelClass = model_classes[model_name] diff --git a/examples/streamlit/sysdyn_sir.py b/examples/streamlit/sysdyn_sir.py index b4b986f..22f95b7 100644 --- a/examples/streamlit/sysdyn_sir.py +++ b/examples/streamlit/sysdyn_sir.py @@ -42,7 +42,7 @@ # --------------------------------------------------------------------------- # 1. Environment -env = Environment(start_time=0, end_time=steps) +env = Environment(start_time=1, end_time=steps) # 2. Model sir = SIR( From 13d951a9a61c772decfc620d14d9f74b50abc889 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9rgio=20Souza=20Costa?= Date: Tue, 20 Oct 2026 20:22:34 -0300 Subject: [PATCH 3/3] test: compare every model with the TerraME sysdyn package tests/test_terrame_parity.py transcribes, step by step, the update rule of each model in terrame/sysdyn (lua/.lua, commit 9e671cb) with the TerraME default parameters and finalTime, and compares the whole trajectory with the DisSModel model run from t=1 to finalTime (rtol = atol = 1e-9). Coffee, ChaoticGrowth, Homeostasis, LimitedGrowth, Lorenz, PopulationGrowth, PredatorPrey, RoomTemperature, SIR, Tub (default and with the periodic inflow) and Yeast match. Two models differ and are marked as strict xfail, with the cause, so that a fix shows up: - MonoLake: TerraME's sci.Spline is a Catmull-Rom curve sampled at 1000 points per segment and read by linear interpolation (transcribed in the test); the model uses scipy CubicSpline. Up to 26 KAF in 50 steps. - Daisyworld: TerraME recomputes the planet albedo for the black daisies after updating the white area. Up to 0.61 of area during the run, same equilibrium at t = 100. RandomWalk is not compared (different random generators). The SIR "maximum" event and the PredatorPrey "period" option are not ported; they are inactive with the default values. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_011GTEeQcZQHHFJCk7Jdqe61 --- CHANGELOG.md | 7 + tests/test_terrame_parity.py | 287 +++++++++++++++++++++++++++++++++++ 2 files changed, 294 insertions(+) create mode 100644 tests/test_terrame_parity.py diff --git a/CHANGELOG.md b/CHANGELOG.md index f69872b..bad9034 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,13 @@ which runs N steps at t = 1..N like a TerraME `Timer` with `finalTime = N`. The default `start_time=0` runs one step more and shifts the clock by one hour in `RoomTemperature`. A note in the docs explains the convention. +- `tests/test_terrame_parity.py` compares the full trajectory of each model + with a step-by-step transcription of the TerraME `sysdyn` package + (`terrame/sysdyn`, commit 9e671cb). Eleven models match within 1e-9. + Two known differences are kept as strict `xfail`: `MonoLake` (TerraME's + sampled Catmull-Rom spline vs scipy `CubicSpline`) and `Daisyworld` + (TerraME recomputes the planet albedo after updating the white area; same + equilibrium, different trajectory). `RandomWalk` is not compared. ## [0.1.1] — 2026-10-08 diff --git a/tests/test_terrame_parity.py b/tests/test_terrame_parity.py new file mode 100644 index 0000000..eeff55d --- /dev/null +++ b/tests/test_terrame_parity.py @@ -0,0 +1,287 @@ +""" +Parity with the TerraME ``sysdyn`` package. + +Each ``ref_*`` function below transcribes, step by step, the update rule of a +model in https://github.com/TerraME/sysdyn (``lua/.lua``, commit +9e671cb, 2020-08-20), with the TerraME default parameters and ``finalTime``. +A TerraME ``Timer`` runs t = 1, ..., finalTime, so the DisSModel model runs in +``Environment(start_time=1, end_time=finalTime)``. + +The test compares the full trajectory, not only the final value. The reference +is a transcription, not output of TerraME itself: it checks that the Python +models follow the Lua equations and the TerraME clock. + +Not covered: ``RandomWalk`` (different random generators); the SIR +``maximum`` event and the PredatorPrey ``period`` option (not ported, and +inactive with the default values). + +Known differences, kept as strict ``xfail`` so that a fix is noticed: + +* ``MonoLake``: TerraME's ``sci.Spline`` is a Catmull-Rom curve sampled at + ``steps`` points per segment and read by linear interpolation; the Python + model uses ``scipy.interpolate.CubicSpline``. Up to 26 KAF of difference in + 50 steps. +* ``Daisyworld``: the Lua model recomputes the planet albedo for the black + daisies after updating the white area; the Python model uses the albedo + from the start of the step. The trajectories differ (up to 0.61 in area) + and reach the same equilibrium. +""" +from __future__ import annotations + +import numpy as np +import pytest + +import dissmodel_sysdyn.models as sd +from dissmodel.core import Environment + +RTOL = 1e-9 +ATOL = 1e-9 + + +def run_dissmodel(cls, final_time, attrs, **params): + """Run ``cls`` for t = 1..final_time and return the state after each step.""" + env = Environment(start_time=1, end_time=final_time) + model = cls(**params) # noqa: F841 (connects to env on creation) + trace = [] + original = cls.execute + + def recording(self): + original(self) + trace.append([getattr(self, a) for a in attrs]) + + # record the state after each execute(), for this run only + try: + cls.execute = recording + env.run() + finally: + cls.execute = original + return np.array(trace, dtype=float) + + +# --------------------------------------------------------------------------- +# Transcriptions of terrame/sysdyn lua/*.lua +# --------------------------------------------------------------------------- + +def ref_coffee(final_time=20): + temperature, room, out = 80.0, 20.0, [] + for _ in range(final_time): + temperature = temperature - (temperature - room) * 0.1 + out.append([temperature]) + return out + + +def ref_chaotic_growth(final_time=300): + pop, rate, out = 0.10, 4.0, [] + for _ in range(final_time): + pop = rate * pop * (1 - pop) + out.append([pop]) + return out + + +def ref_homeostasis(final_time=30): + stock, gain, rate, out = 0.0, 2.0, -0.4, [] + for _ in range(final_time): + stock = stock + gain + rate * stock + out.append([stock]) + return out + + +def ref_limited_growth(final_time=100): + pop, rate, capacity, out = 300.0, 0.105, 20000.0, [] + for _ in range(final_time): + pop = pop + pop * rate * (1 - pop / capacity) + out.append([pop]) + return out + + +def ref_lorenz(final_time=10000): + x = y = z = 1.0 + delta, rho, sigma, beta = 0.01, 28.0, 10.0, 8.0 / 3.0 + out = [] + for _ in range(final_time): + nx = sigma * (y - x) + ny = x * (rho - z) - y + nz = x * y - beta * z + x, y, z = x + delta * nx, y + delta * ny, z + delta * nz + out.append([x, y, z]) + return out + + +def ref_population_growth(final_time=7): + population, growth, growth_change, out = 60.0, 0.5, 1.0, [] + for _ in range(final_time): + population = population * (1 + growth) + growth = growth * growth_change + out.append([population]) + return out + + +def ref_predator_prey(final_time=500): + predator, prey, out = 40.0, 1000.0, [] + for _ in range(final_time): + prey = prey + 0.08 * prey - 0.001 * prey * predator + predator = predator - 0.02 * predator + 0.00002 * prey * predator # uses the new prey + out.append([prey, predator]) + return out + + +def ref_room_temperature(final_time=24): + temp_set, inside, out = 20.0, 15.0, [] + for t in range(1, final_time + 1): + outside = 15 - 0.1 * (t - 12) ** 2 + inflow = 0.33 * (temp_set - inside) + outflow = 0.30 * (inside - outside) + inside = inside + inflow - outflow + out.append([inside, outside]) + return out + + +def ref_sir(final_time=30): + s, i, r = 9998.0, 2.0, 0.0 + duration, contacts, probability = 2, 6, 0.25 + out = [] + for _ in range(final_time): + proportion = s / (s + i + r) + new_infected = i * contacts * probability * proportion + new_recovered = i / duration + s = s - new_infected + r = r + new_recovered + i = i + new_infected - new_recovered + out.append([s, i, r]) + return out + + +def ref_tub(final_time=8, in_flow=0.0): + water, out_flow, out = 40.0, 5.0, [] + for t in range(1, final_time + 1): + water = water - out_flow + if water < 0: + water = 0.0 + if t >= 10 and t % 10 == 0: # Event{start = 10, period = 10} + water = water + in_flow + out.append([water]) + return out + + +def ref_yeast(final_time=9): + cells, capacity, rate, out = 9.6, 665.0, 1.1, [] + for _ in range(final_time): + cells = cells + cells * rate * (1 - cells / capacity) + if cells > capacity: + cells = capacity + out.append([cells]) + return out + + +_SIGMA, _SOLAR_FLUX, _HEAT = 5.67e-08, 3668.0, 2.06e09 + + +def _daisy_growth(temp_k): + temp = temp_k - 273 + return 1 - 0.003265 * (22.5 - temp) ** 2 if 5.0 < temp < 40.0 else 0.0 + + +def ref_daisyworld(final_time=100): + lum, white, black, empty = 0.7, 0.40, 0.273, 0.327 + + def albedo(): + return white * 0.75 + black * 0.25 + empty * 0.50 + + out = [] + for _ in range(final_time): + ave = (lum * _SOLAR_FLUX * (1 - albedo()) / (4 * _SIGMA)) ** 0.25 + t_white = (_HEAT * (albedo() - 0.75) + ave ** 4) ** 0.25 + white = white + white * (_daisy_growth(t_white) * empty - 0.3) + t_black = (_HEAT * (albedo() - 0.25) + ave ** 4) ** 0.25 # albedo with the new white area + black = black + black * (_daisy_growth(t_black) * empty - 0.3) + empty = 1.0 - (black + white) + out.append([white, black, empty]) + return out + + +def _terrame_spline(xs, ys, steps=1000): + """TerraME sci.Spline (commit 42b0b3e): sampled Catmull-Rom + linear interpolation.""" + pts = list(zip(xs, ys)) + count = len(pts) - 1 + values = [] + for i in range(1, count + 1): + if i == 1: + p0, p1, p2, p3 = pts[0], pts[0], pts[1], pts[2] + elif i == count: + p0, p1, p2, p3 = pts[-3], pts[-2], pts[-1], pts[-1] + else: + p0, p1, p2, p3 = pts[i - 2], pts[i - 1], pts[i], pts[i + 1] + t = 0.0 + while t <= 1: + def cr(a, b, c, d): + return 0.5 * (2 * b + (c - a) * t + (2 * a - 5 * b + 4 * c - d) * t * t + + (3 * b - a - 3 * c + d) * t ** 3) + x = cr(p0[0], p1[0], p2[0], p3[0]) + y = cr(p0[1], p1[1], p2[1], p3[1]) + if not values or (values[-1][0] != x and values[-1][1] != y): + values.append((x, y)) + t += 1 / steps + vx, vy = (np.array(v) for v in zip(*values)) + return lambda v: float(np.interp(v, vx, vy)) + + +def ref_mono_lake(final_time=50): + volume = [0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000] + surface = _terrame_spline(volume, [0, 24.7, 35.3, 48.6, 54.3, 57.2, 61.6, 66.0, 69.9]) + elevation = _terrame_spline(volume, [6224, 6335, 6369, 6392, 6412, 6430, 6447, 6463, 6477]) + water, out = 2228.0, [] + for _ in range(final_time): + area = surface(water) + inflow = area * 0.67 + 150 + 47.6 + outflow = 3.75 * area + 100 + 33.6 + water = water + inflow - outflow + out.append([water, elevation(water)]) + return out + + +# --------------------------------------------------------------------------- + +KNOWN_DIFFERENCE = { + "MonoLake": "TerraME uses a sampled Catmull-Rom spline; the model uses scipy CubicSpline", + "Daisyworld": "TerraME recomputes the planet albedo after updating the white area", +} + +CASES = [ + # name, class, finalTime, attributes, reference, parameters + ("Coffee", sd.Coffee, 20, ["temperature"], ref_coffee, {}), + ("ChaoticGrowth", sd.ChaoticGrowth, 300, ["pop"], ref_chaotic_growth, {}), + ("Homeostasis", sd.Homeostasis, 30, ["stock"], ref_homeostasis, {}), + ("LimitedGrowth", sd.LimitedGrowth, 100, ["pop"], ref_limited_growth, {}), + ("Lorenz", sd.Lorenz, 10000, ["x", "y", "z"], ref_lorenz, {}), + ("PopulationGrowth", sd.PopulationGrowth, 7, ["population"], ref_population_growth, {}), + ("PredatorPrey", sd.PredatorPrey, 500, ["prey", "predator"], ref_predator_prey, {}), + ("RoomTemperature", sd.RoomTemperature, 24, ["inside", "outside"], ref_room_temperature, {}), + ("SIR", sd.SIR, 30, ["susceptible", "infected", "recovered"], ref_sir, {}), + ("Tub", sd.Tub, 8, ["water"], ref_tub, {}), + ("Yeast", sd.Yeast, 9, ["cells"], ref_yeast, {}), + ("Daisyworld", sd.Daisyworld, 100, ["white_area", "black_area", "empty_area"], + ref_daisyworld, {}), + ("MonoLake", sd.MonoLake, 50, ["water_in_lake", "level"], ref_mono_lake, {}), +] + + +def _param(case): + name = case[0] + marks = ([pytest.mark.xfail(reason=KNOWN_DIFFERENCE[name], strict=True)] + if name in KNOWN_DIFFERENCE else []) + return pytest.param(*case, id=name, marks=marks) + + +@pytest.mark.parametrize("name, cls, final_time, attrs, reference, params", + [_param(c) for c in CASES]) +def test_trajectory_matches_terrame(name, cls, final_time, attrs, reference, params): + got = run_dissmodel(cls, final_time, attrs, **params) + expected = np.array(reference(final_time), dtype=float) + assert got.shape == expected.shape + np.testing.assert_allclose(got, expected, rtol=RTOL, atol=ATOL) + + +def test_tub_periodic_inflow(): + """Exercises the TerraME Event{start = 10, period = 10} (inactive by default).""" + got = run_dissmodel(sd.Tub, 30, ["water"], in_flow=20.0) + np.testing.assert_allclose(got, np.array(ref_tub(30, in_flow=20.0)), rtol=RTOL, atol=ATOL)