catlog/stdlib/analyses/ode/
mass_action.rs

1/*! Mass-action ODE analysis of models.
2
3Such ODEs are based on the *law of mass action* familiar from chemistry and
4mathematical epidemiology.
5 */
6
7use std::collections::HashMap;
8use std::hash::{BuildHasherDefault, Hash};
9
10use nalgebra::DVector;
11use num_traits::Zero;
12use ustr::{IdentityHasher, Ustr, ustr};
13
14#[cfg(feature = "serde")]
15use serde::{Deserialize, Serialize};
16#[cfg(feature = "serde-wasm")]
17use tsify_next::Tsify;
18
19use super::ODEAnalysis;
20use crate::dbl::{
21    model::{DiscreteTabModel, FgDblModel, TabEdge},
22    theory::{TabMorType, TabObType},
23};
24use crate::one::FgCategory;
25use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem};
26use crate::zero::{alg::Polynomial, rig::Monomial};
27
28/// Data defining a mass-action ODE problem for a model.
29#[derive(Clone)]
30#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
31#[cfg_attr(feature = "serde-wasm", derive(Tsify))]
32#[cfg_attr(
33    feature = "serde-wasm",
34    tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object)
35)]
36pub struct MassActionProblemData<Id>
37where
38    Id: Eq + Hash,
39{
40    /// Map from morphism IDs to rate coefficients (nonnegative reals).
41    rates: HashMap<Id, f32>,
42
43    /// Map from object IDs to initial values (nonnegative reals).
44    #[cfg_attr(feature = "serde", serde(rename = "initialValues"))]
45    initial_values: HashMap<Id, f32>,
46
47    /// Duration of simulation.
48    duration: f32,
49}
50
51type Parameter<Id> = Polynomial<Id, f32, u8>;
52type StockFlowModel<Id> = DiscreteTabModel<Id, Ustr, BuildHasherDefault<IdentityHasher>>;
53
54/** Mass-action ODE analysis for stock-flow models.
55
56Mass action dynamics TODO
57 */
58pub struct StockFlowMassActionAnalysis {
59    /// Object type for stocks.
60    pub stock_ob_type: TabObType<Ustr, Ustr>,
61    /// Morphism types for flows between stocks.
62    pub flow_mor_type: TabMorType<Ustr, Ustr>,
63    /// Morphism types for links for stocks to flows.
64    pub link_mor_type: TabMorType<Ustr, Ustr>,
65}
66
67impl Default for StockFlowMassActionAnalysis {
68    fn default() -> Self {
69        let stock_ob_type = TabObType::Basic(ustr("Object"));
70        let flow_mor_type = TabMorType::Hom(Box::new(stock_ob_type.clone()));
71        Self {
72            stock_ob_type,
73            flow_mor_type,
74            link_mor_type: TabMorType::Basic(ustr("Link")),
75        }
76    }
77}
78
79impl StockFlowMassActionAnalysis {
80    /** Creates a mass-action system from a model.
81
82    The resulting system has symbolic rate coefficients.
83     */
84    pub fn create_system<Id: Eq + Clone + Hash + Ord>(
85        &self,
86        model: &StockFlowModel<Id>,
87    ) -> PolynomialSystem<Id, Parameter<Id>, u8> {
88        let mut terms: HashMap<Id, Monomial<Id, u8>> = model
89            .mor_generators_with_type(&self.flow_mor_type)
90            .map(|flow| {
91                let dom = model.mor_generator_dom(&flow).unwrap_basic();
92                (flow, Monomial::generator(dom))
93            })
94            .collect();
95
96        for link in model.mor_generators_with_type(&self.link_mor_type) {
97            let dom = model.mor_generator_dom(&link).unwrap_basic();
98            let path = model.mor_generator_cod(&link).unwrap_tabulated();
99            let Some(TabEdge::Basic(cod)) = path.only() else {
100                panic!("Codomain of link should be basic morphism");
101            };
102            if let Some(term) = terms.get_mut(&cod) {
103                *term = std::mem::take(term) * Monomial::generator(dom);
104            } else {
105                panic!("Codomain of link does not belong to model");
106            };
107        }
108
109        let terms: Vec<(Id, Polynomial<Id, Parameter<Id>, u8>)> = terms
110            .into_iter()
111            .map(|(flow, term)| {
112                let param = Parameter::generator(flow.clone());
113                (flow, [(param, term)].into_iter().collect())
114            })
115            .collect();
116
117        let mut sys: PolynomialSystem<Id, Parameter<Id>, u8> = PolynomialSystem::new();
118        for ob in model.ob_generators_with_type(&self.stock_ob_type) {
119            sys.add_term(ob, Polynomial::zero());
120        }
121        for (flow, term) in terms.iter() {
122            let dom = model.mor_generator_dom(flow).unwrap_basic();
123            sys.add_term(dom, -term.clone());
124        }
125        for (flow, term) in terms {
126            let cod = model.mor_generator_cod(&flow).unwrap_basic();
127            sys.add_term(cod, term);
128        }
129        sys
130    }
131
132    /** Creates a numerical mass-action system from a model.
133
134    The resulting system has numerical rate coefficients and is ready to solve.
135     */
136    pub fn create_numerical_system<Id: Eq + Clone + Hash + Ord>(
137        &self,
138        model: &StockFlowModel<Id>,
139        data: MassActionProblemData<Id>,
140    ) -> ODEAnalysis<Id, NumericalPolynomialSystem<u8>> {
141        let sys = self.create_system(model);
142
143        let objects: Vec<_> = sys.components.keys().cloned().collect();
144        let initial_values = objects
145            .iter()
146            .map(|ob| data.initial_values.get(ob).copied().unwrap_or_default());
147        let x0 = DVector::from_iterator(objects.len(), initial_values);
148
149        let sys = sys
150            .extend_scalars(|poly| {
151                poly.eval(|flow| data.rates.get(flow).copied().unwrap_or_default())
152            })
153            .to_numerical();
154
155        let problem = ODEProblem::new(sys, x0).end_time(data.duration);
156        let ob_index: HashMap<_, _> =
157            objects.into_iter().enumerate().map(|(i, x)| (x, i)).collect();
158        ODEAnalysis::new(problem, ob_index)
159    }
160}
161
162#[cfg(test)]
163mod tests {
164    use expect_test::expect;
165    use std::sync::Arc;
166
167    use super::*;
168    use crate::stdlib::{models::backward_link, theories::th_category_links};
169
170    #[test]
171    fn backward_link_dynamics() {
172        let th = Arc::new(th_category_links());
173        let model = backward_link(th);
174        let analysis: StockFlowMassActionAnalysis = Default::default();
175        let sys = analysis.create_system(&model);
176        let expected = expect!([r#"
177            dx = ((-1) f) x y
178            dy = f x y
179        "#]);
180        expected.assert_eq(&sys.to_string());
181    }
182}
catlog/stdlib/analyses/ode/mass_action.rs

catlog/stdlib/analyses/ode/
mass_action.rs
