catlog/stdlib/analyses/ode/
lotka_volterra.rs

1//! Lotka-Volterra ODE analysis of models.
2
3use std::{collections::HashMap, hash::Hash};
4
5use nalgebra::{DMatrix, DVector};
6use ustr::Ustr;
7
8#[cfg(feature = "serde")]
9use serde::{Deserialize, Serialize};
10#[cfg(feature = "serde-wasm")]
11use tsify_next::Tsify;
12
13use super::ODEAnalysis;
14use crate::dbl::model::{DiscreteDblModel, FgDblModel};
15use crate::one::{
16    FgCategory,
17    fin_category::{FinMor, UstrFinCategory},
18};
19use crate::simulate::ode::{LotkaVolterraSystem, ODEProblem};
20
21/// Data defining a Lotka-Volterra ODE problem for a model.
22#[derive(Clone)]
23#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
24#[cfg_attr(feature = "serde-wasm", derive(Tsify))]
25#[cfg_attr(
26    feature = "serde-wasm",
27    tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object)
28)]
29pub struct LotkaVolterraProblemData<Id>
30where
31    Id: Eq + Hash,
32{
33    /// Map from morphism IDs to interaction coefficients (nonnegative reals).
34    #[cfg_attr(feature = "serde", serde(rename = "interactionCoefficients"))]
35    interaction_coeffs: HashMap<Id, f32>,
36
37    /// Map from object IDs to growth rates (arbitrary real numbers).
38    #[cfg_attr(feature = "serde", serde(rename = "growthRates"))]
39    growth_rates: HashMap<Id, f32>,
40
41    /// Map from object IDs to initial values (nonnegative reals).
42    #[cfg_attr(feature = "serde", serde(rename = "initialValues"))]
43    initial_values: HashMap<Id, f32>,
44
45    /// Duration of simulation.
46    duration: f32,
47}
48
49type Model<Id> = DiscreteDblModel<Id, UstrFinCategory>;
50
51/** Lotka-Volterra ODE analysis for models of a double theory.
52
53The main situation we have in mind is the Lotka-Volterra ODE semantics for
54regulatory networks (signed graphs) described in our [*Compositionality*
55paper](crate::refs::RegNets).
56*/
57pub struct LotkaVolterraAnalysis {
58    var_ob_type: Ustr,
59    positive_mor_types: Vec<FinMor<Ustr, Ustr>>,
60    negative_mor_types: Vec<FinMor<Ustr, Ustr>>,
61}
62
63impl LotkaVolterraAnalysis {
64    /// Creates a new Lotka-Volterra analysis for the given object type.
65    pub fn new(var_ob_type: Ustr) -> Self {
66        Self {
67            var_ob_type,
68            positive_mor_types: Vec::new(),
69            negative_mor_types: Vec::new(),
70        }
71    }
72
73    /// Adds a morphism type defining a positive interaction between objects.
74    pub fn add_positive(mut self, mor_type: FinMor<Ustr, Ustr>) -> Self {
75        self.positive_mor_types.push(mor_type);
76        self
77    }
78
79    /// Adds a morphism type defining a negative interaction between objects.
80    pub fn add_negative(mut self, mor_type: FinMor<Ustr, Ustr>) -> Self {
81        self.negative_mor_types.push(mor_type);
82        self
83    }
84
85    /// Creates a Lotka-Volterra system from a model.
86    pub fn create_system<Id: Eq + Clone + Hash + Ord>(
87        &self,
88        model: &Model<Id>,
89        data: LotkaVolterraProblemData<Id>,
90    ) -> ODEAnalysis<Id, LotkaVolterraSystem> {
91        let mut objects: Vec<_> = model.ob_generators_with_type(&self.var_ob_type).collect();
92        objects.sort();
93        let ob_index: HashMap<_, _> =
94            objects.iter().cloned().enumerate().map(|(i, x)| (x, i)).collect();
95
96        let n = objects.len();
97
98        let mut A = DMatrix::from_element(n, n, 0.0f32);
99        for mor_type in self.positive_mor_types.iter() {
100            for mor in model.mor_generators_with_type(mor_type) {
101                let i = *ob_index.get(&model.mor_generator_dom(&mor)).unwrap();
102                let j = *ob_index.get(&model.mor_generator_cod(&mor)).unwrap();
103                A[(j, i)] += data.interaction_coeffs.get(&mor).copied().unwrap_or_default();
104            }
105        }
106        for mor_type in self.negative_mor_types.iter() {
107            for mor in model.mor_generators_with_type(mor_type) {
108                let i = *ob_index.get(&model.mor_generator_dom(&mor)).unwrap();
109                let j = *ob_index.get(&model.mor_generator_cod(&mor)).unwrap();
110                A[(j, i)] -= data.interaction_coeffs.get(&mor).copied().unwrap_or_default();
111            }
112        }
113
114        let growth_rates =
115            objects.iter().map(|ob| data.growth_rates.get(ob).copied().unwrap_or_default());
116        let b = DVector::from_iterator(n, growth_rates);
117
118        let initial_values = objects
119            .iter()
120            .map(|ob| data.initial_values.get(ob).copied().unwrap_or_default());
121        let x0 = DVector::from_iterator(n, initial_values);
122
123        let system = LotkaVolterraSystem::new(A, b);
124        let problem = ODEProblem::new(system, x0).end_time(data.duration);
125        ODEAnalysis::new(problem, ob_index)
126    }
127}
128
129#[cfg(test)]
130mod test {
131    use std::sync::Arc;
132    use ustr::ustr;
133
134    use super::*;
135    use crate::{simulate::ode::lotka_volterra, stdlib};
136
137    #[test]
138    fn predator_prey() {
139        let th = Arc::new(stdlib::theories::th_signed_category());
140        let neg_feedback = stdlib::models::negative_feedback(th);
141
142        let (prey, pred) = (ustr("x"), ustr("y"));
143        let (pos, neg) = (ustr("positive"), ustr("negative"));
144        let data = LotkaVolterraProblemData {
145            interaction_coeffs: [(pos, 1.0), (neg, 1.0)].into_iter().collect(),
146            growth_rates: [(prey, 2.0), (pred, -1.0)].into_iter().collect(),
147            initial_values: [(prey, 1.0), (pred, 1.0)].into_iter().collect(),
148            duration: 10.0,
149        };
150        let analysis = LotkaVolterraAnalysis::new(ustr("Object"))
151            .add_positive(FinMor::Id(ustr("Object")))
152            .add_negative(FinMor::Generator(ustr("Negative")))
153            .create_system(&neg_feedback, data);
154        assert_eq!(analysis.problem, lotka_volterra::create_predator_prey());
155    }
156}