catlog/stdlib/analyses/ode/
lotka_volterra.rs

1//! Lotka-Volterra ODE analysis of models.
2//!
3//! The main entry point for this module is
4//! [`lotka_volterra_analysis`](SignedCoefficientBuilder::lotka_volterra_analysis).
5
6use std::collections::HashMap;
7use std::hash::Hash;
8use std::ops::Add;
9
10use indexmap::IndexMap;
11use itertools::Itertools;
12use nalgebra::{DMatrix, DVector, Scalar};
13use num_traits::{One, Zero};
14
15#[cfg(feature = "serde")]
16use serde::{Deserialize, Serialize};
17#[cfg(feature = "serde-wasm")]
18use tsify::Tsify;
19
20use super::{ODEAnalysis, Parameter, SignedCoefficientBuilder};
21use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem};
22use crate::{
23    dbl::model::DiscreteDblModel,
24    one::QualifiedPath,
25    zero::{QualifiedName, alg::Polynomial, rig::Monomial},
26};
27
28/// Data defining a Lotka-Volterra ODE problem for a model.
29#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
30#[cfg_attr(feature = "serde-wasm", derive(Tsify))]
31#[cfg_attr(
32    feature = "serde-wasm",
33    tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object)
34)]
35pub struct LotkaVolterraProblemData {
36    /// Map from morphism IDs to interaction coefficients (nonnegative reals).
37    #[cfg_attr(feature = "serde", serde(rename = "interactionCoefficients"))]
38    interaction_coeffs: HashMap<QualifiedName, f32>,
39
40    /// Map from object IDs to growth rates (arbitrary real numbers).
41    #[cfg_attr(feature = "serde", serde(rename = "growthRates"))]
42    growth_rates: HashMap<QualifiedName, f32>,
43
44    /// Map from object IDs to initial values (nonnegative reals).
45    #[cfg_attr(feature = "serde", serde(rename = "initialValues"))]
46    initial_values: HashMap<QualifiedName, f32>,
47
48    /// Duration of simulation.
49    duration: f32,
50}
51
52/// Construct a Lotka-Volterra dynamical system.
53///
54/// A system of ODEs that is affine in its *logarithmic* derivative. These are
55/// sometimes called the "generalized Lotka-Volterra equations." For more, see
56/// [Wikipedia](https://en.wikipedia.org/wiki/Generalized_Lotka%E2%80%93Volterra_equation).
57pub fn lotka_volterra_system<Var, Coef>(
58    vars: &[Var],
59    interaction_coeffs: DMatrix<Coef>,
60    growth_rates: DVector<Coef>,
61) -> PolynomialSystem<Var, Coef, u8>
62where
63    Var: Clone + Hash + Ord,
64    Coef: Clone + Add<Output = Coef> + One + Scalar + Zero,
65{
66    let system = PolynomialSystem {
67        components: interaction_coeffs
68            .row_iter()
69            .zip(vars)
70            .zip(&growth_rates)
71            .map(|((row, i), r)| {
72                (
73                    i.clone(),
74                    Polynomial::<_, Coef, _>::generator(i.clone())
75                        * (row
76                            .iter()
77                            .zip(vars)
78                            .map(|(a, j)| (a.clone(), Monomial::generator(j.clone())))
79                            .collect::<Polynomial<_, _, _>>()
80                            + r.clone()),
81                )
82            })
83            .collect(),
84    };
85    system.normalize()
86}
87
88impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
89    /// Lotka-Volterra ODE analysis for a model of a double theory.
90    ///
91    /// The main application we have in mind is the Lotka-Volterra ODE semantics for
92    /// signed graphs described in our [paper on regulatory
93    /// networks](crate::refs::RegNets).
94    pub fn lotka_volterra_analysis(
95        &self,
96        model: &DiscreteDblModel,
97        data: LotkaVolterraProblemData,
98    ) -> ODEAnalysis<NumericalPolynomialSystem<u8>> {
99        let (system, ob_index) = self.lotka_volterra_system(model);
100        let n = ob_index.len();
101
102        let initial_values = ob_index
103            .keys()
104            .map(|ob| data.initial_values.get(ob).copied().unwrap_or_default());
105        let x0 = DVector::from_iterator(n, initial_values);
106
107        let system = system
108            .extend_scalars(|poly| {
109                poly.eval(|id| {
110                    data.interaction_coeffs
111                        .get(id)
112                        .or(data.growth_rates.get(id))
113                        .copied()
114                        .unwrap_or_default()
115                })
116            })
117            .to_numerical();
118        let problem = ODEProblem::new(system, x0).end_time(data.duration);
119        ODEAnalysis::new(problem, ob_index)
120    }
121
122    /// Lotka-Volterra ODE system for an model of a double theory.
123    pub fn lotka_volterra_system(
124        &self,
125        model: &DiscreteDblModel,
126    ) -> (
127        PolynomialSystem<QualifiedName, Parameter<QualifiedName>, u8>,
128        IndexMap<QualifiedName, usize>,
129    ) {
130        let (matrix, ob_index) = self.build_matrix(model);
131        let n = ob_index.len();
132
133        let growth_rate_params = ob_index
134            .keys()
135            .map(|ob| [(1.0, Monomial::generator(ob.clone()))].into_iter().collect());
136        let b = DVector::from_iterator(n, growth_rate_params);
137
138        let system = lotka_volterra_system(&ob_index.keys().cloned().collect_vec(), matrix, b);
139        (system, ob_index)
140    }
141}
142
143#[cfg(test)]
144mod test {
145    use expect_test::expect;
146    use std::rc::Rc;
147
148    use super::*;
149    use crate::stdlib;
150    use crate::{one::Path, zero::name};
151
152    fn builder() -> SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
153        SignedCoefficientBuilder::new(name("Object"))
154            .add_positive(Path::Id(name("Object")))
155            .add_negative(Path::single(name("Negative")))
156    }
157
158    #[test]
159    fn predator_prey_symbolic() {
160        let th = Rc::new(stdlib::theories::th_signed_category());
161        let neg_feedback = stdlib::models::negative_feedback(th);
162        let (sys, _) = builder().lotka_volterra_system(&neg_feedback);
163        let sys = sys.extend_scalars(|coef| coef.map_variables(|name| format!("Param({name})")));
164        let expected = expect!([r#"
165            dx = Param(x) x - Param(negative) x y
166            dy = Param(positive) x y + Param(y) y
167        "#]);
168        expected.assert_eq(&sys.to_string());
169    }
170
171    #[test]
172    fn predator_prey_numerical() {
173        let th = Rc::new(stdlib::theories::th_signed_category());
174        let neg_feedback = stdlib::models::negative_feedback(th);
175
176        let data = LotkaVolterraProblemData {
177            interaction_coeffs: [(name("positive"), 1.0), (name("negative"), 1.0)]
178                .into_iter()
179                .collect(),
180            growth_rates: [(name("x"), 2.0), (name("y"), -1.0)].into_iter().collect(),
181            initial_values: [(name("x"), 1.0), (name("y"), 1.0)].into_iter().collect(),
182            duration: 10.0,
183        };
184
185        let sys = builder().lotka_volterra_analysis(&neg_feedback, data).problem.system;
186        let expected = expect!([r#"
187            dx0 = 2 x0 - x0 x1
188            dx1 = x0 x1 - x1
189        "#]);
190        expected.assert_eq(&sys.to_string());
191    }
192}
catlog/stdlib/analyses/ode/lotka_volterra.rs

catlog/stdlib/analyses/ode/
lotka_volterra.rs
