catlog/stdlib/analyses/ode/
linear_ode.rs

1//! Constant-coefficient linear first-order ODE analysis of models.
2//!
3//! The main entry point for this module is
4//! [`linear_ode_analysis`](SignedCoefficientBuilder::linear_ode_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};
13use num_traits::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, rig::Monomial},
26};
27
28/// Data defining a linear 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 LinearODEProblemData {
36    /// Map from morphism IDs to interaction coefficients (nonnegative reals).
37    #[cfg_attr(feature = "serde", serde(rename = "coefficients"))]
38    coefficients: HashMap<QualifiedName, f32>,
39
40    /// Map from object IDs to initial values (nonnegative reals).
41    #[cfg_attr(feature = "serde", serde(rename = "initialValues"))]
42    initial_values: HashMap<QualifiedName, f32>,
43
44    /// Duration of simulation.
45    duration: f32,
46}
47
48/// Construct a linear (first-order) dynamical system;
49/// a semantics for causal loop diagrams.
50pub fn linear_polynomial_system<Var, Coef>(
51    vars: &[Var],
52    coefficients: DMatrix<Coef>,
53) -> PolynomialSystem<Var, Coef, u8>
54where
55    Var: Clone + Hash + Ord,
56    Coef: Clone + Add<Output = Coef> + Zero,
57{
58    let system = PolynomialSystem {
59        components: coefficients
60            .row_iter()
61            .zip(vars)
62            .map(|(row, i)| {
63                (
64                    i.clone(),
65                    row.iter()
66                        .zip(vars)
67                        .map(|(a, j)| (a.clone(), Monomial::generator(j.clone())))
68                        .collect(),
69                )
70            })
71            .collect(),
72    };
73    system.normalize()
74}
75
76impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
77    /// Linear ODE analysis for a model of a double theory.
78    ///
79    /// This analysis is a special case of linear ODE analysis for *extended* causal
80    /// loop diagrams but can serve as a simple/naive semantics for causal loop
81    /// diagrams, hopefully useful for toy models and demonstration purposes.
82    pub fn linear_ode_analysis(
83        &self,
84        model: &DiscreteDblModel,
85        data: LinearODEProblemData,
86    ) -> ODEAnalysis<NumericalPolynomialSystem<u8>> {
87        let (system, ob_index) = self.linear_ode_system(model);
88        let n = ob_index.len();
89
90        let initial_values = ob_index
91            .keys()
92            .map(|ob| data.initial_values.get(ob).copied().unwrap_or_default());
93        let x0 = DVector::from_iterator(n, initial_values);
94
95        let system = system
96            .extend_scalars(|poly| {
97                poly.eval(|id| data.coefficients.get(id).copied().unwrap_or_default())
98            })
99            .to_numerical();
100        let problem = ODEProblem::new(system, x0).end_time(data.duration);
101        ODEAnalysis::new(problem, ob_index)
102    }
103
104    /// Linear ODE system for a model of a double theory.
105    pub fn linear_ode_system(
106        &self,
107        model: &DiscreteDblModel,
108    ) -> (
109        PolynomialSystem<QualifiedName, Parameter<QualifiedName>, u8>,
110        IndexMap<QualifiedName, usize>,
111    ) {
112        let (matrix, ob_index) = self.build_matrix(model);
113        let system = linear_polynomial_system(&ob_index.keys().cloned().collect_vec(), matrix);
114        (system, ob_index)
115    }
116}
117
118#[cfg(test)]
119mod test {
120    use expect_test::expect;
121    use std::rc::Rc;
122
123    use super::*;
124    use crate::stdlib;
125    use crate::{one::Path, zero::name};
126
127    fn builder() -> SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
128        SignedCoefficientBuilder::new(name("Object"))
129            .add_positive(Path::Id(name("Object")))
130            .add_negative(Path::single(name("Negative")))
131    }
132
133    #[test]
134    fn negative_feedback_symbolic() {
135        let th = Rc::new(stdlib::theories::th_signed_category());
136        let neg_feedback = stdlib::models::negative_feedback(th);
137        let (sys, _) = builder().linear_ode_system(&neg_feedback);
138        let expected = expect![[r#"
139            dx = -negative y
140            dy = positive x
141        "#]];
142        expected.assert_eq(&sys.to_string());
143    }
144
145    #[test]
146    fn negative_feedback_numerical() {
147        let th = Rc::new(stdlib::theories::th_signed_category());
148        let neg_feedback = stdlib::models::negative_feedback(th);
149
150        let data = LinearODEProblemData {
151            coefficients: [(name("positive"), 2.0), (name("negative"), 1.0)].into_iter().collect(),
152            initial_values: [(name("x"), 1.0), (name("y"), 1.0)].into_iter().collect(),
153            duration: 10.0,
154        };
155
156        let sys = builder().linear_ode_analysis(&neg_feedback, data).problem.system;
157        let expected = expect![[r#"
158            dx0 = -x1
159            dx1 = 2 x0
160        "#]];
161        expected.assert_eq(&sys.to_string());
162    }
163}
catlog/stdlib/analyses/ode/linear_ode.rs

catlog/stdlib/analyses/ode/
linear_ode.rs
