diff --git a/packages/catlog/src/simulate/ode/linear_ode.rs b/packages/catlog/src/simulate/ode/linear_ode.rs deleted file mode 100644 index 01c74785c..000000000 --- a/packages/catlog/src/simulate/ode/linear_ode.rs +++ /dev/null @@ -1,81 +0,0 @@ -//! Constant-coefficient linear first-order differential equations. - -use nalgebra::{DMatrix, DVector}; - -#[cfg(test)] -use super::ODEProblem; -use super::ODESystem; - -/// A (constant-coefficient) linear (first-order) dynamical system. -/// -/// A system of linear first-order ODEs with constant coefficients; a semantics for -/// causal loop diagrams. -#[derive(Clone, Debug, PartialEq)] -pub struct LinearODESystem { - coefficients: DMatrix, -} - -impl LinearODESystem { - /// Constructs a linear ODE system with the given coefficient matrix. - pub fn new(A: DMatrix) -> Self { - Self { coefficients: A } - } -} - -impl ODESystem for LinearODESystem { - fn vector_field(&self, dx: &mut DVector, x: &DVector, _t: f32) { - let A = &self.coefficients; - *dx = A * x - } -} - -#[cfg(test)] -pub(crate) fn create_neg_loops_pos_connector() -> ODEProblem { - use nalgebra::{dmatrix, dvector}; - - let A = dmatrix![-0.3, 0.0, 0.0; - 0.0, 0.0, 0.5; - 1.0, -2.0, 0.0]; - let system = LinearODESystem::new(A); - let initial = dvector![2.0, 1.0, 1.0]; - ODEProblem::new(system, initial).end_time(10.0) -} - -#[cfg(test)] -mod tests { - use expect_test::expect; - - use super::super::textplot_ode_result; - use super::*; - - #[test] - fn neg_loops_pos_connector() { - let problem = create_neg_loops_pos_connector(); - let result = problem.solve_rk4(0.1).unwrap(); - let expected = expect![[" - ⡑⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ 2.0 - ⠄⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠂⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡁⠀⠀⣀⠤⠚⠲⣒⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠄⡠⠊⠀⠀⠀⠀⠀⠑⠬⣆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠈⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠚⢄⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⡤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡁⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀⠑⡄⠑⠢⢄⡀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⢘⡔⠊⠉⠉⠒⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠄⠀⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠈⠉⠒⠤⢄⣀⠀⠀⡜⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⢱⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀ - ⠂⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠉⢱⠓⠢⠤⢄⣀⡀⠀⡠⠃⠀⠀⠀⢇⠀⠀⠀⠀⠀⠈⢢⠀⠀⠀⠀⠀⠀ - ⡁⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢄⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠈⡝⠉⠑⠒⠒⠢⠼⡤⠤⢄⣀⣀⣀⣀⡱⡀⠀⠀⠀⠀ - ⡄⢀⠀⡀⢀⠘⡄⢀⠀⡀⢀⠀⡀⢀⠀⡀⢈⢆⡀⢀⠀⡀⢀⡸⡀⢀⠀⡀⢀⢀⡎⢀⠀⡀⢀⠀⡀⢀⢣⡀⢀⠀⡀⢀⠀⡈⢙⡍⡉⢉⠁ - ⠂⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢢⠀⠀⠀⢠⠃⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠘⢄⠀⠀ - ⡁⠀⠀⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⡜⠀⠀⠀⢀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠈⢢⠀ - ⠄⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢱⠣⠤⠤⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁ - ⠂⠀⠀⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡁⠀⠀⠀⠀⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀ - ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⢀⠇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⡠⠂ - ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⢀⡰⠁⠀ - ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢄⡀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠁⠀⠀⠀ - ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ -1.8 - 0.0 10.0 - "]]; - expected.assert_eq(&textplot_ode_result(&problem, &result)); - } -} diff --git a/packages/catlog/src/simulate/ode/lotka_volterra.rs b/packages/catlog/src/simulate/ode/lotka_volterra.rs deleted file mode 100644 index 82561725a..000000000 --- a/packages/catlog/src/simulate/ode/lotka_volterra.rs +++ /dev/null @@ -1,82 +0,0 @@ -//! Lotka-Volterra differential equations. - -use nalgebra::{DMatrix, DVector}; - -#[cfg(test)] -use super::ODEProblem; -use super::ODESystem; - -/// A Lotka-Volterra dynamical system. -/// -/// A system of ODEs that is affine in its *logarithmic* derivative. These are -/// sometimes called the "generalized Lotka-Volterra equations." For more, see -/// [Wikipedia](https://en.wikipedia.org/wiki/Generalized_Lotka%E2%80%93Volterra_equation). -#[derive(Clone, Debug, PartialEq)] -pub struct LotkaVolterraSystem { - interaction_coeffs: DMatrix, - growth_rates: DVector, -} - -impl LotkaVolterraSystem { - /// Constructs a new Lokta-Volterra system with the given parameters. - pub fn new(A: DMatrix, b: DVector) -> Self { - Self { interaction_coeffs: A, growth_rates: b } - } -} - -impl ODESystem for LotkaVolterraSystem { - fn vector_field(&self, dx: &mut DVector, x: &DVector, _t: f32) { - let A = &self.interaction_coeffs; - let b = &self.growth_rates; - *dx = (A * x + b).component_mul(x); - } -} - -#[cfg(test)] -pub(crate) fn create_predator_prey() -> ODEProblem { - let A = DMatrix::from_row_slice(2, 2, &[0.0, -1.0, 1.0, 0.0]); - let b = DVector::from_column_slice(&[2.0, -1.0]); - let sys = LotkaVolterraSystem::new(A, b); - - let x0 = DVector::from_column_slice(&[1.0, 1.0]); - ODEProblem::new(sys, x0).end_time(10.0) -} - -#[cfg(test)] -mod tests { - use expect_test::expect; - - use super::super::textplot_ode_result; - use super::*; - - #[test] - fn predator_prey() { - let problem = create_predator_prey(); - let result = problem.solve_rk4(0.1).unwrap(); - let expected = expect![[" - ⡁⠀⠀⠀⠀⠀⠀⠀⢠⠊⢢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ 3.5 - ⠄⠀⠀⠀⠀⠀⠀⠀⡇⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡜⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠂⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠇⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡁⠀⠀⠀⠀⠀⠀⡎⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠄⠀⠀⠀⠀⠀⢀⠇⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠂⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠇⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡁⠀⠀⠀⠀⠀⡎⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠄⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡎⠀⠀⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠂⠀⠀⠀⠀⣸⡀⠀⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡇⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡁⠀⠀⠀⡎⡜⢣⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡰⢹⠸⡀⠀⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠄⠀⠀⡸⠀⡇⠈⡆⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⢰⠁⡜⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠄ - ⠂⠀⢠⠃⢸⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⡎⢀⠇⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⡸⠀ - ⡁⠀⡎⠀⡎⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⡸⠀⡸⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⢠⠃⠀ - ⠄⢰⠁⢰⠁⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⢀⠇⢀⠇⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⠀⡎⠀⠀ - ⢂⠇⢀⠇⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⡜⠀⡜⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢄⠀⠀⠀⢰⠁⡰⠁ - ⡝⡠⠊⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⣀⢰⣁⠜⠀⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⣀⢀⢇⡰⠁⠀ - ⠍⠀⠀⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠋⠀⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡏⠁⠀⠀⠀ - ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡜⠀⠀⠀⠀⠀ - ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢆⠀⠀⠀⠀⠀⠀⠀⠀⡠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀ - ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⡀⠀⠀⢀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⡀⠀⠀⢀⡠⠔⠁⠀⠀⠀⠀⠀⠀⠀ - ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ 0.4 - 0.0 10.0 - "]]; - expected.assert_eq(&textplot_ode_result(&problem, &result)); - } -} diff --git a/packages/catlog/src/simulate/ode/mod.rs b/packages/catlog/src/simulate/ode/mod.rs index 9ff29c5de..6e6e1e5a6 100644 --- a/packages/catlog/src/simulate/ode/mod.rs +++ b/packages/catlog/src/simulate/ode/mod.rs @@ -158,13 +158,7 @@ pub(crate) fn textplot_mapped_ode_result( } pub mod kuramoto; -#[allow(non_snake_case)] -pub mod linear_ode; -#[allow(non_snake_case)] -pub mod lotka_volterra; pub mod polynomial; pub use kuramoto::*; -pub use linear_ode::*; -pub use lotka_volterra::*; pub use polynomial::*; diff --git a/packages/catlog/src/simulate/ode/polynomial.rs b/packages/catlog/src/simulate/ode/polynomial.rs index 98c5a71dc..a448ffae7 100644 --- a/packages/catlog/src/simulate/ode/polynomial.rs +++ b/packages/catlog/src/simulate/ode/polynomial.rs @@ -6,6 +6,7 @@ use std::ops::{Add, Neg}; use derivative::Derivative; use indexmap::IndexMap; +use itertools::Itertools; use nalgebra::DVector; use num_traits::{One, Pow, Zero}; @@ -115,6 +116,21 @@ where }) .collect() } + + /// Converts to a single aligned LaTeX environment. + pub fn to_latex_string(&self) -> String + where + Var: Display, + Coef: Display + PartialEq + One + Neg, + Exp: Display + PartialEq + One, + { + let equations = self + .to_latex_equations() + .into_iter() + .map(|p| p.lhs + " &= " + &p.rhs) + .join("\\\\\n"); + format!("$$\n\\begin{{align*}}\n{}\n\\end{{align*}}\n$$\n", equations) + } } #[derive(Debug, PartialEq, Eq)] @@ -184,6 +200,7 @@ where /// /// Such a system is ready for use in numerical solvers: the coefficients are /// floating point numbers and the variables are consecutive integer indices. +#[derive(PartialEq, Debug)] pub struct NumericalPolynomialSystem { /// Components of the vector field. pub components: Vec>, diff --git a/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs b/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs index 40c06ed17..4b4d8f924 100644 --- a/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs +++ b/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs @@ -4,8 +4,11 @@ //! [`linear_ode_analysis`](SignedCoefficientBuilder::linear_ode_analysis). use std::collections::HashMap; +use std::hash::Hash; +use std::ops::Add; -use nalgebra::DVector; +use itertools::Itertools; +use nalgebra::{DMatrix, DVector}; #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; @@ -14,8 +17,11 @@ use tsify::Tsify; use super::{ODEAnalysis, SignedCoefficientBuilder}; use crate::dbl::model::DiscreteDblModel; -use crate::simulate::ode::{LinearODESystem, ODEProblem}; -use crate::{one::QualifiedPath, zero::QualifiedName}; +use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem}; +use crate::{ + one::QualifiedPath, + zero::{QualifiedName, rig::Monomial}, +}; /// Data defining a linear ODE problem for a model. #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] @@ -37,6 +43,33 @@ pub struct LinearODEProblemData { duration: f32, } +/// Construct a linear (first-order) dynamical system; +/// a semantics for causal loop diagrams. +pub fn linear_polynomial_system( + vars: &[Var], + coefficients: DMatrix, +) -> PolynomialSystem +where + Var: Clone + Hash + Ord, + Coef: Clone + Add, +{ + PolynomialSystem { + components: coefficients + .row_iter() + .zip(vars) + .map(|(row, i)| { + ( + i.clone(), + row.iter() + .zip(vars) + .map(|(a, j)| (a.clone(), Monomial::generator(j.clone()))) + .collect(), + ) + }) + .collect(), + } +} + impl SignedCoefficientBuilder { /// Linear ODE analysis for a model of a double theory. /// @@ -47,8 +80,8 @@ impl SignedCoefficientBuilder { &self, model: &DiscreteDblModel, data: LinearODEProblemData, - ) -> ODEAnalysis { - let (matrix, ob_index) = self.build_matrix(model, &data.coefficients); + ) -> ODEAnalysis> { + let (matrix, ob_index) = self.build_matrix(model); let n = ob_index.len(); let initial_values = ob_index @@ -56,23 +89,32 @@ impl SignedCoefficientBuilder { .map(|ob| data.initial_values.get(ob).copied().unwrap_or_default()); let x0 = DVector::from_iterator(n, initial_values); - let system = LinearODESystem::new(matrix); + let system = linear_polynomial_system(&ob_index.clone().into_keys().collect_vec(), matrix) + .extend_scalars(|poly| { + poly.eval(|mor| data.coefficients.get(mor).copied().unwrap_or_default()) + }) + .map(|p| p.normalize()) + .to_numerical(); let problem = ODEProblem::new(system, x0).end_time(data.duration); ODEAnalysis::new(problem, ob_index) } } #[cfg(test)] +#[allow(non_snake_case)] mod test { + use expect_test::expect; use std::rc::Rc; use super::*; use crate::dbl::model::MutDblModel; + use crate::simulate::ode::textplot_ode_result; + use crate::stdlib; + use crate::stdlib::analyses::ode::Parameter; use crate::{one::Path, zero::name}; - use crate::{simulate::ode::linear_ode, stdlib}; + use nalgebra::{dmatrix, dvector}; - #[test] - fn neg_loops_pos_connector() { + fn neg_loops_pos_connector_from_theory() -> ODEProblem> { let th = Rc::new(stdlib::theories::th_signed_category()); let mut test_model = DiscreteDblModel::new(th); @@ -92,10 +134,98 @@ mod test { .collect(), duration: 10.0, }; - let analysis = SignedCoefficientBuilder::new(name("Object")) + SignedCoefficientBuilder::new(name("Object")) .add_positive(Path::Id(name("Object"))) .add_negative(Path::single(name("Negative"))) - .linear_ode_analysis(&test_model, data); - assert_eq!(analysis.problem, linear_ode::create_neg_loops_pos_connector()); + .linear_ode_analysis(&test_model, data) + .problem + } + + fn neg_loops_pos_connector_from_matrix() -> ODEProblem> { + ODEProblem::new(matrix_example().to_numerical(), dvector![2.0, 1.0, 1.0]).end_time(10.0) + } + fn matrix_symb_coeff_example() -> PolynomialSystem, u8> + { + let A = dmatrix!["aa", "ba", "xa"; + "ab", "bb", "xb"; + "ax", "bx", "xx"] + .map(|v| [(1.0, Monomial::generator(QualifiedName::from([v])))].into_iter().collect()); + linear_polynomial_system( + &vec!["A", "B", "X"].into_iter().map(|v| QualifiedName::from([v])).collect_vec(), + A, + ) + } + fn matrix_example() -> PolynomialSystem { + let coeffs: HashMap<_, _> = [("aa", -0.3), ("ax", 1.0), ("bx", -2.0), ("xb", 0.5)] + .into_iter() + .map(|(n, v)| (QualifiedName::from([n]), v)) + .collect(); + matrix_symb_coeff_example() + .extend_scalars(|coeff| coeff.eval(|v| coeffs.get(v).copied().unwrap_or_default())) + .map(|p| p.normalize()) + } + + #[test] + fn matrix_agrees_with_theory() { + assert_eq!(neg_loops_pos_connector_from_theory(), neg_loops_pos_connector_from_matrix()); + } + + #[test] + fn linear_solve() { + let problem = neg_loops_pos_connector_from_matrix(); + let result = problem.solve_rk4(0.1).unwrap(); + let expected = expect![[" + ⡑⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ 2.0 + ⠄⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠀⣀⠤⠚⠲⣒⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⡠⠊⠀⠀⠀⠀⠀⠑⠬⣆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠈⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠚⢄⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⡤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀⠑⡄⠑⠢⢄⡀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⢘⡔⠊⠉⠉⠒⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠈⠉⠒⠤⢄⣀⠀⠀⡜⠀⠀⠀⠀⠀⠀⠀⢀⠔⠁⢱⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠉⢱⠓⠢⠤⢄⣀⡀⠀⡠⠃⠀⠀⠀⢇⠀⠀⠀⠀⠀⠈⢢⠀⠀⠀⠀⠀⠀ + ⡁⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢄⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠈⡝⠉⠑⠒⠒⠢⠼⡤⠤⢄⣀⣀⣀⣀⡱⡀⠀⠀⠀⠀ + ⡄⢀⠀⡀⢀⠘⡄⢀⠀⡀⢀⠀⡀⢀⠀⡀⢈⢆⡀⢀⠀⡀⢀⡸⡀⢀⠀⡀⢀⢀⡎⢀⠀⡀⢀⠀⡀⢀⢣⡀⢀⠀⡀⢀⠀⡈⢙⡍⡉⢉⠁ + ⠂⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢢⠀⠀⠀⢠⠃⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠘⢄⠀⠀ + ⡁⠀⠀⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⡜⠀⠀⠀⢀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠈⢢⠀ + ⠄⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢱⠣⠤⠤⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁ + ⠂⠀⠀⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠀⠀⠀⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⢀⠇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⡠⠂ + ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⢀⡰⠁⠀ + ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢄⡀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠁⠀⠀⠀ + ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ -1.8 + 0.0 10.0 + "]]; + expected.assert_eq(&textplot_ode_result(&problem, &result)); + } + + #[test] + fn latex_symbolic() { + let expected = expect![[r#" + $$ + \begin{align*} + \frac{\mathrm{d}}{\mathrm{d}t} A &= aa A + ba B + xa X\\ + \frac{\mathrm{d}}{\mathrm{d}t} B &= ab A + bb B + xb X\\ + \frac{\mathrm{d}}{\mathrm{d}t} X &= ax A + bx B + xx X + \end{align*} + $$ + "#]]; + expected.assert_eq(&matrix_symb_coeff_example().to_latex_string()); + } + + #[test] + fn latex_numerical() { + let expected = expect![[r#" + $$ + \begin{align*} + \frac{\mathrm{d}}{\mathrm{d}t} A &= (-0.3) A\\ + \frac{\mathrm{d}}{\mathrm{d}t} B &= 0.5 X\\ + \frac{\mathrm{d}}{\mathrm{d}t} X &= A + (-2) B + \end{align*} + $$ + "#]]; + expected.assert_eq(&matrix_example().to_latex_string()); } } diff --git a/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs b/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs index 323e32701..6996ef4f1 100644 --- a/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs +++ b/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs @@ -4,9 +4,13 @@ //! [`lotka_volterra_analysis`](SignedCoefficientBuilder::lotka_volterra_analysis). use std::collections::HashMap; +use std::hash::Hash; +use std::ops::Add; -use nalgebra::DVector; +use itertools::Itertools; +use nalgebra::{DMatrix, DVector, Scalar}; +use num_traits::One; #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; #[cfg(feature = "serde-wasm")] @@ -14,8 +18,12 @@ use tsify::Tsify; use super::{ODEAnalysis, SignedCoefficientBuilder}; use crate::dbl::model::DiscreteDblModel; -use crate::simulate::ode::{LotkaVolterraSystem, ODEProblem}; -use crate::{one::QualifiedPath, zero::QualifiedName}; +use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem}; +use crate::zero::alg::Polynomial; +use crate::{ + one::QualifiedPath, + zero::{QualifiedName, rig::Monomial}, +}; /// Data defining a Lotka-Volterra ODE problem for a model. #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] @@ -41,6 +49,41 @@ pub struct LotkaVolterraProblemData { duration: f32, } +/// Construct a Lotka-Volterra dynamical system. +/// +/// A system of ODEs that is affine in its *logarithmic* derivative. These are +/// sometimes called the "generalized Lotka-Volterra equations." For more, see +/// [Wikipedia](https://en.wikipedia.org/wiki/Generalized_Lotka%E2%80%93Volterra_equation). +pub fn lotka_volterra_system( + vars: &[Var], + interaction_coeffs: DMatrix, + growth_rates: DVector, +) -> PolynomialSystem +where + Var: Clone + Hash + Ord, + Coef: Clone + Add + One + Scalar, +{ + PolynomialSystem { + components: interaction_coeffs + .row_iter() + .zip(vars) + .zip(&growth_rates) + .map(|((row, i), r)| { + ( + i.clone(), + Polynomial::<_, Coef, _>::generator(i.clone()) + * (row + .iter() + .zip(vars) + .map(|(a, j)| (a.clone(), Monomial::generator(j.clone()))) + .collect::>() + + r.clone()), + ) + }) + .collect(), + } +} + impl SignedCoefficientBuilder { /// Lotka-Volterra ODE analysis for a model of a double theory. /// @@ -51,20 +94,32 @@ impl SignedCoefficientBuilder { &self, model: &DiscreteDblModel, data: LotkaVolterraProblemData, - ) -> ODEAnalysis { - let (matrix, ob_index) = self.build_matrix(model, &data.interaction_coeffs); + ) -> ODEAnalysis> { + let (matrix, ob_index) = self.build_matrix(model); let n = ob_index.len(); - let growth_rates = - ob_index.keys().map(|ob| data.growth_rates.get(ob).copied().unwrap_or_default()); - let b = DVector::from_iterator(n, growth_rates); + let growth_rate_params = ob_index + .keys() + .map(|ob| [(1.0, Monomial::generator(ob.clone()))].into_iter().collect()); + let b = DVector::from_iterator(n, growth_rate_params); let initial_values = ob_index .keys() .map(|ob| data.initial_values.get(ob).copied().unwrap_or_default()); let x0 = DVector::from_iterator(n, initial_values); - let system = LotkaVolterraSystem::new(matrix, b); + let system = lotka_volterra_system(&ob_index.clone().into_keys().collect_vec(), matrix, b) + .extend_scalars(|poly| { + poly.eval(|id| { + data.interaction_coeffs + .get(id) + .or(data.growth_rates.get(id)) + .copied() + .unwrap_or_default() + }) + }) + .map(|p| p.normalize()) + .to_numerical(); let problem = ODEProblem::new(system, x0).end_time(data.duration); ODEAnalysis::new(problem, ob_index) } @@ -75,11 +130,15 @@ mod test { use std::rc::Rc; use super::*; + use crate::simulate::ode::textplot_ode_result; + use crate::stdlib; + use crate::stdlib::analyses::ode::Parameter; use crate::{one::Path, zero::name}; - use crate::{simulate::ode::lotka_volterra, stdlib}; + use expect_test::expect; + use itertools::Itertools; + use nalgebra::{dmatrix, dvector}; - #[test] - fn predator_prey() { + fn predator_prey_from_theory() -> ODEProblem> { let th = Rc::new(stdlib::theories::th_signed_category()); let neg_feedback = stdlib::models::negative_feedback(th); @@ -91,10 +150,108 @@ mod test { initial_values: [(name("x"), 1.0), (name("y"), 1.0)].into_iter().collect(), duration: 10.0, }; - let analysis = SignedCoefficientBuilder::new(name("Object")) + SignedCoefficientBuilder::new(name("Object")) .add_positive(Path::Id(name("Object"))) .add_negative(Path::single(name("Negative"))) - .lotka_volterra_analysis(&neg_feedback, data); - assert_eq!(analysis.problem, lotka_volterra::create_predator_prey()); + .lotka_volterra_analysis(&neg_feedback, data) + .problem + } + + fn predator_prey_from_matrix() -> ODEProblem> { + ODEProblem::new(matrix_example().to_numerical(), dvector![1.0, 1.0]).end_time(10.0) + } + fn matrix_symb_coeff_example() -> PolynomialSystem, u8> + { + let interaction_coeffs = dmatrix!["a11", "a12"; + "a21", "a22"] + .map(|v| [(1.0, Monomial::generator(QualifiedName::from([v])))].into_iter().collect()); + let growth_rates = dvector!["r1", "r2"] + .map(|v| [(1.0, Monomial::generator(QualifiedName::from([v])))].into_iter().collect()); + lotka_volterra_system( + &vec!["x", "y"].into_iter().map(|v| QualifiedName::from([v])).collect_vec(), + interaction_coeffs, + growth_rates, + ) + } + fn matrix_example() -> PolynomialSystem { + let coeffs: HashMap<_, _> = [("a12", -1.0), ("a21", 1.0), ("r1", 2.0), ("r2", -1.0)] + .into_iter() + .map(|(n, v)| (QualifiedName::from([n]), v)) + .collect(); + matrix_symb_coeff_example() + .extend_scalars(|coeff| coeff.eval(|v| coeffs.get(v).copied().unwrap_or_default())) + .map(|p| p.normalize()) + } + + #[test] + fn matrix_agrees_with_theory() { + assert_eq!(predator_prey_from_theory(), predator_prey_from_matrix()); + } + + #[test] + fn lv_solve() { + let problem = predator_prey_from_matrix(); + let result = problem.solve_rk4(0.1).unwrap(); + let expected = expect![[" + ⡁⠀⠀⠀⠀⠀⠀⠀⢠⠊⢢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ 3.5 + ⠄⠀⠀⠀⠀⠀⠀⠀⡇⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡜⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠇⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠀⠀⠀⠀⠀⡎⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⠀⠀⠀⢀⠇⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠇⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠀⠀⠀⠀⡎⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡎⠀⠀⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⠀⠀⣸⡀⠀⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡇⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠀⠀⡎⡜⢣⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡰⢹⠸⡀⠀⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⡸⠀⡇⠈⡆⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⢰⠁⡜⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠄ + ⠂⠀⢠⠃⢸⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⡎⢀⠇⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⡸⠀ + ⡁⠀⡎⠀⡎⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⡸⠀⡸⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⢠⠃⠀ + ⠄⢰⠁⢰⠁⠀⠀⠀⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⢀⠇⢀⠇⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⠀⡎⠀⠀ + ⢂⠇⢀⠇⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⡜⠀⡜⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢄⠀⠀⠀⢰⠁⡰⠁ + ⡝⡠⠊⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⣀⢰⣁⠜⠀⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⣀⢀⢇⡰⠁⠀ + ⠍⠀⠀⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠋⠀⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡏⠁⠀⠀⠀ + ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠸⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡜⠀⠀⠀⠀⠀ + ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢆⠀⠀⠀⠀⠀⠀⠀⠀⡠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⡀⠀⠀⢀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⡀⠀⠀⢀⡠⠔⠁⠀⠀⠀⠀⠀⠀⠀ + ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠉⠉⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ 0.4 + 0.0 10.0 + "]]; + expected.assert_eq(&textplot_ode_result(&problem, &result)); + } + + #[test] + fn latex_symbolic() { + let expected = expect![[r#" + $$ + \begin{align*} + \frac{\mathrm{d}}{\mathrm{d}t} x &= (r_{1}) x + (a_{12}) x y + (a_{11}) x^2\\ + \frac{\mathrm{d}}{\mathrm{d}t} y &= (a_{21}) x y + (r_{2}) y + (a_{22}) y^2 + \end{align*} + $$ + "#]]; + expected.assert_eq( + &matrix_symb_coeff_example() + .extend_scalars(|p| { + p.map_variables(|n| { + let s = n.to_string(); + let (a, b) = s.split_at(1); + QualifiedName::from(format!("{}_{{{}}}", a, b).as_ref()) + }) + }) + .to_latex_string(), + ); + } + + #[test] + fn latex_numerical() { + let expected = expect![[r#" + $$ + \begin{align*} + \frac{\mathrm{d}}{\mathrm{d}t} x &= 2 x + -x y\\ + \frac{\mathrm{d}}{\mathrm{d}t} y &= x y + -y + \end{align*} + $$ + "#]]; + expected.assert_eq(&matrix_example().to_latex_string()); } } diff --git a/packages/catlog/src/stdlib/analyses/ode/signed_coefficients.rs b/packages/catlog/src/stdlib/analyses/ode/signed_coefficients.rs index fd9b605b2..1efae5556 100644 --- a/packages/catlog/src/stdlib/analyses/ode/signed_coefficients.rs +++ b/packages/catlog/src/stdlib/analyses/ode/signed_coefficients.rs @@ -1,11 +1,15 @@ //! Helper module to build analyses based on signed coefficient matrices. -use std::{collections::HashMap, hash::Hash}; +use std::{fmt::Debug, hash::Hash}; use indexmap::IndexMap; use nalgebra::DMatrix; +use num_traits::Zero; -use crate::dbl::model::FgDblModel; +use crate::{ + dbl::model::FgDblModel, + zero::{alg::Polynomial, rig::Monomial}, +}; /// Builder for signed coefficient matrices and analyses based on them. /// @@ -17,6 +21,9 @@ pub struct SignedCoefficientBuilder { negative_mor_types: Vec, } +/// Symbolic parameter in polynomial system. +pub type Parameter = Polynomial; + impl SignedCoefficientBuilder { /// Creates a new builder for the given object type. pub fn new(var_ob_type: ObType) -> Self { @@ -39,17 +46,16 @@ impl SignedCoefficientBuilder { self } - /// Builds the matrix of coefficients for the given model. + /// Builds the matrix of symbolic coefficients for the given model. /// /// Returns the coefficient matrix along with an ordered map from object /// generators to integer indices. pub fn build_matrix( &self, model: &impl FgDblModel, - coeffs: &HashMap, - ) -> (DMatrix, IndexMap) + ) -> (DMatrix>, IndexMap) where - Id: Eq + Clone + Hash + Ord, + Id: Eq + Clone + Hash + Ord + Debug + 'static, { let ob_index: IndexMap<_, _> = model .ob_generators_with_type(&self.var_ob_type) @@ -58,19 +64,19 @@ impl SignedCoefficientBuilder { .collect(); let n = ob_index.len(); - let mut mat = DMatrix::from_element(n, n, 0.0f32); + let mut mat = DMatrix::from_element(n, n, Parameter::::zero()); for mor_type in self.positive_mor_types.iter() { for mor in model.mor_generators_with_type(mor_type) { let i = *ob_index.get(&model.mor_generator_dom(&mor)).unwrap(); let j = *ob_index.get(&model.mor_generator_cod(&mor)).unwrap(); - mat[(j, i)] += coeffs.get(&mor).copied().unwrap_or_default(); + mat[(j, i)] += (1.0, Monomial::generator(mor)); } } for mor_type in self.negative_mor_types.iter() { for mor in model.mor_generators_with_type(mor_type) { let i = *ob_index.get(&model.mor_generator_dom(&mor)).unwrap(); let j = *ob_index.get(&model.mor_generator_cod(&mor)).unwrap(); - mat[(j, i)] -= coeffs.get(&mor).copied().unwrap_or_default(); + mat[(j, i)] += (-1.0, Monomial::generator(mor)); } } diff --git a/packages/catlog/src/zero/alg.rs b/packages/catlog/src/zero/alg.rs index c994653cb..8c422cf72 100644 --- a/packages/catlog/src/zero/alg.rs +++ b/packages/catlog/src/zero/alg.rs @@ -1,9 +1,9 @@ //! Commutative algebra and polynomials. -use num_traits::{One, Pow, Zero}; +use num_traits::{One, Pow, Zero, one, zero}; use std::collections::BTreeMap; use std::fmt::Display; -use std::iter::{Product, Sum}; +use std::iter::{self, Product, Sum}; use std::ops::{Add, AddAssign, Mul, Neg}; use derivative::Derivative; @@ -35,7 +35,7 @@ pub trait CommAlg: CommRing + Module { /// In abstract terms, polynomials with coefficients valued in a [commutative /// ring](super::rig::CommRing) *R* are the free [commutative algebra](CommAlg) /// over *R*. -#[derive(Clone, PartialEq, Eq, Derivative)] +#[derive(Clone, PartialEq, Eq, Derivative, Debug)] #[derivative(Default(bound = ""))] pub struct Polynomial(Combination, Coef>); @@ -183,6 +183,29 @@ where // XXX: Lots of boilerplate to delegate the module structure of `Polynomial` to // `Combination` because these traits cannot be derived automatically. +impl Sum for Polynomial +where + Var: Ord, + Coef: AdditiveMonoid, + Exp: Ord, +{ + fn sum>(iter: I) -> Self { + iter.fold(zero(), |acc, x| acc + x) + } +} + +impl Add for Polynomial +where + Var: Ord, + Coef: Add, + Exp: Ord + Zero, +{ + type Output = Polynomial; + fn add(self, a: Coef) -> Self::Output { + Polynomial(self.0 + iter::once((a, one::>())).collect()) + } +} + impl AddAssign<(Coef, Monomial)> for Polynomial where Var: Ord, @@ -210,7 +233,7 @@ where impl Zero for Polynomial where Var: Ord, - Coef: Add + Zero, + Coef: Zero, Exp: Ord, { fn zero() -> Self {