From 10fcf6c0fb375aa626888774cf4039441b591717 Mon Sep 17 00:00:00 2001 From: George Thomas Date: Wed, 4 Feb 2026 13:08:36 +0000 Subject: [PATCH 1/7] Derive some useful instances --- packages/catlog/src/simulate/ode/polynomial.rs | 1 + packages/catlog/src/zero/alg.rs | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/packages/catlog/src/simulate/ode/polynomial.rs b/packages/catlog/src/simulate/ode/polynomial.rs index 98c5a71dc..518371987 100644 --- a/packages/catlog/src/simulate/ode/polynomial.rs +++ b/packages/catlog/src/simulate/ode/polynomial.rs @@ -184,6 +184,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/zero/alg.rs b/packages/catlog/src/zero/alg.rs index c994653cb..ec218c1a2 100644 --- a/packages/catlog/src/zero/alg.rs +++ b/packages/catlog/src/zero/alg.rs @@ -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>); From 77a436d6a5838e37230c5a80e3bc64774cd4a598 Mon Sep 17 00:00:00 2001 From: George Thomas Date: Wed, 4 Feb 2026 16:29:48 +0000 Subject: [PATCH 2/7] Remove redundant constraint --- packages/catlog/src/zero/alg.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/catlog/src/zero/alg.rs b/packages/catlog/src/zero/alg.rs index ec218c1a2..0687b4ccc 100644 --- a/packages/catlog/src/zero/alg.rs +++ b/packages/catlog/src/zero/alg.rs @@ -210,7 +210,7 @@ where impl Zero for Polynomial where Var: Ord, - Coef: Add + Zero, + Coef: Zero, Exp: Ord, { fn zero() -> Self { From 64cec9364e5744afdd686a1e2dba86f76112837f Mon Sep 17 00:00:00 2001 From: George Thomas Date: Wed, 4 Feb 2026 16:30:05 +0000 Subject: [PATCH 3/7] Add more useful instances --- packages/catlog/src/zero/alg.rs | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/packages/catlog/src/zero/alg.rs b/packages/catlog/src/zero/alg.rs index 0687b4ccc..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; @@ -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, From 122e0f981e8a093b2d3a67264ae73e30b937402a Mon Sep 17 00:00:00 2001 From: George Thomas Date: Wed, 4 Feb 2026 16:30:50 +0000 Subject: [PATCH 4/7] Convert Lotka-Voltera ODEs to polynomials --- .../catlog/src/simulate/ode/lotka_volterra.rs | 27 +++++++++++++++++-- .../src/stdlib/analyses/ode/lotka_volterra.rs | 6 ++--- 2 files changed, 28 insertions(+), 5 deletions(-) diff --git a/packages/catlog/src/simulate/ode/lotka_volterra.rs b/packages/catlog/src/simulate/ode/lotka_volterra.rs index 82561725a..e31bb761c 100644 --- a/packages/catlog/src/simulate/ode/lotka_volterra.rs +++ b/packages/catlog/src/simulate/ode/lotka_volterra.rs @@ -1,5 +1,7 @@ //! Lotka-Volterra differential equations. +use crate::simulate::ode::NumericalPolynomialSystem; +use crate::zero::alg::Polynomial; use nalgebra::{DMatrix, DVector}; #[cfg(test)] @@ -22,6 +24,27 @@ impl LotkaVolterraSystem { pub fn new(A: DMatrix, b: DVector) -> Self { Self { interaction_coeffs: A, growth_rates: b } } + + /// Converts to a numerical polynomial system. + pub fn to_polynomial(self) -> NumericalPolynomialSystem { + NumericalPolynomialSystem { + components: self + .interaction_coeffs + .row_iter() + .enumerate() + .zip(self.growth_rates.into_iter()) + .map(|((i, row), rate)| { + Polynomial::<_, f32, _>::generator(i) + * (row + .iter() + .enumerate() + .map(|(j, a)| Polynomial::generator(j) * *a) + .sum::>() + + *rate) + }) + .collect(), + } + } } impl ODESystem for LotkaVolterraSystem { @@ -33,10 +56,10 @@ impl ODESystem for LotkaVolterraSystem { } #[cfg(test)] -pub(crate) fn create_predator_prey() -> ODEProblem { +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 sys = LotkaVolterraSystem::new(A, b).to_polynomial(); let x0 = DVector::from_column_slice(&[1.0, 1.0]); ODEProblem::new(sys, x0).end_time(10.0) diff --git a/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs b/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs index 323e32701..f68fa10b2 100644 --- a/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs +++ b/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs @@ -14,7 +14,7 @@ use tsify::Tsify; use super::{ODEAnalysis, SignedCoefficientBuilder}; use crate::dbl::model::DiscreteDblModel; -use crate::simulate::ode::{LotkaVolterraSystem, ODEProblem}; +use crate::simulate::ode::{LotkaVolterraSystem, NumericalPolynomialSystem, ODEProblem}; use crate::{one::QualifiedPath, zero::QualifiedName}; /// Data defining a Lotka-Volterra ODE problem for a model. @@ -51,7 +51,7 @@ impl SignedCoefficientBuilder { &self, model: &DiscreteDblModel, data: LotkaVolterraProblemData, - ) -> ODEAnalysis { + ) -> ODEAnalysis> { let (matrix, ob_index) = self.build_matrix(model, &data.interaction_coeffs); let n = ob_index.len(); @@ -64,7 +64,7 @@ impl SignedCoefficientBuilder { .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 = LotkaVolterraSystem::new(matrix, b).to_polynomial(); let problem = ODEProblem::new(system, x0).end_time(data.duration); ODEAnalysis::new(problem, ob_index) } From ac9d9deb7ff1919bfbaedb013450fe4959bc5cb2 Mon Sep 17 00:00:00 2001 From: Patrick Aldis Date: Thu, 5 Feb 2026 11:14:47 +0000 Subject: [PATCH 5/7] Convert linear ODEs to polynomials --- .../catlog/src/simulate/ode/linear_ode.rs | 22 +++++++++++++++++-- .../src/stdlib/analyses/ode/linear_ode.rs | 6 ++--- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/packages/catlog/src/simulate/ode/linear_ode.rs b/packages/catlog/src/simulate/ode/linear_ode.rs index 01c74785c..7e9a3cae2 100644 --- a/packages/catlog/src/simulate/ode/linear_ode.rs +++ b/packages/catlog/src/simulate/ode/linear_ode.rs @@ -2,6 +2,8 @@ use nalgebra::{DMatrix, DVector}; +use crate::{simulate::ode::NumericalPolynomialSystem, zero::alg::Polynomial}; + #[cfg(test)] use super::ODEProblem; use super::ODESystem; @@ -20,6 +22,22 @@ impl LinearODESystem { pub fn new(A: DMatrix) -> Self { Self { coefficients: A } } + + /// Converts to a numerical polynomial system. + pub fn to_polynomial(self) -> NumericalPolynomialSystem { + NumericalPolynomialSystem { + components: self + .coefficients + .row_iter() + .map(|row| { + row.iter() + .enumerate() + .map(|(j, a)| Polynomial::generator(j) * *a) + .sum::>() + }) + .collect(), + } + } } impl ODESystem for LinearODESystem { @@ -30,13 +48,13 @@ impl ODESystem for LinearODESystem { } #[cfg(test)] -pub(crate) fn create_neg_loops_pos_connector() -> ODEProblem { +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 system = LinearODESystem::new(A).to_polynomial(); let initial = dvector![2.0, 1.0, 1.0]; ODEProblem::new(system, initial).end_time(10.0) } diff --git a/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs b/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs index 40c06ed17..5125f3f42 100644 --- a/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs +++ b/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs @@ -14,7 +14,7 @@ use tsify::Tsify; use super::{ODEAnalysis, SignedCoefficientBuilder}; use crate::dbl::model::DiscreteDblModel; -use crate::simulate::ode::{LinearODESystem, ODEProblem}; +use crate::simulate::ode::{LinearODESystem, NumericalPolynomialSystem, ODEProblem}; use crate::{one::QualifiedPath, zero::QualifiedName}; /// Data defining a linear ODE problem for a model. @@ -47,7 +47,7 @@ impl SignedCoefficientBuilder { &self, model: &DiscreteDblModel, data: LinearODEProblemData, - ) -> ODEAnalysis { + ) -> ODEAnalysis> { let (matrix, ob_index) = self.build_matrix(model, &data.coefficients); let n = ob_index.len(); @@ -56,7 +56,7 @@ 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 = LinearODESystem::new(matrix).to_polynomial(); let problem = ODEProblem::new(system, x0).end_time(data.duration); ODEAnalysis::new(problem, ob_index) } From 3002eb0e4a5eae6bd9be8f7f9f741302ccecfb0a Mon Sep 17 00:00:00 2001 From: George Thomas Date: Thu, 5 Feb 2026 11:41:16 +0000 Subject: [PATCH 6/7] Remove Lotka-Voltera type --- .../catlog/src/simulate/ode/lotka_volterra.rs | 61 ++++++------------- .../src/stdlib/analyses/ode/lotka_volterra.rs | 5 +- 2 files changed, 21 insertions(+), 45 deletions(-) diff --git a/packages/catlog/src/simulate/ode/lotka_volterra.rs b/packages/catlog/src/simulate/ode/lotka_volterra.rs index e31bb761c..b041ab295 100644 --- a/packages/catlog/src/simulate/ode/lotka_volterra.rs +++ b/packages/catlog/src/simulate/ode/lotka_volterra.rs @@ -6,52 +6,31 @@ use nalgebra::{DMatrix, DVector}; #[cfg(test)] use super::ODEProblem; -use super::ODESystem; -/// A Lotka-Volterra dynamical system. +/// 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). -#[derive(Clone, Debug, PartialEq)] -pub struct LotkaVolterraSystem { +pub fn lotka_volterra_system( 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 } - } - - /// Converts to a numerical polynomial system. - pub fn to_polynomial(self) -> NumericalPolynomialSystem { - NumericalPolynomialSystem { - components: self - .interaction_coeffs - .row_iter() - .enumerate() - .zip(self.growth_rates.into_iter()) - .map(|((i, row), rate)| { - Polynomial::<_, f32, _>::generator(i) - * (row - .iter() - .enumerate() - .map(|(j, a)| Polynomial::generator(j) * *a) - .sum::>() - + *rate) - }) - .collect(), - } - } -} - -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); +) -> NumericalPolynomialSystem { + NumericalPolynomialSystem { + components: interaction_coeffs + .row_iter() + .enumerate() + .zip(growth_rates.into_iter()) + .map(|((i, row), rate)| { + Polynomial::<_, f32, _>::generator(i) + * (row + .iter() + .enumerate() + .map(|(j, a)| Polynomial::generator(j) * *a) + .sum::>() + + *rate) + }) + .collect(), } } @@ -59,10 +38,8 @@ impl ODESystem for LotkaVolterraSystem { 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).to_polynomial(); - let x0 = DVector::from_column_slice(&[1.0, 1.0]); - ODEProblem::new(sys, x0).end_time(10.0) + ODEProblem::new(lotka_volterra_system(A, b), x0).end_time(10.0) } #[cfg(test)] diff --git a/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs b/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs index f68fa10b2..dd61c5ce2 100644 --- a/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs +++ b/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs @@ -14,7 +14,7 @@ use tsify::Tsify; use super::{ODEAnalysis, SignedCoefficientBuilder}; use crate::dbl::model::DiscreteDblModel; -use crate::simulate::ode::{LotkaVolterraSystem, NumericalPolynomialSystem, ODEProblem}; +use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, lotka_volterra_system}; use crate::{one::QualifiedPath, zero::QualifiedName}; /// Data defining a Lotka-Volterra ODE problem for a model. @@ -64,8 +64,7 @@ impl SignedCoefficientBuilder { .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).to_polynomial(); - let problem = ODEProblem::new(system, x0).end_time(data.duration); + let problem = ODEProblem::new(lotka_volterra_system(matrix, b), x0).end_time(data.duration); ODEAnalysis::new(problem, ob_index) } } From f2ae653841f1cf30899f542e99daae5a10587fff Mon Sep 17 00:00:00 2001 From: George Thomas Date: Thu, 5 Feb 2026 11:43:34 +0000 Subject: [PATCH 7/7] Remove Linear ODE type --- .../catlog/src/simulate/ode/linear_ode.rs | 51 +++++-------------- .../src/stdlib/analyses/ode/linear_ode.rs | 4 +- 2 files changed, 16 insertions(+), 39 deletions(-) diff --git a/packages/catlog/src/simulate/ode/linear_ode.rs b/packages/catlog/src/simulate/ode/linear_ode.rs index 7e9a3cae2..12eeed43a 100644 --- a/packages/catlog/src/simulate/ode/linear_ode.rs +++ b/packages/catlog/src/simulate/ode/linear_ode.rs @@ -1,49 +1,27 @@ //! Constant-coefficient linear first-order differential equations. -use nalgebra::{DMatrix, DVector}; +use nalgebra::DMatrix; use crate::{simulate::ode::NumericalPolynomialSystem, zero::alg::Polynomial}; #[cfg(test)] use super::ODEProblem; -use super::ODESystem; -/// A (constant-coefficient) linear (first-order) dynamical system. +/// Construct 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 } - } - - /// Converts to a numerical polynomial system. - pub fn to_polynomial(self) -> NumericalPolynomialSystem { - NumericalPolynomialSystem { - components: self - .coefficients - .row_iter() - .map(|row| { - row.iter() - .enumerate() - .map(|(j, a)| Polynomial::generator(j) * *a) - .sum::>() - }) - .collect(), - } - } -} - -impl ODESystem for LinearODESystem { - fn vector_field(&self, dx: &mut DVector, x: &DVector, _t: f32) { - let A = &self.coefficients; - *dx = A * x +pub fn linear_polynomial_system(coefficients: DMatrix) -> NumericalPolynomialSystem { + NumericalPolynomialSystem { + components: coefficients + .row_iter() + .map(|row| { + row.iter() + .enumerate() + .map(|(j, a)| Polynomial::generator(j) * *a) + .sum::>() + }) + .collect(), } } @@ -54,9 +32,8 @@ pub(crate) fn create_neg_loops_pos_connector() -> ODEProblem { .map(|ob| data.initial_values.get(ob).copied().unwrap_or_default()); let x0 = DVector::from_iterator(n, initial_values); - let system = LinearODESystem::new(matrix).to_polynomial(); + let system = linear_polynomial_system(matrix); let problem = ODEProblem::new(system, x0).end_time(data.duration); ODEAnalysis::new(problem, ob_index) }