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::{BTreeMap, HashMap};
8
9use nalgebra::DVector;
10use num_traits::Zero;
11
12#[cfg(feature = "serde")]
13use serde::{Deserialize, Serialize};
14#[cfg(feature = "serde-wasm")]
15use tsify::Tsify;
16
17use super::ODEAnalysis;
18use crate::dbl::{
19    model::{DiscreteTabModel, FgDblModel, ModalDblModel, MutDblModel, TabEdge},
20    theory::{ModalMorType, ModalObType, TabMorType, TabObType},
21};
22use crate::one::FgCategory;
23use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem};
24use crate::zero::{QualifiedName, alg::Polynomial, name, rig::Monomial};
25
26/// Data defining a mass-action ODE problem for a model.
27#[derive(Clone)]
28#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
29#[cfg_attr(feature = "serde-wasm", derive(Tsify))]
30#[cfg_attr(
31    feature = "serde-wasm",
32    tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object)
33)]
34pub struct MassActionProblemData {
35    /// Map from morphism IDs to rate coefficients (nonnegative reals).
36    rates: HashMap<QualifiedName, f32>,
37
38    /// Map from object IDs to initial values (nonnegative reals).
39    #[cfg_attr(feature = "serde", serde(rename = "initialValues"))]
40    initial_values: HashMap<QualifiedName, f32>,
41
42    /// Duration of simulation.
43    duration: f32,
44}
45
46/// Symbolic parameter in mass-action polynomial system.
47type Parameter<Id> = Polynomial<Id, f32, u8>;
48
49/** Mass-action ODE analysis for Petri nets.
50
51This struct implements the object part of the functorial semantics for reaction
52networks (aka, Petri nets) due to [Baez & Pollard](crate::refs::ReactionNets).
53 */
54pub struct PetriNetMassActionAnalysis {
55    /// Object type for places.
56    pub place_ob_type: ModalObType,
57    /// Morphism type for transitions.
58    pub transition_mor_type: ModalMorType,
59}
60
61impl Default for PetriNetMassActionAnalysis {
62    fn default() -> Self {
63        let ob_type = ModalObType::new(name("Object"));
64        Self {
65            place_ob_type: ob_type.clone(),
66            transition_mor_type: ModalMorType::Zero(ob_type),
67        }
68    }
69}
70
71impl PetriNetMassActionAnalysis {
72    /// Creates a mass-action system with symbolic rate coefficients.
73    pub fn build_system(
74        &self,
75        model: &ModalDblModel,
76    ) -> PolynomialSystem<QualifiedName, Parameter<QualifiedName>, u8> {
77        let mut sys = PolynomialSystem::new();
78        for ob in model.ob_generators_with_type(&self.place_ob_type) {
79            sys.add_term(ob, Polynomial::zero());
80        }
81        for mor in model.mor_generators_with_type(&self.transition_mor_type) {
82            let inputs = model
83                .get_dom(&mor)
84                .and_then(|ob| ob.clone().collect_product(None))
85                .unwrap_or_default();
86            let outputs = model
87                .get_cod(&mor)
88                .and_then(|ob| ob.clone().collect_product(None))
89                .unwrap_or_default();
90
91            let term: Monomial<_, _> =
92                inputs.iter().map(|ob| (ob.clone().unwrap_generator(), 1)).collect();
93            let term: Polynomial<_, _, _> =
94                [(Parameter::generator(mor), term)].into_iter().collect();
95            for input in inputs {
96                sys.add_term(input.unwrap_generator(), -term.clone());
97            }
98            for output in outputs {
99                sys.add_term(output.unwrap_generator(), term.clone());
100            }
101        }
102
103        // Normalize since terms commonly cancel out in mass-action dynamics.
104        sys.normalize()
105    }
106
107    /// Creates a mass-action system with numerical rate coefficients.
108    pub fn build_numerical_system(
109        &self,
110        model: &ModalDblModel,
111        data: MassActionProblemData,
112    ) -> ODEAnalysis<NumericalPolynomialSystem<u8>> {
113        into_numerical_system(self.build_system(model), data)
114    }
115}
116
117/// Mass-action ODE analysis for stock-flow models.
118pub struct StockFlowMassActionAnalysis {
119    /// Object type for stocks.
120    pub stock_ob_type: TabObType,
121    /// Morphism types for flows between stocks.
122    pub flow_mor_type: TabMorType,
123    /// Morphism types for links for stocks to flows.
124    pub link_mor_type: TabMorType,
125}
126
127impl Default for StockFlowMassActionAnalysis {
128    fn default() -> Self {
129        let stock_ob_type = TabObType::Basic(name("Object"));
130        let flow_mor_type = TabMorType::Hom(Box::new(stock_ob_type.clone()));
131        Self {
132            stock_ob_type,
133            flow_mor_type,
134            link_mor_type: TabMorType::Basic(name("Link")),
135        }
136    }
137}
138
139impl StockFlowMassActionAnalysis {
140    /// Creates a mass-action system with symbolic rate coefficients.
141    pub fn build_system(
142        &self,
143        model: &DiscreteTabModel,
144    ) -> PolynomialSystem<QualifiedName, Parameter<QualifiedName>, u8> {
145        let mut terms: HashMap<QualifiedName, Monomial<QualifiedName, u8>> = model
146            .mor_generators_with_type(&self.flow_mor_type)
147            .map(|flow| {
148                let dom = model.mor_generator_dom(&flow).unwrap_basic();
149                (flow, Monomial::generator(dom))
150            })
151            .collect();
152
153        for link in model.mor_generators_with_type(&self.link_mor_type) {
154            let dom = model.mor_generator_dom(&link).unwrap_basic();
155            let path = model.mor_generator_cod(&link).unwrap_tabulated();
156            let Some(TabEdge::Basic(cod)) = path.only() else {
157                panic!("Codomain of link should be basic morphism");
158            };
159            if let Some(term) = terms.get_mut(&cod) {
160                *term = std::mem::take(term) * Monomial::generator(dom);
161            } else {
162                panic!("Codomain of link does not belong to model");
163            };
164        }
165
166        let terms: Vec<_> = terms
167            .into_iter()
168            .map(|(flow, term)| {
169                let param = Parameter::generator(flow.clone());
170                let poly: Polynomial<_, _, _> = [(param, term)].into_iter().collect();
171                (flow, poly)
172            })
173            .collect();
174
175        let mut sys = PolynomialSystem::new();
176        for ob in model.ob_generators_with_type(&self.stock_ob_type) {
177            sys.add_term(ob, Polynomial::zero());
178        }
179        for (flow, term) in terms.iter() {
180            let dom = model.mor_generator_dom(flow).unwrap_basic();
181            sys.add_term(dom, -term.clone());
182        }
183        for (flow, term) in terms {
184            let cod = model.mor_generator_cod(&flow).unwrap_basic();
185            sys.add_term(cod, term);
186        }
187        sys
188    }
189
190    /// Creates a mass-action system with numerical rate coefficients.
191    pub fn build_numerical_system(
192        &self,
193        model: &DiscreteTabModel,
194        data: MassActionProblemData,
195    ) -> ODEAnalysis<NumericalPolynomialSystem<u8>> {
196        into_numerical_system(self.build_system(model), data)
197    }
198}
199
200fn into_numerical_system(
201    sys: PolynomialSystem<QualifiedName, Parameter<QualifiedName>, u8>,
202    data: MassActionProblemData,
203) -> ODEAnalysis<NumericalPolynomialSystem<u8>> {
204    let ob_index: BTreeMap<_, _> =
205        sys.components.keys().cloned().enumerate().map(|(i, x)| (x, i)).collect();
206    let n = ob_index.len();
207
208    let initial_values = ob_index
209        .keys()
210        .map(|ob| data.initial_values.get(ob).copied().unwrap_or_default());
211    let x0 = DVector::from_iterator(n, initial_values);
212
213    let sys = sys
214        .extend_scalars(|poly| poly.eval(|flow| data.rates.get(flow).copied().unwrap_or_default()))
215        .to_numerical();
216
217    let problem = ODEProblem::new(sys, x0).end_time(data.duration);
218    ODEAnalysis::new(problem, ob_index)
219}
220
221#[cfg(test)]
222mod tests {
223    use expect_test::expect;
224    use std::rc::Rc;
225
226    use super::*;
227    use crate::stdlib::{models::*, theories::*};
228
229    #[test]
230    fn backward_link_dynamics() {
231        let th = Rc::new(th_category_links());
232        let model = backward_link(th);
233        let sys = StockFlowMassActionAnalysis::default().build_system(&model);
234        let expected = expect!([r#"
235            dx = ((-1) f) x y
236            dy = f x y
237        "#]);
238        expected.assert_eq(&sys.to_string());
239    }
240
241    #[test]
242    fn catalysis_dynamics() {
243        let th = Rc::new(th_sym_monoidal_category());
244        let model = catalyzed_reaction(th);
245        let sys = PetriNetMassActionAnalysis::default().build_system(&model);
246        let expected = expect!([r#"
247            dc = 0
248            dx = ((-1) f) c x
249            dy = f c x
250        "#]);
251        expected.assert_eq(&sys.to_string());
252    }
253}
catlog/stdlib/analyses/ode/mass_action.rs

catlog/stdlib/analyses/ode/
mass_action.rs
