Skip to content

Factor linear/Lotka-Voltera ODE analyses through PolynomialSystem #1007

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
georgefst wants to merge 7 commits into
base: main
Choose a base branch
from
Open

Factor linear/Lotka-Voltera ODE analyses through PolynomialSystem #1007

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
10fcf6c
Derive some useful instances
georgefst Feb 4, 2026
77a436d
Remove redundant constraint
georgefst Feb 4, 2026
64cec93
Add more useful instances
georgefst Feb 4, 2026
122e0f9
Convert Lotka-Voltera ODEs to polynomials
georgefst Feb 4, 2026
ac9d9de
Convert linear ODEs to polynomials
Feb 5, 2026
3002eb0
Remove Lotka-Voltera type
georgefst Feb 5, 2026
f2ae653
Remove Linear ODE type
georgefst Feb 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
39 changes: 17 additions & 22 deletions packages/catlog/src/simulate/ode/linear_ode.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,44 +1,39 @@
//! 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<f32>,
}

impl LinearODESystem {
/// Constructs a linear ODE system with the given coefficient matrix.
pub fn new(A: DMatrix<f32>) -> Self {
Self { coefficients: A }
}
}

impl ODESystem for LinearODESystem {
fn vector_field(&self, dx: &mut DVector<f32>, x: &DVector<f32>, _t: f32) {
let A = &self.coefficients;
*dx = A * x
pub fn linear_polynomial_system(coefficients: DMatrix<f32>) -> NumericalPolynomialSystem<u8> {
NumericalPolynomialSystem {
components: coefficients
.row_iter()
.map(|row| {
row.iter()
.enumerate()
.map(|(j, a)| Polynomial::generator(j) * *a)
.sum::<Polynomial<_, _, _>>()
})
.collect(),
}
}

#[cfg(test)]
pub(crate) fn create_neg_loops_pos_connector() -> ODEProblem<LinearODESystem> {
pub(crate) fn create_neg_loops_pos_connector() -> ODEProblem<NumericalPolynomialSystem<u8>> {
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)
ODEProblem::new(linear_polynomial_system(A), initial).end_time(10.0)
}

#[cfg(test)]
Expand Down
44 changes: 22 additions & 22 deletions packages/catlog/src/simulate/ode/lotka_volterra.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,45 +1,45 @@
//! Lotka-Volterra differential equations.

use crate::simulate::ode::NumericalPolynomialSystem;
use crate::zero::alg::Polynomial;
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<f32>,
growth_rates: DVector<f32>,
}

impl LotkaVolterraSystem {
/// Constructs a new Lokta-Volterra system with the given parameters.
pub fn new(A: DMatrix<f32>, b: DVector<f32>) -> Self {
Self { interaction_coeffs: A, growth_rates: b }
}
}

impl ODESystem for LotkaVolterraSystem {
fn vector_field(&self, dx: &mut DVector<f32>, x: &DVector<f32>, _t: f32) {
let A = &self.interaction_coeffs;
let b = &self.growth_rates;
*dx = (A * x + b).component_mul(x);
) -> NumericalPolynomialSystem<u8> {
NumericalPolynomialSystem {
components: interaction_coeffs
.row_iter()
.enumerate()
.zip(growth_rates.into_iter())

Check failure on line 23 in packages/catlog/src/simulate/ode/lotka_volterra.rs

View workflow job for this annotation

GitHub Actions / rust lints 

explicit call to `.into_iter()` in function argument accepting `IntoIterator`
Copy link
Member

@epatters epatters Feb 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To avoid triggering the linter in CI, you can run cargo clippy locally before pushing.

.map(|((i, row), rate)| {
Polynomial::<_, f32, _>::generator(i)
* (row
.iter()
.enumerate()
.map(|(j, a)| Polynomial::generator(j) * *a)
.sum::<Polynomial<_, _, _>>()
+ *rate)
})
.collect(),
}
}

#[cfg(test)]
pub(crate) fn create_predator_prey() -> ODEProblem<LotkaVolterraSystem> {
pub(crate) fn create_predator_prey() -> ODEProblem<NumericalPolynomialSystem<u8>> {
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)
ODEProblem::new(lotka_volterra_system(A, b), x0).end_time(10.0)
}

#[cfg(test)]
Expand Down
1 change: 1 addition & 0 deletions packages/catlog/src/simulate/ode/polynomial.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<Exp> {
/// Components of the vector field.
pub components: Vec<Polynomial<usize, f32, Exp>>,
Expand Down
6 changes: 3 additions & 3 deletions packages/catlog/src/stdlib/analyses/ode/linear_ode.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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::{NumericalPolynomialSystem, ODEProblem, linear_polynomial_system};
use crate::{one::QualifiedPath, zero::QualifiedName};

/// Data defining a linear ODE problem for a model.
Expand Down Expand Up @@ -47,7 +47,7 @@ impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
&self,
model: &DiscreteDblModel,
data: LinearODEProblemData,
) -> ODEAnalysis<LinearODESystem> {
) -> ODEAnalysis<NumericalPolynomialSystem<u8>> {
Copy link
Author

@georgefst georgefst Feb 5, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The exponent for linear systems is obviously always 1. So it's tempting to use or define a singleton type, but that might be overkill. Not sure how easy or conventional it is in Rust.

Copy link
Contributor

@patrickaldis patrickaldis Feb 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure this accounts for constants (zero exponents) though. We'd need a type with 0 and 1

Copy link
Author

@georgefst georgefst Feb 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh yeah, good point 🤦‍♂️

let (matrix, ob_index) = self.build_matrix(model, &data.coefficients);
let n = ob_index.len();

Expand All @@ -56,7 +56,7 @@ impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
.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(matrix);
let problem = ODEProblem::new(system, x0).end_time(data.duration);
ODEAnalysis::new(problem, ob_index)
}
Expand Down
7 changes: 3 additions & 4 deletions packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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::{NumericalPolynomialSystem, ODEProblem, lotka_volterra_system};
use crate::{one::QualifiedPath, zero::QualifiedName};

/// Data defining a Lotka-Volterra ODE problem for a model.
Expand Down Expand Up @@ -51,7 +51,7 @@ impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
&self,
model: &DiscreteDblModel,
data: LotkaVolterraProblemData,
) -> ODEAnalysis<LotkaVolterraSystem> {
) -> ODEAnalysis<NumericalPolynomialSystem<u8>> {
let (matrix, ob_index) = self.build_matrix(model, &data.interaction_coeffs);
let n = ob_index.len();

Expand All @@ -64,8 +64,7 @@ impl SignedCoefficientBuilder<QualifiedName, QualifiedPath> {
.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 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)
}
}
Expand Down
31 changes: 27 additions & 4 deletions packages/catlog/src/zero/alg.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;
Expand Down Expand Up @@ -35,7 +35,7 @@ pub trait CommAlg: CommRing + Module<Ring = Self::R> {
/// 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<Var, Coef, Exp>(Combination<Monomial<Var, Exp>, Coef>);

Expand Down Expand Up @@ -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<Var, Coef, Exp> Sum for Polynomial<Var, Coef, Exp>
Copy link
Author

@georgefst georgefst Feb 5, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's nothing here that's actually specific to Polynomial. The exact same code could be used for any AdditiveMonoid. But I don't think Rust allows instances with that sort of flexibility. We could define some reusable helper, but probably not worth it.

Note that Polynomial is an AdditiveMonoid precisely when its coefficient type is, due to the instances:

  • impl<Var, Coef, Exp> Add for Polynomial<Var, Coef, Exp>
  • impl<Var, Coef, Exp> Zero for Polynomial<Var, Coef, Exp>

Copy link
Author

@georgefst georgefst Feb 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, and also, looking at the comment above this, maybe we should actually move the main Impl code to Combination, and just delegate to that here?

where
Var: Ord,
Coef: AdditiveMonoid,
Exp: Ord,
{
fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
iter.fold(zero(), |acc, x| acc + x)
}
}

impl<Var, Coef, Exp> Add<Coef> for Polynomial<Var, Coef, Exp>
where
Var: Ord,
Coef: Add<Output = Coef>,
Exp: Ord + Zero,
{
type Output = Polynomial<Var, Coef, Exp>;
fn add(self, a: Coef) -> Self::Output {
Polynomial(self.0 + iter::once((a, one::<Monomial<_, _>>())).collect())
}
}

impl<Var, Coef, Exp> AddAssign<(Coef, Monomial<Var, Exp>)> for Polynomial<Var, Coef, Exp>
where
Var: Ord,
Expand Down Expand Up @@ -210,7 +233,7 @@ where
impl<Var, Coef, Exp> Zero for Polynomial<Var, Coef, Exp>
where
Var: Ord,
Coef: Add<Output = Coef> + Zero,
Coef: Zero,
Copy link
Author

@georgefst georgefst Feb 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This slight generalisation was useful during debugging. I'm not sure if it was necessary in the end, but there's no harm.

Coming from Haskell, I am a bit surprised that the compiler doesn't warn about redundant constraints.

Copy link
Member

@epatters epatters Feb 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, I've found myself wishing for such warnings but AFAIK, the Rust compiler doesn't warn about either redundant constraints or unnecessary constraints.

Exp: Ord,
{
fn zero() -> Self {
Expand Down
Loading