catlog/stdlib/analyses/ode/
polynomial_ode.rs

1//! ODE analysis of models of the logic of systems of polynomial ODEs.
2use std::collections::HashMap;
3
4use indexmap::IndexMap;
5use nalgebra::DVector;
6use num_traits::Zero;
7#[cfg(feature = "serde")]
8use serde::{Deserialize, Serialize};
9#[cfg(feature = "serde-wasm")]
10use tsify::Tsify;
11
12use crate::{
13    dbl::{
14        modal::{ModalMorType, ModalObType, ModeApp},
15        model::{FpDblModel, ModalDblModel, ModalOb, MutDblModel},
16        theory::NonUnital,
17    },
18    simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem},
19    zero::{QualifiedName, alg::Polynomial, name, rig::Monomial},
20};
21
22use super::{ODEAnalysis, Parameter};
23
24/// Data defining an unbalanced mass-action ODE problem for a model.
25#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
26#[cfg_attr(feature = "serde-wasm", derive(Tsify))]
27#[cfg_attr(
28    feature = "serde-wasm",
29    tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object)
30)]
31pub struct PolynomialODEProblemData {
32    /// Map from morphism IDs to coefficients (nonnegative reals).
33    coefficients: HashMap<QualifiedName, f32>,
34
35    /// Map from object IDs to initial values (nonnegative reals).
36    #[cfg_attr(feature = "serde", serde(rename = "initialValues"))]
37    pub initial_values: HashMap<QualifiedName, f32>,
38
39    /// Duration of simulation.
40    pub duration: f32,
41}
42
43/// Polynomial ODE analysis.
44///
45/// The "canonical" analysis for system of polynomial ODEs, namely interpreting
46/// them as actual systems of polynomial ODEs.
47pub struct PolynomialODEAnalysis {
48    /// Object type for variables.
49    pub variable_ob_type: ModalObType,
50    /// Morphism type for positive contributions.
51    pub positive_contribution_mor_type: ModalMorType,
52    /// Morphism type for negative contributions.
53    pub negative_contribution_mor_type: ModalMorType,
54}
55
56impl Default for PolynomialODEAnalysis {
57    fn default() -> Self {
58        Self {
59            variable_ob_type: ModalObType::new(name("State")),
60            positive_contribution_mor_type: ModeApp::new(name("Contribution")).into(),
61            negative_contribution_mor_type: ModeApp::new(name("NegativeContribution")).into(),
62        }
63    }
64}
65
66impl PolynomialODEAnalysis {
67    /// Creates a system with symbolic coefficients.
68    pub fn build_system(
69        &self,
70        model: &ModalDblModel<NonUnital>,
71    ) -> PolynomialSystem<QualifiedName, Parameter<QualifiedName>, i8> {
72        let mut sys = PolynomialSystem::new();
73
74        // Create a variable for each object.
75        for ob in model.ob_generators_with_type(&self.variable_ob_type) {
76            sys.add_term(ob, Polynomial::zero());
77        }
78
79        let make_term = |mor: QualifiedName| {
80            let (Some(ModalOb::List(_, inputs)), Some(output)) =
81                (model.get_dom(&mor), model.get_cod(&mor))
82            else {
83                return None;
84            };
85
86            let term: Monomial<_, _> =
87                inputs.iter().cloned().map(|ob| (ob.unwrap_generator(), 1)).collect();
88            let term: Polynomial<_, _, _> =
89                [(Parameter::generator(mor), term.clone())].into_iter().collect();
90
91            Some((output.clone().unwrap_generator(), term))
92        };
93
94        // Add a monomial with positive sign for each positive contribution.
95        for mor in model.mor_generators_with_type(&self.positive_contribution_mor_type) {
96            if let Some((var, term)) = make_term(mor) {
97                sys.add_term(var, term);
98            }
99        }
100
101        // Add a monomial with negative sign for each negative contribution.
102        for mor in model.mor_generators_with_type(&self.negative_contribution_mor_type) {
103            if let Some((var, term)) = make_term(mor) {
104                sys.add_term(var, -term);
105            }
106        }
107
108        sys.normalize()
109    }
110}
111
112/// Substitutes numerical rate coefficients into a symbolic mass-action system.
113pub fn extend_polynomial_ode_scalars(
114    sys: PolynomialSystem<QualifiedName, Parameter<QualifiedName>, i8>,
115    data: &PolynomialODEProblemData,
116) -> PolynomialSystem<QualifiedName, f32, i8> {
117    let sys = sys.extend_scalars(|poly| {
118        poly.eval(|mor| data.coefficients.get(mor).cloned().unwrap_or_default())
119    });
120
121    sys.normalize()
122}
123
124/// Builds the numerical ODE analysis for a system of polynomial ODEs whose scalars have been substituted.
125pub fn polynomial_ode_analysis(
126    sys: PolynomialSystem<QualifiedName, f32, i8>,
127    data: PolynomialODEProblemData,
128) -> ODEAnalysis<NumericalPolynomialSystem<i8>> {
129    let ob_index: IndexMap<_, _> =
130        sys.components.keys().cloned().enumerate().map(|(i, x)| (x, i)).collect();
131    let n = ob_index.len();
132
133    let initial_values = ob_index
134        .keys()
135        .map(|ob| data.initial_values.get(ob).copied().unwrap_or_default());
136    let x0 = DVector::from_iterator(n, initial_values);
137
138    let num_sys = sys.to_numerical();
139    let problem = ODEProblem::new(num_sys, x0).end_time(data.duration);
140
141    ODEAnalysis::new(problem, ob_index)
142}
143
144#[cfg(test)]
145mod tests {
146    use expect_test::expect;
147    use std::rc::Rc;
148
149    use super::*;
150    use crate::{
151        simulate::ode::LatexEquation,
152        stdlib::{models::*, theories::*},
153        tt,
154    };
155
156    // (Unsigned) Lotka–Volterra dynamics on a two-level model.
157    #[test]
158    fn lotka_volterra_equations() {
159        let th = Rc::new(th_polynomial_ode_system());
160        let model = lotka_volterra_dynamics(th);
161        let sys = PolynomialODEAnalysis::default().build_system(&model);
162        let expected = expect!([r#"
163            dA = A_growth A + BA_interaction A B
164            dB = AB_interaction A B + B_growth B + CB_interaction B C
165            dC = BC_interaction B C + C_growth C
166        "#]);
167        expected.assert_eq(&sys.to_string());
168    }
169
170    // (Unsigned) Lotka–Volterra dynamics on a two-level model with LaTeX.
171    #[test]
172    fn lotka_volterra_equations_latex() {
173        let th = Rc::new(th_polynomial_ode_system());
174        let model = lotka_volterra_dynamics(th);
175        let sys = PolynomialODEAnalysis::default().build_system(&model);
176        let expected = vec![
177            LatexEquation {
178                lhs: "\\frac{\\mathrm{d}}{\\mathrm{d}t} A".to_string(),
179                rhs: "A_growth \\cdot A + BA_interaction \\cdot A \\cdot B".to_string(),
180            },
181            LatexEquation {
182                lhs: "\\frac{\\mathrm{d}}{\\mathrm{d}t} B".to_string(),
183                rhs: "AB_interaction \\cdot A \\cdot B + B_growth \\cdot B + CB_interaction \\cdot B \\cdot C"
184                    .to_string(),
185            },
186            LatexEquation {
187                lhs: "\\frac{\\mathrm{d}}{\\mathrm{d}t} C".to_string(),
188                rhs: "BC_interaction \\cdot B \\cdot C + C_growth \\cdot C".to_string(),
189            },
190        ];
191        assert_eq!(expected, sys.to_latex_equations());
192    }
193
194    // DoubleTT elaboration from text.
195    #[test]
196    fn ode_system_from_text() {
197        let th = Rc::new(th_polynomial_ode_system());
198        let model = tt::modelgen::Model::from_text(
199            &th.into(),
200            "[
201                X : State,
202                Y : State,
203                A : State,
204                f : Contribution[[X, Y, Y], A],
205                g : Contribution[[X, X], Y],
206                h : Contribution[[A], X],
207            ]",
208        );
209        let model = model.unwrap().as_modal_non_unital().unwrap();
210        let sys = PolynomialODEAnalysis::default().build_system(&model);
211        let expected = expect!([r#"
212            dX = h A
213            dY = g X^2
214            dA = f X Y^2
215        "#]);
216        expected.assert_eq(&sys.to_string());
217    }
218}
catlog/stdlib/analyses/ode/polynomial_ode.rs

catlog/stdlib/analyses/ode/
polynomial_ode.rs
