Add seed specification from python

.pre-commit-config.yaml

@@ -12,6 +12,8 @@ repos:
       - id: end-of-file-fixer
       - id: trailing-whitespace
       - id: mixed-line-ending
+      - id: name-tests-test
+        args: [--pytest-test-first]
   - repo: https://github.com/crate-ci/typos
     rev: typos-dict-v0.11.20
     hooks:

CHANGELOG.md

@@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+
+* It is now possible to specify a random seed from Python to run the simulator.
+
 ### Changed
 
 * Upgraded PyO3 to v0.21.2.

pyproject.toml

@@ -30,6 +30,7 @@ ignore = [
 ]
 
 [tool.ruff.lint.per-file-ignores]
+"tests/test_*.py" = ["D", "S101", "PLR2004", "ANN201"]
 "examples/*" = ["INP001"]
 "examples/sir.py" = ["T201"]
 "python/rebop/__init__.py" = ["D104"]

python/rebop/__init__.py

@@ -1,3 +1,5 @@
+from __future__ import annotations
+
 import xarray as xr
 
 from .rebop import Gillespie, __version__
@@ -12,13 +14,14 @@ def run_xarray(
     init: dict[str, int],
     tmax: float,
     nb_steps: int,
+    seed: int | None = None,
 ) -> xr.Dataset:
     """Run the system until `tmax` with `nb_steps` steps.
 
     The initial configuration is specified in the dictionary `init`.
     Returns an xarray Dataset.
     """
-    times, result = og_run(self, init, tmax, nb_steps)
+    times, result = og_run(self, init, tmax, nb_steps, seed)
     ds = xr.Dataset(
         data_vars={
             name: xr.DataArray(values, dims="time", coords={"time": times})

src/lib.rs

@@ -299,20 +299,25 @@ impl Gillespie {
     /// The initial configuration is specified in the dictionary `init`.
     /// Returns `times, vars` where `times` is an array of `nb_steps + 1` uniformly spaced time
     /// points between `0` and `tmax`, and `vars` is a dictionary of species name to array of
-    /// values at the given time points.
+    /// values at the given time points.  One can specify a random `seed` for reproducibility.
     fn run(
         &self,
         init: HashMap<String, usize>,
         tmax: f64,
         nb_steps: usize,
+        seed: Option<u64>,
     ) -> PyResult<(Vec<f64>, HashMap<String, Vec<isize>>)> {
         let mut x0 = vec![0; self.species.len()];
         for (name, &value) in &init {
             if let Some(&id) = self.species.get(name) {
                 x0[id] = value as isize;
             }
         }
-        let mut g = gillespie::Gillespie::new(x0);
+        let mut g = match seed {
+            Some(seed) => gillespie::Gillespie::new_with_seed(x0, seed),
+            None => gillespie::Gillespie::new(x0),
+        };
+
         for (rate, reactants, products) in self.reactions.iter() {
             let mut vreactants = vec![0; self.species.len()];
             for reactant in reactants {

tests/__init__.py

@@ -0,0 +1 @@
+"""Tests for the python bindings."""

tests/test_rebop.py

@@ -0,0 +1,36 @@
+import numpy as np
+import numpy.testing as npt
+import pytest
+import rebop
+import xarray as xr
+
+
+def sir_model(transmission: float = 1e-4, recovery: float = 0.01) -> rebop.Gillespie:
+    sir = rebop.Gillespie()
+    sir.add_reaction(transmission, ["S", "I"], ["I", "I"])
+    sir.add_reaction(recovery, ["I"], ["R"])
+    return sir
+
+
+@pytest.mark.parametrize("seed", [None, *range(10)])
+def test_sir(seed: int):
+    sir = sir_model()
+    ds = sir.run({"S": 999, "I": 1}, tmax=250, nb_steps=250, seed=seed)
+    assert isinstance(ds, xr.Dataset)
+    npt.assert_array_equal(ds.time, np.arange(251))
+    assert all(ds.S >= 0)
+    assert all(ds.I >= 0)
+    assert all(ds.R >= 0)
+    assert all(ds.S <= 1000)
+    assert all(ds.I <= 1000)
+    assert all(ds.R <= 1000)
+    npt.assert_array_equal(ds.S + ds.I + ds.R, [1000] * 251)
+
+
+def test_fixed_seed():
+    sir = sir_model()
+    ds = sir.run({"S": 999, "I": 1}, tmax=250, nb_steps=250, seed=42)
+
+    assert ds.S[-1] == 0
+    assert ds.I[-1] == 166
+    assert ds.R[-1] == 834