Parameters

python/rebop/gillespie.py

@@ -2,6 +2,7 @@
 
 from __future__ import annotations
 
+import warnings
 from collections.abc import Sequence
 from typing import TypeAlias
 
@@ -35,6 +36,27 @@ def add_reaction(
     ) -> None:
         """Add a reaction to the system.
 
+        Reaction rates can be specified with a number (for a reaction obeying
+        the Law of Mass Action) or a string (for an arbitrary reaction rate).
+        The string can involve any species involved in reactions, and also
+        parameters defined with the `params` argument of the `run` method.
+
+        If you can, use the law of mass action, which is probably going to be
+        more efficient and can be more correct. For example, for a dimerisation
+        equation:
+
+        ```python
+        s = Gillespie()
+        # Correct, and recommended
+        s.add_reaction(4.2, ["A", "A"], ["AA"])
+        # Correct, but not recommended
+        s.add_reaction("4.2 * A * (A-1)", ["A", "A"], ["AA"])
+        # Incorrect (this would be a constant propensity)
+        s.add_reaction("4.2", ["A", "A"], ["AA"])
+        # Incorrect (incorrect expression)
+        s.add_reaction("4.2 * A^2", ["A", "A"], ["AA"])
+        ```
+
         Example
         -------
         ```python
@@ -45,7 +67,11 @@ def add_reaction(
         # with forward rate 0.1 and reverse rate 0.01
         s.add_reaction(0.1, ["A", "B"], ["C"], 0.01)
         # Add the transformation B -> C with a non-LMA rate
-        s.add_reaction("2.1 * B * C / (5 + C)", ["B"], ["C"])
+        s.add_reaction("2.1 * B * A / (5 + A)", ["B"], ["C"])
+        # Add the transformation B -> C with a non-LMA rate and parameters
+        # that need to be defined by passing `params={"k": 2.1, "K_A": 5}` to
+        # the method `run`
+        s.add_reaction("k * B * A / (K_A + A)", ["B"], ["C"])
         ```
 
         Parameters
@@ -74,6 +100,7 @@ def run(  # noqa: PLR0913 too many parameters in function definition
         tmax: float,
         nb_steps: int,
         *,
+        params: dict[str, float] | None = None,
         rng: RNGLike | SeedLike | None = None,
         sparse: bool | None = None,
         var_names: Sequence[str] | None = None,
@@ -90,6 +117,8 @@ def run(  # noqa: PLR0913 too many parameters in function definition
         nb_steps : int
             Number of steps to return, equally distributed between 0 and `tmax`.
             If 0, then all reactions are returned.
+        params : dict[str, float] | None
+            Dictionary of values for the parameters that appear in the rates.
         rng : RNGLike | SeedLike | None
             Numpy `Generator`, `BitGenerator` or seed.
         sparse : bool | None
@@ -107,12 +136,17 @@ def run(  # noqa: PLR0913 too many parameters in function definition
             being chemical species (all of them, or just the subset defined
             by `var_names`).
         """
+        params = {} if params is None else params
         rng_ = np.random.default_rng(rng)
         seed = rng_.integers(np.iinfo(np.uint64).max, dtype=np.uint64)
+        try:
+            self.gillespie.set_init(init)
+        except UserWarning as e:
+            warnings.warn(e, stacklevel=2)
         times, result = self.gillespie.run(
-            init,
             tmax,
             nb_steps,
+            params=params,
             seed=seed,
             sparse=sparse,
             var_names=var_names,

src/expr.rs

@@ -42,31 +42,44 @@ pub(crate) enum PExpr {
 }
 
 impl PExpr {
-    pub fn to_expr(&self, species: &HashMap<String, usize>) -> Result<Expr, String> {
+    pub fn to_expr(
+        &self,
+        species: &HashMap<String, usize>,
+        params: &HashMap<String, f64>,
+    ) -> Result<Expr, String> {
         let expr = match self {
             PExpr::Constant(c) => Expr::Constant(*c),
-            PExpr::Variable(s) => Expr::Concentration(
-                *species.get(s).ok_or(
-                    format!("species `{s}` is in a rate but is not involved in any reaction")
-                        .to_string(),
-                )?,
-            ),
-            PExpr::Add(a, b) => {
-                Expr::Add(Box::new(a.to_expr(species)?), Box::new(b.to_expr(species)?))
-            }
-            PExpr::Sub(a, b) => {
-                Expr::Sub(Box::new(a.to_expr(species)?), Box::new(b.to_expr(species)?))
-            }
-            PExpr::Mul(a, b) => {
-                Expr::Mul(Box::new(a.to_expr(species)?), Box::new(b.to_expr(species)?))
+            PExpr::Variable(s) => {
+                if let Some(i) = species.get(s) {
+                    Expr::Concentration(*i)
+                } else {
+                    let v = params
+                        .get(s)
+                        .ok_or(format!("Parameter {s} should have a value").to_string())?;
+                    Expr::Constant(*v)
+                }
             }
-            PExpr::Div(a, b) => {
-                Expr::Div(Box::new(a.to_expr(species)?), Box::new(b.to_expr(species)?))
-            }
-            PExpr::Pow(a, b) => {
-                Expr::Pow(Box::new(a.to_expr(species)?), Box::new(b.to_expr(species)?))
-            }
-            PExpr::Exp(a) => Expr::Exp(Box::new(a.to_expr(species)?)),
+            PExpr::Add(a, b) => Expr::Add(
+                Box::new(a.to_expr(species, params)?),
+                Box::new(b.to_expr(species, params)?),
+            ),
+            PExpr::Sub(a, b) => Expr::Sub(
+                Box::new(a.to_expr(species, params)?),
+                Box::new(b.to_expr(species, params)?),
+            ),
+            PExpr::Mul(a, b) => Expr::Mul(
+                Box::new(a.to_expr(species, params)?),
+                Box::new(b.to_expr(species, params)?),
+            ),
+            PExpr::Div(a, b) => Expr::Div(
+                Box::new(a.to_expr(species, params)?),
+                Box::new(b.to_expr(species, params)?),
+            ),
+            PExpr::Pow(a, b) => Expr::Pow(
+                Box::new(a.to_expr(species, params)?),
+                Box::new(b.to_expr(species, params)?),
+            ),
+            PExpr::Exp(a) => Expr::Exp(Box::new(a.to_expr(species, params)?)),
         };
         Ok(expr)
     }
@@ -321,27 +334,32 @@ mod tests {
             )),
         );
         // fmt: on
-        let mut h = HashMap::new();
-        h.insert("A".to_string(), 0);
-        h.insert("B".to_string(), 1);
-        h.insert("C".to_string(), 2);
-        h.insert("D".to_string(), 3);
-        h.insert("E".to_string(), 4);
-        h.insert("F".to_string(), 5);
-        assert_eq!(pe.to_expr(&h), Ok(e));
+        let mut init = HashMap::new();
+        let mut params = HashMap::new();
+        init.insert("A".to_string(), 0);
+        init.insert("B".to_string(), 1);
+        init.insert("C".to_string(), 2);
+        init.insert("D".to_string(), 3);
+        init.insert("E".to_string(), 4);
+        init.insert("F".to_string(), 5);
+        assert_eq!(pe.to_expr(&init, &params), Ok(e));
     }
 
     #[test]
     fn test_eval() {
         let pe: PExpr = "1.21 * C + B - A / D ^ E * (F + exp(D))".parse().unwrap();
-        let mut h = HashMap::new();
-        h.insert("A".to_string(), 0);
-        h.insert("B".to_string(), 1);
-        h.insert("C".to_string(), 2);
-        h.insert("D".to_string(), 3);
-        h.insert("E".to_string(), 4);
-        h.insert("F".to_string(), 5);
+        let mut init = HashMap::new();
+        let mut params = HashMap::new();
+        init.insert("A".to_string(), 0);
+        init.insert("B".to_string(), 1);
+        init.insert("C".to_string(), 2);
+        init.insert("D".to_string(), 3);
+        init.insert("E".to_string(), 4);
+        init.insert("F".to_string(), 5);
         let mut species = vec![2, 3, 5, 7, 11, 13];
-        assert_eq!(pe.to_expr(&h).unwrap().eval(&species), 9.049998877643098);
+        assert_eq!(
+            pe.to_expr(&init, &params).unwrap().eval(&species),
+            9.049998877643098
+        );
     }
 }

src/lib.rs

@@ -230,7 +230,7 @@
 //! * [SmartCell](http://software.crg.es/smartcell/)
 //! * [NFsim](http://michaelsneddon.net/nfsim/)
 
-use pyo3::exceptions::{PyKeyError, PyValueError};
+use pyo3::exceptions::{PyUserWarning, PyValueError};
 use pyo3::prelude::*;
 use std::collections::HashMap;
 use std::fmt::{self, Display, Formatter};
@@ -248,6 +248,7 @@ use expr::PExpr;
 #[pyclass]
 struct Gillespie {
     species: HashMap<String, usize>,
+    init: HashMap<String, usize>,
     reactions: Vec<(PRate, Vec<String>, Vec<String>)>,
 }
 
@@ -282,6 +283,7 @@ impl PRate {
     fn to_gillespie_rate(
         &self,
         species: &HashMap<String, usize>,
+        params: &HashMap<String, f64>,
     ) -> Result<gillespie::Rate, String> {
         let rate = match self {
             PRate::Lma(rate, reactants) => {
@@ -291,7 +293,7 @@ impl PRate {
                 }
                 gillespie::Rate::lma(*rate, rate_reactants)
             }
-            PRate::Expr(nomexpr) => gillespie::Rate::expr(nomexpr.to_expr(species)?),
+            PRate::Expr(pexpr) => gillespie::Rate::expr(pexpr.to_expr(species, params)?),
         };
         Ok(rate)
     }
@@ -323,6 +325,7 @@ impl Gillespie {
     fn new() -> Self {
         Gillespie {
             species: HashMap::new(),
+            init: HashMap::new(),
             reactions: Vec::new(),
         }
     }
@@ -369,6 +372,24 @@ impl Gillespie {
     fn nb_reactions(&self) -> PyResult<usize> {
         Ok(self.reactions.len())
     }
+    /// Set the initial count of species
+    fn set_init(&mut self, init: HashMap<String, usize>) -> PyResult<()> {
+        let mut warning = false;
+        for species in init.keys() {
+            if !self.species.contains_key(species) {
+                warning = true;
+                self.add_species(species);
+            }
+        }
+        self.init = init;
+        if warning {
+            let message =
+                "Some species are not involved in any reactions. You should probably instead use parameters.";
+            Err(PyUserWarning::new_err(message))
+        } else {
+            Ok(())
+        }
+    }
     /// Run the system until `tmax` with `nb_steps` steps.
     ///
     /// The initial configuration is specified in the dictionary `init`.
@@ -377,19 +398,26 @@ impl Gillespie {
     /// values at the given time points.  One can specify a random `seed` for reproducibility.
     /// If `nb_steps` is `0`, then returns all reactions, ending with the first that happens at
     /// or after `tmax`.
-    #[pyo3(signature = (init, tmax, nb_steps, seed=None, sparse=None, var_names=None))]
+    #[pyo3(signature = (tmax, nb_steps, params, seed=None, sparse=None, var_names=None))]
     fn run(
         &self,
-        init: HashMap<String, usize>,
         tmax: f64,
         nb_steps: usize,
+        params: HashMap<String, f64>,
         seed: Option<u64>,
         sparse: Option<bool>,
         var_names: Option<Vec<String>>,
     ) -> PyResult<(Vec<f64>, HashMap<String, Vec<isize>>)> {
+        for p in params.keys() {
+            if self.species.contains_key(p) {
+                return Err(PyValueError::new_err(format!(
+                    "Species {p} cannot also be a parameter."
+                )));
+            }
+        }
         let sparse = sparse.unwrap_or(false);
         let mut x0 = vec![0; self.species.len()];
-        for (name, &value) in &init {
+        for (name, &value) in &self.init {
             if let Some(&id) = self.species.get(name) {
                 x0[id] = value as isize;
             }
@@ -418,9 +446,9 @@ impl Gillespie {
             for product in products {
                 actions[self.species[product]] += 1;
             }
-            match rate.to_gillespie_rate(&self.species) {
+            match rate.to_gillespie_rate(&self.species, &params) {
                 Ok(rate) => g.add_reaction(rate, actions),
-                Err(e) => return Err(PyKeyError::new_err(e)),
+                Err(e) => return Err(PyValueError::new_err(e)),
             }
         }
         let mut times = Vec::new();
@@ -463,13 +491,6 @@ impl Gillespie {
                 }
             }
         }
-        // Add the species that have been given in the initial condition
-        // but are not involved in the reactions
-        for (species, value) in init.iter() {
-            if !self.species.contains_key(species) {
-                result.insert(species.clone(), vec![*value as isize; nb_steps + 1]);
-            }
-        }
         Ok((times, result))
     }
     fn __str__(&self) -> PyResult<String> {

tests/test_rebop.py

@@ -98,17 +98,13 @@ def test_arbitrary_rates() -> None:
 
     s = rebop.Gillespie()
     s.add_reaction("B", [], ["A"])
-    with pytest.raises(
-        KeyError, match="species `B` is in a rate but is not involved in any reaction"
+    with pytest.warns(
+        UserWarning,
+        match="species are not involved in any reaction",
     ):
-        s.run({}, tmax=10, nb_steps=100)
-
-    s = rebop.Gillespie()
-    s.add_reaction("B", [], ["A"])
-    with pytest.raises(
-        KeyError, match="species `B` is in a rate but is not involved in any reaction"
-    ):
-        s.run({"B": 1}, tmax=10, nb_steps=100)
+        ds = s.run({"B": 1}, tmax=10, nb_steps=100)
+    assert ds.A[-1] >= 1
+    np.testing.assert_equal(ds.B.values, [1] * 101)
 
     s = rebop.Gillespie()
     s.add_reaction("B", [], ["A"])
@@ -154,3 +150,17 @@ def test_run_empty() -> None:
     ds = s.run({}, tmax=10, nb_steps=10)
     expected = xr.Dataset(coords={"time": np.linspace(0, 10, 11)})
     xr.testing.assert_identical(ds, expected)
+
+
+def test_parameters() -> None:
+    s = rebop.Gillespie()
+    s.add_reaction(4, ["A"], ["B"])
+    with pytest.raises(ValueError, match="Species B cannot also be a parameter"):
+        s.run({}, 10, 10, params={"B": 4.2})
+
+    s = rebop.Gillespie()
+    s.add_reaction("k", [], ["A"])
+    with pytest.raises(ValueError, match="Parameter k should have a value"):
+        s.run({}, 10, 10)
+    ds = s.run({}, 10, 10, params={"k": 0.4})
+    assert ds.A[-1] > 0