Skip to content

Commit e7f9913

Browse files
georgefstgeorgefst
committed
Refactor linear ODE systems to use symbolic coefficients
1 parent 6567a16 commit e7f9913

6 files changed

Lines changed: 175 additions & 100 deletions

File tree

‎packages/catlog/src/simulate/ode/linear_ode.rs‎

Lines changed: 0 additions & 76 deletions
This file was deleted.

‎packages/catlog/src/simulate/ode/mod.rs‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -159,12 +159,9 @@ pub(crate) fn textplot_mapped_ode_result<Sys>(
159159

160160
pub mod kuramoto;
161161
#[allow(non_snake_case)]
162-
pub mod linear_ode;
163-
#[allow(non_snake_case)]
164162
pub mod lotka_volterra;
165163
pub mod polynomial;
166164

167165
pub use kuramoto::*;
168-
pub use linear_ode::*;
169166
pub use lotka_volterra::*;
170167
pub use polynomial::*;

‎packages/catlog/src/simulate/ode/polynomial.rs‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ use std::ops::{Add, Neg};
66

77
use derivative::Derivative;
88
use indexmap::IndexMap;
9+
use itertools::Itertools;
910
use nalgebra::DVector;
1011
use num_traits::{One, Pow, Zero};
1112

@@ -115,6 +116,21 @@ where
115116
})
116117
.collect()
117118
}
119+
120+
/// Converts to a single aligned LaTeX environment.
121+
pub fn to_latex_string(&self) -> String
122+
where
123+
Var: Display,
124+
Coef: Display + PartialEq + One + Neg<Output = Coef>,
125+
Exp: Display + PartialEq + One,
126+
{
127+
let equations = self
128+
.to_latex_equations()
129+
.into_iter()
130+
.map(|p| p.lhs + " &= " + &p.rhs)
131+
.join("\\\\\n");
132+
format!("$$\n\\begin{{align*}}\n{}\n\\end{{align*}}\n$$\n", equations)
133+
}
118134
}
119135

120136
#[derive(Debug, PartialEq, Eq)]

‎packages/catlog/src/stdlib/analyses/ode/linear_ode.rs‎

Lines changed: 141 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -4,8 +4,11 @@
44
//! [`linear_ode_analysis`](SignedCoefficientBuilder::linear_ode_analysis).
55
66
use std::collections::HashMap;
7+
use std::hash::Hash;
8+
use std::ops::Add;
79

8-
use nalgebra::DVector;
10+
use itertools::Itertools;
11+
use nalgebra::{DMatrix, DVector};
912

1013
#[cfg(feature = "serde")]
1114
use serde::{Deserialize, Serialize};
@@ -14,8 +17,11 @@ use tsify::Tsify;
1417

1518
use super::{ODEAnalysis, SignedCoefficientBuilder};
1619
use crate::dbl::model::DiscreteDblModel;
17-
use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, linear_polynomial_system};
18-
use crate::{one::QualifiedPath, zero::QualifiedName};
20+
use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem};
21+
use crate::{
22+
one::QualifiedPath,
23+
zero::{QualifiedName, rig::Monomial},
24+
};
1925

2026
/// Data defining a linear ODE problem for a model.
2127
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
@@ -37,6 +43,33 @@ pub struct LinearODEProblemData {
3743
duration: f32,
3844
}
3945

46+
/// Construct a linear (first-order) dynamical system;
47+
/// a semantics for causal loop diagrams.
48+
pub fn linear_polynomial_system<Var, Coef>(
49+
vars: &[Var],
50+
coefficients: DMatrix<Coef>,
51+
) -> PolynomialSystem<Var, Coef, u8>
52+
where
53+
Var: Clone + Hash + Ord,
54+
Coef: Clone + Add<Output = Coef>,
55+
{
56+
PolynomialSystem {
57+
components: coefficients
58+
.row_iter()
59+
.zip(vars)
60+
.map(|(row, i)| {
61+
(
62+
i.clone(),
63+
row.iter()
64+
.zip(vars)
65+
.map(|(a, j)| (a.clone(), Monomial::generator(j.clone())))
66+
.collect(),
67+
)
68+
})
69+
.collect(),
70+
}
71+
}
72+
4073
impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
4174
/// Linear ODE analysis for a model of a double theory.
4275
///
@@ -48,31 +81,40 @@ impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
4881
model: &DiscreteDblModel,
4982
data: LinearODEProblemData,
5083
) -> ODEAnalysis<NumericalPolynomialSystem<u8>> {
51-
let (matrix, ob_index) = self.build_matrix(model, &data.coefficients);
84+
let (matrix, ob_index) = self.build_matrix(model);
5285
let n = ob_index.len();
5386

5487
let initial_values = ob_index
5588
.keys()
5689
.map(|ob| data.initial_values.get(ob).copied().unwrap_or_default());
5790
let x0 = DVector::from_iterator(n, initial_values);
5891

59-
let system = linear_polynomial_system(matrix);
92+
let system = linear_polynomial_system(&ob_index.clone().into_keys().collect_vec(), matrix)
93+
.extend_scalars(|poly| {
94+
poly.eval(|mor| data.coefficients.get(mor).copied().unwrap_or_default())
95+
})
96+
.map(|p| p.normalize())
97+
.to_numerical();
6098
let problem = ODEProblem::new(system, x0).end_time(data.duration);
6199
ODEAnalysis::new(problem, ob_index)
62100
}
63101
}
64102

65103
#[cfg(test)]
104+
#[allow(non_snake_case)]
66105
mod test {
106+
use expect_test::expect;
67107
use std::rc::Rc;
68108

69109
use super::*;
70110
use crate::dbl::model::MutDblModel;
111+
use crate::simulate::ode::textplot_ode_result;
112+
use crate::stdlib;
113+
use crate::stdlib::analyses::ode::Parameter;
71114
use crate::{one::Path, zero::name};
72-
use crate::{simulate::ode::linear_ode, stdlib};
115+
use nalgebra::{dmatrix, dvector};
73116

74-
#[test]
75-
fn neg_loops_pos_connector() {
117+
fn neg_loops_pos_connector_from_theory() -> ODEProblem<NumericalPolynomialSystem<u8>> {
76118
let th = Rc::new(stdlib::theories::th_signed_category());
77119

78120
let mut test_model = DiscreteDblModel::new(th);
@@ -92,10 +134,98 @@ mod test {
92134
.collect(),
93135
duration: 10.0,
94136
};
95-
let analysis = SignedCoefficientBuilder::new(name("Object"))
137+
SignedCoefficientBuilder::new(name("Object"))
96138
.add_positive(Path::Id(name("Object")))
97139
.add_negative(Path::single(name("Negative")))
98-
.linear_ode_analysis(&test_model, data);
99-
assert_eq!(analysis.problem, linear_ode::create_neg_loops_pos_connector());
140+
.linear_ode_analysis(&test_model, data)
141+
.problem
142+
}
143+
144+
fn neg_loops_pos_connector_from_matrix() -> ODEProblem<NumericalPolynomialSystem<u8>> {
145+
ODEProblem::new(matrix_example().to_numerical(), dvector![2.0, 1.0, 1.0]).end_time(10.0)
146+
}
147+
fn matrix_symb_coeff_example() -> PolynomialSystem<QualifiedName, Parameter<QualifiedName>, u8>
148+
{
149+
let A = dmatrix!["aa", "ba", "xa";
150+
"ab", "bb", "xb";
151+
"ax", "bx", "xx"]
152+
.map(|v| [(1.0, Monomial::generator(QualifiedName::from([v])))].into_iter().collect());
153+
linear_polynomial_system(
154+
&vec!["A", "B", "X"].into_iter().map(|v| QualifiedName::from([v])).collect_vec(),
155+
A,
156+
)
157+
}
158+
fn matrix_example() -> PolynomialSystem<QualifiedName, f32, u8> {
159+
let coeffs: HashMap<_, _> = [("aa", -0.3), ("ax", 1.0), ("bx", -2.0), ("xb", 0.5)]
160+
.into_iter()
161+
.map(|(n, v)| (QualifiedName::from([n]), v))
162+
.collect();
163+
matrix_symb_coeff_example()
164+
.extend_scalars(|coeff| coeff.eval(|v| coeffs.get(v).copied().unwrap_or_default()))
165+
.map(|p| p.normalize())
166+
}
167+
168+
#[test]
169+
fn matrix_agrees_with_theory() {
170+
assert_eq!(neg_loops_pos_connector_from_theory(), neg_loops_pos_connector_from_matrix());
171+
}
172+
173+
#[test]
174+
fn linear_solve() {
175+
let problem = neg_loops_pos_connector_from_matrix();
176+
let result = problem.solve_rk4(0.1).unwrap();
177+
let expected = expect![["
178+
⡑⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ 2.0
179+
⠄⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
180+
⠂⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
181+
⡁⠀⠀⣀⠤⠚⠲⣒⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
182+
⠄⡠⠊⠀⠀⠀⠀⠀⠑⠬⣆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠈⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
183+
⠚⢄⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⡤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
184+
⡁⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀⠑⡄⠑⠢⢄⡀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⢘⡔⠊⠉⠉⠒⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀
185+
⠄⠀⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠈⠉⠒⠤⢄⣀⠀⠀⡜⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⢱⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀
186+
⠂⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠉⢱⠓⠢⠤⢄⣀⡀⠀⡠⠃⠀⠀⠀⢇⠀⠀⠀⠀⠀⠈⢢⠀⠀⠀⠀⠀⠀
187+
⡁⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢄⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠈⡝⠉⠑⠒⠒⠢⠼⡤⠤⢄⣀⣀⣀⣀⡱⡀⠀⠀⠀⠀
188+
⡄⢀⠀⡀⢀⠘⡄⢀⠀⡀⢀⠀⡀⢀⠀⡀⢈⢆⡀⢀⠀⡀⢀⡸⡀⢀⠀⡀⢀⢀⡎⢀⠀⡀⢀⠀⡀⢀⢣⡀⢀⠀⡀⢀⠀⡈⢙⡍⡉⢉⠁
189+
⠂⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢢⠀⠀⠀⢠⠃⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠘⢄⠀⠀
190+
⡁⠀⠀⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⡜⠀⠀⠀⢀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠈⢢⠀
191+
⠄⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢱⠣⠤⠤⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁
192+
⠂⠀⠀⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀
193+
⡁⠀⠀⠀⠀⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀
194+
⠄⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀
195+
⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⢀⠇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⡠⠂
196+
⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⢀⡰⠁⠀
197+
⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢄⡀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠁⠀⠀⠀
198+
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ -1.8
199+
0.0 10.0
200+
"]];
201+
expected.assert_eq(&textplot_ode_result(&problem, &result));
202+
}
203+
204+
#[test]
205+
fn latex_symbolic() {
206+
let expected = expect![[r#"
207+
$$
208+
\begin{align*}
209+
\frac{\mathrm{d}}{\mathrm{d}t} A &= aa A + ba B + xa X\\
210+
\frac{\mathrm{d}}{\mathrm{d}t} B &= ab A + bb B + xb X\\
211+
\frac{\mathrm{d}}{\mathrm{d}t} X &= ax A + bx B + xx X
212+
\end{align*}
213+
$$
214+
"#]];
215+
expected.assert_eq(&matrix_symb_coeff_example().to_latex_string());
216+
}
217+
218+
#[test]
219+
fn latex_numerical() {
220+
let expected = expect![[r#"
221+
$$
222+
\begin{align*}
223+
\frac{\mathrm{d}}{\mathrm{d}t} A &= (-0.3) A\\
224+
\frac{\mathrm{d}}{\mathrm{d}t} B &= 0.5 X\\
225+
\frac{\mathrm{d}}{\mathrm{d}t} X &= A + (-2) B
226+
\end{align*}
227+
$$
228+
"#]];
229+
expected.assert_eq(&matrix_example().to_latex_string());
100230
}
101231
}

‎packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,9 @@ impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
5252
model: &DiscreteDblModel,
5353
data: LotkaVolterraProblemData,
5454
) -> ODEAnalysis<NumericalPolynomialSystem<u8>> {
55-
let (matrix, ob_index) = self.build_matrix(model, &data.interaction_coeffs);
55+
let (symbolic_matrix, ob_index) = self.build_matrix(model);
56+
let matrix = symbolic_matrix
57+
.map(|p| p.eval(|id| data.interaction_coeffs.get(id).copied().unwrap_or_default()));
5658
let n = ob_index.len();
5759

5860
let growth_rates =

0 commit comments

Comments
 (0)