diff --git a/packages/catlog-wasm/src/analyses.rs b/packages/catlog-wasm/src/analyses.rs index 57fee9811..2b9513433 100644 --- a/packages/catlog-wasm/src/analyses.rs +++ b/packages/catlog-wasm/src/analyses.rs @@ -4,9 +4,11 @@ use serde::{Deserialize, Serialize}; use tsify::Tsify; use catlog::simulate::ode::PolynomialSystem; -use catlog::stdlib::analyses::ode; +use catlog::stdlib::analyses::ode::{self, ODESemanticsAnalysis, ODESemanticsProblemData}; use catlog::zero::QualifiedName; +use crate::latex::{latex_mor_names_linear_ode, latex_mor_names_lotka_volterra}; + use super::latex::{LatexEquations, latex_mor_names, latex_mor_names_mass_action, latex_ob_names}; use super::model::DblModel; use super::result::JsResult; @@ -74,8 +76,8 @@ pub(crate) fn polynomial_ode_simulation( }) } -/// The mass-action analysis is currently implemented for Petri nets and stock-flow -/// diagrams, and we can avoid some code reduplication by making this explicit. +/// Mass-action analysis is currently implemented for Petri nets and stock-flow diagrams +/// and we can avoid some code reduplication by making this explicit. pub enum MassActionAnalysisLogic { /// The modal theory of Petri nets. PetriNet, @@ -88,17 +90,23 @@ fn mass_action_system( model: &DblModel, mass_conservation_type: ode::MassConservationType, logic: MassActionAnalysisLogic, -) -> Result, i8>, String> { +) -> Result, i8>, String> { match logic { MassActionAnalysisLogic::PetriNet => { let realised_model = model.modal_unital()?; - let analysis = ode::PetriNetMassActionAnalysis::default(); - Ok(analysis.build_system(realised_model, mass_conservation_type)) + let analysis = ode::PetriNetMassActionAnalysis { + mass_conservation_type, + ..ode::PetriNetMassActionAnalysis::default() + }; + Ok(analysis.build_system(realised_model)) } MassActionAnalysisLogic::StockFlow => { let realised_model = model.discrete_tab()?; - let analysis = ode::StockFlowMassActionAnalysis::default(); - Ok(analysis.build_system(realised_model, mass_conservation_type)) + let analysis = ode::StockFlowMassActionAnalysis { + mass_conservation_type, + ..ode::StockFlowMassActionAnalysis::default() + }; + Ok(analysis.build_system(realised_model)) } } } @@ -133,10 +141,99 @@ pub(crate) fn mass_action_simulation( logic: MassActionAnalysisLogic, ) -> Result { let sys = mass_action_system(model, data.mass_conservation_type, logic); - let sys_extended_scalars = ode::extend_mass_action_scalars(sys?, &data); + let sys_extended_scalars = data.extend_scalars(sys?); + let latex_equations = + sys_extended_scalars.map_variables(latex_ob_names(model)).to_latex_equations(); + let analysis = data.build_analysis(sys_extended_scalars); + let solution = analysis.solve_with_defaults().map_err(|err| format!("{err:?}")); + Ok(ODEResultWithEquations { + solution: ODEResult(solution.into()), + latex_equations: LatexEquations(latex_equations), + }) +} + +/// Generates the PolynomialSystem for Lotka-Volterra dynamics. +fn lotka_volterra_system( + model: &DblModel, +) -> Result, i8>, String> +{ + let realised_model = model.discrete()?; + let analysis = ode::LotkaVolterraAnalysis::default(); + Ok(analysis.build_system(realised_model)) +} + +/// The analysis data for polynomial ODE equations. +#[derive(Serialize, Deserialize, Tsify)] +#[tsify(into_wasm_abi, from_wasm_abi)] +pub struct LotkaVolterraEquationsData { + #[serde(rename = "trivialData")] + trivial_data: bool, +} + +/// Generates Lotka-Volterra equations for the system. +pub(crate) fn lotka_volterra_equations(model: &DblModel) -> Result { + let sys = lotka_volterra_system(model); + let equations = sys? + .map_variables(latex_ob_names(model)) + .extend_scalars(|param| param.map_variables(latex_mor_names_lotka_volterra(model))) + .to_latex_equations(); + Ok(LatexEquations(equations)) +} + +/// Simulates Lotka-Volterra ODEs. +pub(crate) fn lotka_volterra_simulation( + model: &DblModel, + data: ode::LotkaVolterraProblemData, +) -> Result { + let sys = lotka_volterra_system(model); + let sys_extended_scalars = data.extend_scalars(sys?); + let latex_equations = + sys_extended_scalars.map_variables(latex_ob_names(model)).to_latex_equations(); + let analysis = data.build_analysis(sys_extended_scalars); + let solution = analysis.solve_with_defaults().map_err(|err| format!("{err:?}")); + Ok(ODEResultWithEquations { + solution: ODEResult(solution.into()), + latex_equations: LatexEquations(latex_equations), + }) +} + +/// Generates the PolynomialSystem for linear ODE dynamics. +fn linear_ode_system( + model: &DblModel, +) -> Result, i8>, String> { + let realised_model = model.discrete()?; + let analysis = ode::LinearODEAnalysis::default(); + Ok(analysis.build_system(realised_model)) +} + +/// The analysis data for polynomial ODE equations. +#[derive(Serialize, Deserialize, Tsify)] +#[tsify(into_wasm_abi, from_wasm_abi)] +pub struct LinearODEEquationsData { + #[serde(rename = "trivialData")] + trivial_data: bool, +} + +/// Generates linear ODE equations for the system. +pub(crate) fn linear_ode_equations(model: &DblModel) -> Result { + let sys = linear_ode_system(model); + let equations = sys? + .map_variables(latex_ob_names(model)) + .extend_scalars(|param| param.map_variables(latex_mor_names_linear_ode(model))) + .to_latex_equations(); + Ok(LatexEquations(equations)) +} + +/// Simulates linear ODE equations. +pub(crate) fn linear_ode_simulation( + model: &DblModel, + data: ode::LinearODEProblemData, +) -> Result { + let sys = linear_ode_system(model); + let sys_extended_scalars = data.extend_scalars(sys?); let latex_equations = sys_extended_scalars.map_variables(latex_ob_names(model)).to_latex_equations(); - let analysis = ode::into_mass_action_analysis(sys_extended_scalars, data); + let analysis = data.build_analysis(sys_extended_scalars); let solution = analysis.solve_with_defaults().map_err(|err| format!("{err:?}")); Ok(ODEResultWithEquations { solution: ODEResult(solution.into()), diff --git a/packages/catlog-wasm/src/latex.rs b/packages/catlog-wasm/src/latex.rs index 848494f22..772df0785 100644 --- a/packages/catlog-wasm/src/latex.rs +++ b/packages/catlog-wasm/src/latex.rs @@ -57,7 +57,7 @@ pub(crate) fn latex_mor_names(model: &DblModel) -> impl Fn(&QualifiedName) -> St /// falls back to the domain→codomain format (e.g., `X \to Y`). pub(crate) fn latex_mor_names_mass_action( model: &DblModel, -) -> impl Fn(&ode::FlowParameter) -> String { +) -> impl Fn(&ode::MassActionParameter) -> String { // Returns a LaTeX fragment for a transition, suitable for use as a subscript. // Named morphisms produce `\text{name}`, unnamed ones produce // `\text{dom} \to \text{cod}` so that `\to` is in math mode. @@ -72,31 +72,106 @@ pub(crate) fn latex_mor_names_mass_action( } }; - move |id: &ode::FlowParameter| match id { - ode::FlowParameter::Balanced { transition } => { + move |id: &ode::MassActionParameter| match id { + ode::MassActionParameter::Balanced { flow: transition } => { let sub = transition_subscript(transition); format!("r_{{{sub}}}") } - ode::FlowParameter::Unbalanced { direction, parameter } => match (direction, parameter) { - (ode::Direction::IncomingFlow, ode::RateParameter::PerTransition { transition }) => { - let sub = transition_subscript(transition); - format!("\\rho_{{{sub}}}") + ode::MassActionParameter::Unbalanced { direction, parameter } => { + match (direction, parameter) { + ( + ode::Direction::IncomingFlow, + ode::RateParameter::PerTransition { flow: transition }, + ) => { + let sub = transition_subscript(transition); + format!("\\rho_{{{sub}}}") + } + ( + ode::Direction::OutgoingFlow, + ode::RateParameter::PerTransition { flow: transition }, + ) => { + let sub = transition_subscript(transition); + format!("\\kappa_{{{sub}}}") + } + ( + ode::Direction::IncomingFlow, + ode::RateParameter::PerPlace { flow: transition, stock: place }, + ) => { + let sub = transition_subscript(transition); + let output_place_label = model.ob_namespace.label_string(place); + format!("\\rho_{{{sub}}}^{{\\text{{{output_place_label}}}}}") + } + ( + ode::Direction::OutgoingFlow, + ode::RateParameter::PerPlace { flow: transition, stock: place }, + ) => { + let sub = transition_subscript(transition); + let input_place_label = model.ob_namespace.label_string(place); + format!("\\kappa_{{{sub}}}^{{\\text{{{input_place_label}}}}}") + } } - (ode::Direction::OutgoingFlow, ode::RateParameter::PerTransition { transition }) => { - let sub = transition_subscript(transition); - format!("\\kappa_{{{sub}}}") - } - (ode::Direction::IncomingFlow, ode::RateParameter::PerPlace { transition, place }) => { - let sub = transition_subscript(transition); - let output_place_label = model.ob_namespace.label_string(place); - format!("\\rho_{{{sub}}}^{{\\text{{{output_place_label}}}}}") - } - (ode::Direction::OutgoingFlow, ode::RateParameter::PerPlace { transition, place }) => { - let sub = transition_subscript(transition); - let input_place_label = model.ob_namespace.label_string(place); - format!("\\kappa_{{{sub}}}^{{\\text{{{input_place_label}}}}}") - } - }, + } + } +} + +/// Creates a closure that formats morphism names for Lotka-Volterra LaTeX output. +/// +/// When a morphism has a label, it is used directly. When unnamed, the label +/// falls back to the domain→codomain format (e.g., `X \to Y`). +pub(crate) fn latex_mor_names_lotka_volterra( + model: &DblModel, +) -> impl Fn(&ode::LotkaVolterraParameter) -> String { + // Returns a LaTeX fragment for a transition, suitable for use as a subscript. + // Named morphisms produce `\text{name}`, unnamed ones produce + // `\text{dom} \to \text{cod}` so that `\to` is in math mode. + let transition_subscript = |transition: &QualifiedName| -> String { + if let Some(label) = model.mor_namespace.label(transition) { + format!("\\text{{{label}}}") + } else { + let (dom, cod) = model + .mor_generator_dom_cod_label_strings(transition) + .expect("Morphism in equation system should have domain and codomain"); + format!("\\text{{{dom}}} \\to \\text{{{cod}}}") + } + }; + + move |id: &ode::LotkaVolterraParameter| match id { + ode::LotkaVolterraParameter::Growth { variable } => { + format!("g_{{{variable}}}") + } + ode::LotkaVolterraParameter::Interaction { link } => { + let sub = transition_subscript(link); + format!("k_{{{sub}}}") + } + } +} + +/// Creates a closure that formats morphism names for mass-action LaTeX output. +/// +/// When a morphism has a label, it is used directly. When unnamed, the label +/// falls back to the domain→codomain format (e.g., `X \to Y`). +pub(crate) fn latex_mor_names_linear_ode( + model: &DblModel, +) -> impl Fn(&ode::LinearODEParameter) -> String { + // Returns a LaTeX fragment for a transition, suitable for use as a subscript. + // Named morphisms produce `\text{name}`, unnamed ones produce + // `\text{dom} \to \text{cod}` so that `\to` is in math mode. + let transition_subscript = |transition: &QualifiedName| -> String { + if let Some(label) = model.mor_namespace.label(transition) { + format!("\\text{{{label}}}") + } else { + let (dom, cod) = model + .mor_generator_dom_cod_label_strings(transition) + .expect("Morphism in equation system should have domain and codomain"); + format!("\\text{{{dom}}} \\to \\text{{{cod}}}") + } + }; + + move |id: &ode::LinearODEParameter| match id { + ode::LinearODEParameter::Parameter { morphism } => { + let sub = transition_subscript(morphism); + format!("\\lambda_{{{sub}}}") + } } } @@ -105,6 +180,7 @@ mod tests { use catlog::dbl::modal::{List, ModalMorType, ModalOb, ModalObType}; use catlog::dbl::model::{ModalDblModel, MutDblModel}; use catlog::simulate::ode::LatexEquation; + use catlog::stdlib::analyses::ode::{StockFlowMassActionAnalysis, ode_semantics::*}; use catlog::stdlib::{analyses::ode, theories}; use catlog::zero::{LabelSegment, Namespace, QualifiedName}; use std::rc::Rc; @@ -117,11 +193,13 @@ mod tests { fn unbalanced_mass_action_latex_equations() { let model = backward_link("xxx", "yyy", "fff"); let tab_model = model.discrete_tab().unwrap(); - let analysis = ode::StockFlowMassActionAnalysis::default(); - let sys = analysis.build_system( - tab_model, - ode::MassConservationType::Unbalanced(ode::RateGranularity::PerTransition), - ); + let analysis = StockFlowMassActionAnalysis { + mass_conservation_type: ode::MassConservationType::Unbalanced( + ode::RateGranularity::PerTransition, + ), + ..StockFlowMassActionAnalysis::default() + }; + let sys = analysis.build_system(tab_model); let equations = sys .map_variables(latex_ob_names(&model)) .extend_scalars(|param| param.map_variables(latex_mor_names_mass_action(&model))) @@ -144,11 +222,13 @@ mod tests { fn unnamed_mor_uses_dom_cod_in_equations() { let model = backward_link("xxx", "yyy", ""); let tab_model = model.discrete_tab().unwrap(); - let analysis = ode::StockFlowMassActionAnalysis::default(); - let sys = analysis.build_system( - tab_model, - ode::MassConservationType::Unbalanced(ode::RateGranularity::PerTransition), - ); + let analysis = StockFlowMassActionAnalysis { + mass_conservation_type: ode::MassConservationType::Unbalanced( + ode::RateGranularity::PerTransition, + ), + ..StockFlowMassActionAnalysis::default() + }; + let sys = analysis.build_system(tab_model); let equations = sys .map_variables(latex_ob_names(&model)) .extend_scalars(|param| param.map_variables(latex_mor_names_mass_action(&model))) diff --git a/packages/catlog-wasm/src/theories.rs b/packages/catlog-wasm/src/theories.rs index f01b6c7a1..a2b0c1ad0 100644 --- a/packages/catlog-wasm/src/theories.rs +++ b/packages/catlog-wasm/src/theories.rs @@ -150,16 +150,14 @@ impl ThSignedCategory { &self, model: &DblModel, data: analyses::ode::LotkaVolterraProblemData, - ) -> Result { - Ok(ODEResult( - analyses::ode::SignedCoefficientBuilder::new(name("Object")) - .add_positive(Path::Id(name("Object"))) - .add_negative(name("Negative").into()) - .lotka_volterra_analysis(model.discrete()?, data) - .solve_with_defaults() - .map_err(|err| format!("{err:?}")) - .into(), - )) + ) -> Result { + lotka_volterra_simulation(model, data) + } + + /// Show the equations of the Lotka-Volterra system derived from a model. + #[wasm_bindgen(js_name = "lotkaVolterraEquations")] + pub fn lotka_volterra_equations(&self, model: &DblModel) -> Result { + lotka_volterra_equations(model) } /// Simulate the linear ODE system derived from a model. @@ -168,16 +166,14 @@ impl ThSignedCategory { &self, model: &DblModel, data: analyses::ode::LinearODEProblemData, - ) -> Result { - Ok(ODEResult( - analyses::ode::SignedCoefficientBuilder::new(name("Object")) - .add_positive(Path::Id(name("Object"))) - .add_negative(name("Negative").into()) - .linear_ode_analysis(model.discrete()?, data) - .solve_with_defaults() - .map_err(|err| format!("{err:?}")) - .into(), - )) + ) -> Result { + linear_ode_simulation(model, data) + } + + /// Show the equations of the linear ODE system derived from a model. + #[wasm_bindgen(js_name = "linearODEEquations")] + pub fn linear_ode_equations(&self, model: &DblModel) -> Result { + linear_ode_equations(model) } } @@ -340,26 +336,6 @@ impl ThCategorySignedLinks { pub fn theory(&self) -> DblTheory { DblTheory(self.0.clone().into()) } - - /// Simulates the mass-action ODE system derived from a model. - #[wasm_bindgen(js_name = "massAction")] - pub fn mass_action( - &self, - model: &DblModel, - data: analyses::ode::MassActionProblemData, - ) -> Result { - mass_action_simulation(model, data, MassActionAnalysisLogic::StockFlow) - } - - /// Returns the symbolic mass-action equations in LaTeX format. - #[wasm_bindgen(js_name = "massActionEquations")] - pub fn mass_action_equations( - &self, - model: &DblModel, - data: MassActionEquationsData, - ) -> Result { - mass_action_equations(model, data, MassActionAnalysisLogic::StockFlow) - } } /// The theory of strict symmetric monoidal categories. diff --git a/packages/catlog/src/stdlib/analyses/mod.rs b/packages/catlog/src/stdlib/analyses/mod.rs index 861f3ef8c..2fffc0f45 100644 --- a/packages/catlog/src/stdlib/analyses/mod.rs +++ b/packages/catlog/src/stdlib/analyses/mod.rs @@ -1,6 +1,7 @@ //! Various analyses that can be performed on models. pub(crate) mod petri; +pub(crate) mod stock_flow; #[cfg(feature = "ode")] pub mod ode; diff --git a/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs b/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs index 83f4a6c22..2b6a8f212 100644 --- a/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs +++ b/packages/catlog/src/stdlib/analyses/ode/linear_ode.rs @@ -1,29 +1,141 @@ -//! Constant-coefficient linear first-order ODE analysis of models. +//! Linear constant-coefficient first-order ODE analysis of models. //! -//! The main entry point for this module is -//! [`linear_ode_analysis`](SignedCoefficientBuilder::linear_ode_analysis). +//! This follows the structure of [`ode::ode_semantics`], implementing `ODESemantics` for the struct +//! `LinearODESemantics`. +//! +//! [`ode::ode_semantics`]: crate::stdlib::analyses::ode::ode_semantics use std::collections::HashMap; -use std::hash::Hash; -use std::ops::Add; - -use indexmap::IndexMap; -use itertools::Itertools; -use nalgebra::{DMatrix, DVector}; -use num_traits::Zero; +use std::fmt; #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; #[cfg(feature = "serde-wasm")] use tsify::Tsify; -use super::{ODEAnalysis, Parameter, SignedCoefficientBuilder}; -use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem}; -use crate::{ - dbl::model::DiscreteDblModel, - one::QualifiedPath, - zero::{QualifiedName, rig::Monomial}, +use super::Parameter; +use crate::dbl::model::{FpDblModel, MutDblModel}; +use crate::one::Path; +use crate::simulate::ode::PolynomialSystem; +use crate::stdlib::analyses::ode::ode_semantics::{ + ContributionSign, ODEParameterType, ODESemantics, ODESemanticsAnalysis, + ODESemanticsProblemData, PolynomialODESystemBuilder, }; +use crate::zero::name; +use crate::{dbl::model::DiscreteDblModel, one::QualifiedPath, zero::QualifiedName}; + +/// Implementing LinearODE as an ODE semantics for models of type `DiscreteDblModel`. +pub struct LinearODESemantics; + +impl ODESemantics for LinearODESemantics { + type ModelType = DiscreteDblModel; + type ParameterType = LinearODEParameter; + type AnalysisType = LinearODEAnalysis; + type ProblemDataType = LinearODEProblemData; +} + +/// Parameters in the linear equations correspond only to morphisms. +#[derive(PartialEq, Eq, PartialOrd, Ord, Clone)] +pub enum LinearODEParameter { + /// The parameter associated to a morphism. + Parameter { + /// The morphism. + morphism: QualifiedName, + }, +} + +impl fmt::Display for LinearODEParameter { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + Self::Parameter { morphism } => { + write!(f, "Parameter({})", morphism) + } + } + } +} + +impl ODEParameterType for LinearODEParameter {} + +/// Linear ODE analysis for causal loop diagrams (CLDs). +pub struct LinearODEAnalysis { + /// Object type for variables. + pub var_ob_type: QualifiedName, + /// Morphism type for positive links. + pub pos_link_type: QualifiedPath, + /// Morphism type for negative links. + pub neg_link_type: QualifiedPath, +} + +impl Default for LinearODEAnalysis { + fn default() -> Self { + let ob_type = name("Object"); + Self { + var_ob_type: ob_type.clone(), + pos_link_type: Path::Id(ob_type.clone()), + neg_link_type: Path::single(name("Negative")), + } + } +} + +impl + ODESemanticsAnalysis< + ::ModelType, + ::ParameterType, + > for LinearODEAnalysis +{ + /// Creates a linear system with symbolic rate coefficients. + /// + /// A system of ODEs for building arbitrary LinearODE ODEs from CLDs. + fn build_system_builder( + &self, + model: &::ModelType, + ) -> PolynomialODESystemBuilder<::ParameterType> { + let mut builder = PolynomialODESystemBuilder::new(); + + for var in model.ob_generators_with_type(&self.var_ob_type) { + // For each object, we create a variable. + builder.add_variable(var.clone()); + } + + for mor in model.mor_generators_with_type(&self.pos_link_type) { + let (Some(dom), Some(cod)) = (model.get_dom(&mor), model.get_cod(&mor)) else { + continue; + }; + + // The morphism + // f: x -> y + // becomes the contribution + // \dot{y} += Parameter_f x + builder.add_contribution( + mor.clone(), + cod.clone(), + ContributionSign::Positive, + LinearODEParameter::Parameter { morphism: mor }, + [dom.clone()], + ); + } + + for mor in model.mor_generators_with_type(&self.neg_link_type) { + let (Some(dom), Some(cod)) = (model.get_dom(&mor), model.get_cod(&mor)) else { + continue; + }; + + // The morphism + // f: x -> y + // becomes the contribution + // \dot{y} -= Parameter_f x + builder.add_contribution( + mor.clone(), + cod.clone(), + ContributionSign::Negative, + LinearODEParameter::Parameter { morphism: mor }, + [dom.clone()], + ); + } + + builder + } +} /// Data defining a linear ODE problem for a model. #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] @@ -45,73 +157,34 @@ 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 + Zero, +impl ODESemanticsProblemData<::ParameterType> + for LinearODEProblemData { - let system = 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(), - }; - system.normalize() -} + fn initial_values(&self) -> HashMap { + self.initial_values.clone() + } -impl SignedCoefficientBuilder { - /// Linear ODE analysis for a model of a double theory. - /// - /// This analysis is a special case of linear ODE analysis for *extended* causal - /// loop diagrams but can serve as a simple/naive semantics for causal loop - /// diagrams, hopefully useful for toy models and demonstration purposes. - pub fn linear_ode_analysis( - &self, - model: &DiscreteDblModel, - data: LinearODEProblemData, - ) -> ODEAnalysis> { - let (system, ob_index) = self.linear_ode_system(model); - let n = ob_index.len(); - - 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 = system - .extend_scalars(|poly| { - poly.eval(|id| data.coefficients.get(id).copied().unwrap_or_default()) - }) - .to_numerical(); - let problem = ODEProblem::new(system, x0).end_time(data.duration); - ODEAnalysis::new(problem, ob_index) + fn duration(&self) -> f32 { + self.duration } - /// Linear ODE system for a model of a double theory. - pub fn linear_ode_system( + fn extend_scalars( &self, - model: &DiscreteDblModel, - ) -> ( - PolynomialSystem, u8>, - IndexMap, - ) { - let (matrix, ob_index) = self.build_matrix(model); - let system = linear_polynomial_system(&ob_index.keys().cloned().collect_vec(), matrix); - (system, ob_index) + sys: PolynomialSystem< + QualifiedName, + Parameter<::ParameterType>, + i8, + >, + ) -> PolynomialSystem { + let sys = sys.extend_scalars(|poly| { + poly.eval(|param| match param { + LinearODEParameter::Parameter { morphism } => { + self.coefficients.get(morphism).cloned().unwrap_or_default() + } + }) + }); + + sys.normalize() } } @@ -121,43 +194,89 @@ mod test { use std::rc::Rc; use super::*; - use crate::stdlib; - use crate::{one::Path, zero::name}; + use crate::{ + dbl::model::MutDblModel, + simulate::ode::LatexEquation, + stdlib::{models::*, theories::*}, + }; + + // Symbolic tests. - fn builder() -> SignedCoefficientBuilder { - SignedCoefficientBuilder::new(name("Object")) - .add_positive(Path::Id(name("Object"))) - .add_negative(Path::single(name("Negative"))) + #[test] + fn predator_prey_symbolic() { + let th = Rc::new(th_signed_category()); + let model = negative_feedback(th); + let sys = LinearODEAnalysis::default().build_system(&model); + let expected = expect!([r#" + dx = -Parameter(negative) y + dy = Parameter(positive) x + "#]); + expected.assert_eq(&sys.to_string()); } #[test] - fn negative_feedback_symbolic() { - let th = Rc::new(stdlib::theories::th_signed_category()); - let neg_feedback = stdlib::models::negative_feedback(th); - let (sys, _) = builder().linear_ode_system(&neg_feedback); - let expected = expect![[r#" - dx = -negative y - dy = positive x - "#]]; + fn complicated_symbolic() { + let th = Rc::new(th_signed_category()); + let mut model = DiscreteDblModel::new(th); + model.add_ob(name("a"), name("Object")); + model.add_ob(name("b"), name("Object")); + model.add_ob(name("c"), name("Object")); + model.add_ob(name("d"), name("Object")); + model.add_mor(name("f"), name("a"), name("b"), Path::Id(name("Object"))); + model.add_mor(name("g"), name("b"), name("a"), Path::Id(name("Object"))); + model.add_mor(name("h"), name("b"), name("a"), name("Negative").into()); + model.add_mor(name("i"), name("a"), name("c"), name("Negative").into()); + model.add_mor(name("j"), name("c"), name("d"), Path::Id(name("Object"))); + model.add_mor(name("k"), name("d"), name("b"), name("Negative").into()); + let sys = LinearODEAnalysis::default().build_system(&model); + let expected = expect!([r#" + da = (Parameter(g) - Parameter(h)) b + db = Parameter(f) a - Parameter(k) d + dc = -Parameter(i) a + dd = Parameter(j) c + "#]); expected.assert_eq(&sys.to_string()); } + // Test for LaTeX. + #[test] - fn negative_feedback_numerical() { - let th = Rc::new(stdlib::theories::th_signed_category()); - let neg_feedback = stdlib::models::negative_feedback(th); + fn to_latex() { + let th = Rc::new(th_signed_category()); + let model = negative_feedback(th); + let sys = LinearODEAnalysis::default().build_system(&model); + let expected = vec![ + LatexEquation { + lhs: "\\frac{\\mathrm{d}}{\\mathrm{d}t} x".to_string(), + rhs: "-Parameter(negative) \\cdot y".to_string(), + }, + LatexEquation { + lhs: "\\frac{\\mathrm{d}}{\\mathrm{d}t} y".to_string(), + rhs: "Parameter(positive) \\cdot x".to_string(), + }, + ]; + assert_eq!(expected, sys.to_latex_equations()); + } + + // Numerical test. + + #[test] + fn predator_prey_numerical() { + let th = Rc::new(th_signed_category()); + let model = negative_feedback(th); let data = LinearODEProblemData { - coefficients: [(name("positive"), 2.0), (name("negative"), 1.0)].into_iter().collect(), + coefficients: [(name("positive"), 3.0), (name("negative"), 2.0)].into_iter().collect(), initial_values: [(name("x"), 1.0), (name("y"), 1.0)].into_iter().collect(), duration: 10.0, }; - let sys = builder().linear_ode_analysis(&neg_feedback, data).problem.system; - let expected = expect![[r#" - dx0 = -x1 - dx1 = 2 x0 - "#]]; - expected.assert_eq(&sys.to_string()); + let sys = LinearODEAnalysis::default().build_system(&model); + let analysis = data.extend_scalars(sys); + let expected = expect!([r#" + dx = -2 y + dy = 3 x + "#]); + expected.assert_eq(&analysis.to_string()); } } diff --git a/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs b/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs index 9a1853f16..f335ca12d 100644 --- a/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs +++ b/packages/catlog/src/stdlib/analyses/ode/lotka_volterra.rs @@ -1,29 +1,165 @@ //! Lotka-Volterra ODE analysis of models. //! -//! The main entry point for this module is -//! [`lotka_volterra_analysis`](SignedCoefficientBuilder::lotka_volterra_analysis). +//! This follows the structure of [`ode::ode_semantics`], implementing `ODESemantics` for +//! the struct `LotkaVolterraSemantics`. +//! +//! [`ode::ode_semantics`]: crate::stdlib::analyses::ode::ode_semantics use std::collections::HashMap; -use std::hash::Hash; -use std::ops::Add; - -use indexmap::IndexMap; -use itertools::Itertools; -use nalgebra::{DMatrix, DVector, Scalar}; -use num_traits::{One, Zero}; +use std::fmt; #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; #[cfg(feature = "serde-wasm")] use tsify::Tsify; -use super::{ODEAnalysis, Parameter, SignedCoefficientBuilder}; -use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem}; -use crate::{ - dbl::model::DiscreteDblModel, - one::QualifiedPath, - zero::{QualifiedName, alg::Polynomial, rig::Monomial}, +use super::Parameter; +use crate::dbl::model::{FpDblModel, MutDblModel}; +use crate::one::Path; +use crate::simulate::ode::PolynomialSystem; +use crate::stdlib::analyses::ode::ode_semantics::{ + ContributionSign, ODEParameterType, ODESemantics, ODESemanticsAnalysis, + ODESemanticsProblemData, PolynomialODESystemBuilder, }; +use crate::zero::name; +use crate::{dbl::model::DiscreteDblModel, one::QualifiedPath, zero::QualifiedName}; + +/// Implementing Lotka-Volterra as an ODE semantics for models of type `DiscreteDblModel`. +pub struct LotkaVolterraSemantics; + +impl ODESemantics for LotkaVolterraSemantics { + type ModelType = DiscreteDblModel; + type ParameterType = LotkaVolterraParameter; + type AnalysisType = LotkaVolterraAnalysis; + type ProblemDataType = LotkaVolterraProblemData; +} + +/// Parameters in the Lotka-Volterra equations come in two flavours, corresponding to +/// either variables or links. +#[derive(PartialEq, Eq, PartialOrd, Ord, Clone)] +pub enum LotkaVolterraParameter { + /// The parameter associated to a variable. + Growth { + /// The variable. + variable: QualifiedName, + }, + /// The parameter associated to a link. + Interaction { + /// The link. + link: QualifiedName, + }, +} + +impl fmt::Display for LotkaVolterraParameter { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match &self { + Self::Growth { variable } => { + write!(f, "Growth({})", variable) + } + Self::Interaction { link } => { + write!(f, "Interaction({})", link) + } + } + } +} + +impl ODEParameterType for LotkaVolterraParameter {} + +/// This Lotka-Volterra ODE analysis is intended for application to CLDs. +pub struct LotkaVolterraAnalysis { + /// Object type for variables. + pub var_ob_type: QualifiedName, + /// Morphism type for positive links. + pub pos_link_type: QualifiedPath, + /// Morphism type for negative links. + pub neg_link_type: QualifiedPath, +} + +impl Default for LotkaVolterraAnalysis { + fn default() -> Self { + let ob_type = name("Object"); + Self { + var_ob_type: ob_type.clone(), + pos_link_type: Path::Id(ob_type.clone()), + neg_link_type: Path::single(name("Negative")), + } + } +} + +impl + ODESemanticsAnalysis< + ::ModelType, + ::ParameterType, + > for LotkaVolterraAnalysis +{ + /// Creates a Lotka-Volterra system with symbolic rate coefficients. + /// + /// 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) + /// and [our paper on regulatory networks](crate::refs::RegNets). + fn build_system_builder( + &self, + model: &::ModelType, + ) -> PolynomialODESystemBuilder<::ParameterType> { + let mut builder = PolynomialODESystemBuilder::new(); + + for var in model.ob_generators_with_type(&self.var_ob_type) { + // For each object, we create a variable. + builder.add_variable(var.clone()); + + // The object + // x + // becomes the contribution + // \dot{x} += Growth_x \cdot x + builder.add_contribution( + var.clone(), + var.clone(), + ContributionSign::Positive, + LotkaVolterraParameter::Growth { variable: var.clone() }, + [var], + ); + } + + for mor in model.mor_generators_with_type(&self.pos_link_type) { + let (Some(dom), Some(cod)) = (model.get_dom(&mor), model.get_cod(&mor)) else { + continue; + }; + + // The morphism + // f: x -> y + // becomes the contribution + // \dot{y} += Interaction_f \cdot xy + builder.add_contribution( + mor.clone(), + cod.clone(), + ContributionSign::Positive, + LotkaVolterraParameter::Interaction { link: mor }, + [dom.clone(), cod.clone()], + ); + } + + for mor in model.mor_generators_with_type(&self.neg_link_type) { + let (Some(dom), Some(cod)) = (model.get_dom(&mor), model.get_cod(&mor)) else { + continue; + }; + + // The morphism + // f: x -> y + // becomes the contribution + // \dot{y} -= Interaction_f \cdot xy + builder.add_contribution( + mor.clone(), + cod.clone(), + ContributionSign::Negative, + LotkaVolterraParameter::Interaction { link: mor }, + [dom.clone(), cod.clone()], + ); + } + + builder + } +} /// Data defining a Lotka-Volterra ODE problem for a model. #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] @@ -49,94 +185,38 @@ 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 + Zero, +impl ODESemanticsProblemData<::ParameterType> + for LotkaVolterraProblemData { - let system = 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(), - }; - system.normalize() -} + fn initial_values(&self) -> HashMap { + self.initial_values.clone() + } -impl SignedCoefficientBuilder { - /// Lotka-Volterra ODE analysis for a model of a double theory. - /// - /// The main application we have in mind is the Lotka-Volterra ODE semantics for - /// signed graphs described in our [paper on regulatory - /// networks](crate::refs::RegNets). - pub fn lotka_volterra_analysis( - &self, - model: &DiscreteDblModel, - data: LotkaVolterraProblemData, - ) -> ODEAnalysis> { - let (system, ob_index) = self.lotka_volterra_system(model); - let n = ob_index.len(); - - 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 = system - .extend_scalars(|poly| { - poly.eval(|id| { - data.interaction_coeffs - .get(id) - .or(data.growth_rates.get(id)) - .copied() - .unwrap_or_default() - }) - }) - .to_numerical(); - let problem = ODEProblem::new(system, x0).end_time(data.duration); - ODEAnalysis::new(problem, ob_index) + fn duration(&self) -> f32 { + self.duration } - /// Lotka-Volterra ODE system for an model of a double theory. - pub fn lotka_volterra_system( + fn extend_scalars( &self, - model: &DiscreteDblModel, - ) -> ( - PolynomialSystem, u8>, - IndexMap, - ) { - let (matrix, ob_index) = self.build_matrix(model); - let n = ob_index.len(); - - 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 system = lotka_volterra_system(&ob_index.keys().cloned().collect_vec(), matrix, b); - (system, ob_index) + sys: PolynomialSystem< + QualifiedName, + Parameter<::ParameterType>, + i8, + >, + ) -> PolynomialSystem { + let sys = sys.extend_scalars(|poly| { + poly.eval(|param| match param { + LotkaVolterraParameter::Growth { variable } => { + // FIXME: this won't work, because `variable` will now be `Growth.variable` + self.growth_rates.get(variable).cloned().unwrap_or_default() + } + LotkaVolterraParameter::Interaction { link } => { + self.interaction_coeffs.get(link).cloned().unwrap_or_default() + } + }) + }); + + sys.normalize() } } @@ -146,32 +226,76 @@ mod test { use std::rc::Rc; use super::*; - use crate::stdlib; - use crate::{one::Path, zero::name}; + use crate::{ + dbl::model::MutDblModel, + simulate::ode::LatexEquation, + stdlib::{models::*, theories::*}, + }; - fn builder() -> SignedCoefficientBuilder { - SignedCoefficientBuilder::new(name("Object")) - .add_positive(Path::Id(name("Object"))) - .add_negative(Path::single(name("Negative"))) - } + // Symbolic tests. #[test] fn predator_prey_symbolic() { - let th = Rc::new(stdlib::theories::th_signed_category()); - let neg_feedback = stdlib::models::negative_feedback(th); - let (sys, _) = builder().lotka_volterra_system(&neg_feedback); - let sys = sys.extend_scalars(|coef| coef.map_variables(|name| format!("Param({name})"))); + let th = Rc::new(th_signed_category()); + let model = negative_feedback(th); + let sys = LotkaVolterraAnalysis::default().build_system(&model); let expected = expect!([r#" - dx = Param(x) x - Param(negative) x y - dy = Param(positive) x y + Param(y) y + dx = Growth(x) x - Interaction(negative) x y + dy = Interaction(positive) x y + Growth(y) y "#]); expected.assert_eq(&sys.to_string()); } + #[test] + fn complicated_symbolic() { + let th = Rc::new(th_signed_category()); + let mut model = DiscreteDblModel::new(th); + model.add_ob(name("a"), name("Object")); + model.add_ob(name("b"), name("Object")); + model.add_ob(name("c"), name("Object")); + model.add_ob(name("d"), name("Object")); + model.add_mor(name("f"), name("a"), name("b"), Path::Id(name("Object"))); + model.add_mor(name("g"), name("b"), name("a"), Path::Id(name("Object"))); + model.add_mor(name("h"), name("b"), name("a"), name("Negative").into()); + model.add_mor(name("i"), name("a"), name("c"), name("Negative").into()); + model.add_mor(name("j"), name("c"), name("d"), Path::Id(name("Object"))); + model.add_mor(name("k"), name("d"), name("b"), name("Negative").into()); + let sys = LotkaVolterraAnalysis::default().build_system(&model); + let expected = expect!([r#" + da = Growth(a) a + (Interaction(g) - Interaction(h)) a b + db = Interaction(f) a b + Growth(b) b - Interaction(k) b d + dc = -Interaction(i) a c + Growth(c) c + dd = Interaction(j) c d + Growth(d) d + "#]); + expected.assert_eq(&sys.to_string()); + } + + // Test for LaTeX. + + #[test] + fn to_latex() { + let th = Rc::new(th_signed_category()); + let model = negative_feedback(th); + let sys = LotkaVolterraAnalysis::default().build_system(&model); + let expected = vec![ + LatexEquation { + lhs: "\\frac{\\mathrm{d}}{\\mathrm{d}t} x".to_string(), + rhs: "Growth(x) \\cdot x - Interaction(negative) \\cdot x \\cdot y".to_string(), + }, + LatexEquation { + lhs: "\\frac{\\mathrm{d}}{\\mathrm{d}t} y".to_string(), + rhs: "Interaction(positive) \\cdot x \\cdot y + Growth(y) \\cdot y".to_string(), + }, + ]; + assert_eq!(expected, sys.to_latex_equations()); + } + + // Numerical test. + #[test] fn predator_prey_numerical() { - let th = Rc::new(stdlib::theories::th_signed_category()); - let neg_feedback = stdlib::models::negative_feedback(th); + let th = Rc::new(th_signed_category()); + let model = negative_feedback(th); let data = LotkaVolterraProblemData { interaction_coeffs: [(name("positive"), 1.0), (name("negative"), 1.0)] @@ -182,11 +306,12 @@ mod test { duration: 10.0, }; - let sys = builder().lotka_volterra_analysis(&neg_feedback, data).problem.system; + let sys = LotkaVolterraAnalysis::default().build_system(&model); + let analysis = data.extend_scalars(sys); let expected = expect!([r#" - dx0 = 2 x0 - x0 x1 - dx1 = x0 x1 - x1 + dx = 2 x - x y + dy = x y - y "#]); - expected.assert_eq(&sys.to_string()); + expected.assert_eq(&analysis.to_string()); } } diff --git a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs index 9daac606b..b70761dc0 100644 --- a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs +++ b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs @@ -7,24 +7,43 @@ use std::{collections::HashMap, fmt}; -use indexmap::IndexMap; -use nalgebra::DVector; -use num_traits::Zero; - #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; #[cfg(feature = "serde-wasm")] use tsify::Tsify; -use super::{ODEAnalysis, Parameter}; -use crate::dbl::{ - model::{DiscreteTabModel, FpDblModel, ModalDblModel, TabEdge}, - theory::{ModalMorType, ModalObType, TabMorType, TabObType, Unital}, -}; -use crate::one::FgCategory; -use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem}; +use super::Parameter; +use crate::simulate::ode::PolynomialSystem; +use crate::stdlib::analyses::ode::ode_semantics::*; use crate::stdlib::analyses::petri::transition_interface; -use crate::zero::{QualifiedName, alg::Polynomial, name, rig::Monomial}; +use crate::stdlib::analyses::stock_flow::flow_interface; +use crate::zero::{QualifiedName, name}; +use crate::{ + dbl::{ + model::{DiscreteTabModel, FpDblModel, ModalDblModel}, + theory::{ModalMorType, ModalObType, TabMorType, TabObType, Unital}, + }, + zero::name_seg, +}; + +/// Mass-action semantics for Petri nets. +pub struct PetriNetMassActionSemantics; +/// Mass-action semantics for stock-flow diagrams. +pub struct StockFlowMassActionSemantics; + +impl ODESemantics for PetriNetMassActionSemantics { + type ModelType = ModalDblModel; + type ParameterType = MassActionParameter; + type AnalysisType = PetriNetMassActionAnalysis; + type ProblemDataType = MassActionProblemData; +} + +impl ODESemantics for StockFlowMassActionSemantics { + type ModelType = DiscreteTabModel; + type ParameterType = MassActionParameter; + type AnalysisType = StockFlowMassActionAnalysis; + type ProblemDataType = MassActionProblemData; +} /// There are three types of mass-action semantics, each more expressive than the previous: /// - balanced @@ -49,22 +68,23 @@ pub enum MassConservationType { #[cfg_attr(feature = "serde-wasm", derive(Tsify))] #[cfg_attr(feature = "serde-wasm", tsify(into_wasm_abi, from_wasm_abi))] pub enum RateGranularity { - /// Each transition gets assigned a single consumption and single production rate. + /// Each flow gets assigned a single consumption and single production rate. PerTransition, - /// Each transition gets assigned a consumption rate for each input place and - /// a production rate for each output place. + /// Each flow gets assigned a consumption rate for each input stock and + /// a production rate for each output stock. PerPlace, } -/// Parameters in the generated polynomial equations are *undirected* in the -/// balanced case and *directed* in the unbalanced case. +/// Now, corresponding to each term of `MassConvervationType`, we have different +/// terms for `MassActionParameter`. Parameters in the generated polynomial equations +/// are *undirected* in the balanced case and *directed* in the unbalanced case. #[derive(PartialEq, Eq, PartialOrd, Ord, Clone)] -pub enum FlowParameter { +pub enum MassActionParameter { /// If mass is conserved, we don't need to worry whether a flow is incoming or outgoing. Balanced { /// Since there is no direction, the rate parameter corresponds to a single transition. - transition: QualifiedName, + flow: QualifiedName, }, /// If mass is not conserved, then we need to know whether a flow is incoming or outgoing. Unbalanced { @@ -78,19 +98,19 @@ pub enum FlowParameter { /// Depending on the rate granularity, the parameters are specified by different structures. #[derive(PartialEq, Eq, PartialOrd, Ord, Clone)] pub enum RateParameter { - /// For per transition rates, we simply need to know the associated transition. + /// For per flow rates, we simply need to know the associated flow. PerTransition { - /// The transition to which we associate the rate parameter. - transition: QualifiedName, + /// The flow to which we associate the rate parameter. + flow: QualifiedName, }, - /// For per place rates, we need to know both the transition and the corresponding - /// input/output place. + /// For per stock rates, we need to know both the transition and the corresponding + /// input/output stock. PerPlace { - /// The transition whose input/output objects we wish to associate rate parameters. - transition: QualifiedName, - /// The input/output object to which we associate the rate parameter. - place: QualifiedName, + /// The flow whose input/output objects we wish to associate rate parameters. + flow: QualifiedName, + /// The input/output stock to which we associate the rate parameter. + stock: QualifiedName, }, } @@ -106,33 +126,33 @@ pub enum Direction { OutgoingFlow, } -impl fmt::Display for FlowParameter { +impl fmt::Display for MassActionParameter { fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { match &self { - FlowParameter::Balanced { transition: trans } => { + Self::Balanced { flow: trans } => { write!(f, "{}", trans) } - FlowParameter::Unbalanced { + Self::Unbalanced { direction: Direction::IncomingFlow, - parameter: RateParameter::PerTransition { transition: trans }, + parameter: RateParameter::PerTransition { flow: trans }, } => { write!(f, "Incoming({})", trans) } - FlowParameter::Unbalanced { + Self::Unbalanced { direction: Direction::IncomingFlow, - parameter: RateParameter::PerPlace { transition: trans, place: output }, + parameter: RateParameter::PerPlace { flow: trans, stock: output }, } => { write!(f, "([{}]->{})", trans, output) } - FlowParameter::Unbalanced { + Self::Unbalanced { direction: Direction::OutgoingFlow, - parameter: RateParameter::PerTransition { transition: trans }, + parameter: RateParameter::PerTransition { flow: trans }, } => { write!(f, "Outgoing({})", trans) } - FlowParameter::Unbalanced { + Self::Unbalanced { direction: Direction::OutgoingFlow, - parameter: RateParameter::PerPlace { transition: trans, place: input }, + parameter: RateParameter::PerPlace { flow: trans, stock: input }, } => { write!(f, "({}->[{}])", input, trans) } @@ -140,51 +160,7 @@ impl fmt::Display for FlowParameter { } } -/// Data defining an unbalanced mass-action ODE problem for a model. -#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] -#[cfg_attr(feature = "serde-wasm", derive(Tsify))] -#[cfg_attr( - feature = "serde-wasm", - tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object) -)] -pub struct MassActionProblemData { - /// Whether or not mass is conserved. - #[cfg_attr(feature = "serde", serde(rename = "massConservationType"))] - pub mass_conservation_type: MassConservationType, - - /// Map from morphism IDs to consumption rate coefficients (nonnegative reals), - /// for the balanced per transition case. - /// N.B. This is renamed to "rates" in catlog-wasm for backwards compatibility. - #[cfg_attr(feature = "serde", serde(rename = "rates"))] - transition_rates: HashMap, - - /// Map from morphism IDs to consumption rate coefficients (nonnegative reals), - /// for the unbalanced per transition case. - #[cfg_attr(feature = "serde", serde(rename = "transitionConsumptionRates"))] - transition_consumption_rates: HashMap, - - /// Map from morphism IDs to production rate coefficients (nonnegative reals), - /// for the unbalanced per transition case. - #[cfg_attr(feature = "serde", serde(rename = "transitionProductionRates"))] - transition_production_rates: HashMap, - - /// Map from morphism IDs to (map from input objects to consumption rate coefficients), - /// for the unbalanced per place case (nonnegative reals). - #[cfg_attr(feature = "serde", serde(rename = "placeConsumptionRates"))] - place_consumption_rates: HashMap>, - - /// Map from morphism IDs to (map from output objects to production rate coefficients), - /// for the unbalanced per place case (nonnegative reals). - #[cfg_attr(feature = "serde", serde(rename = "placeProductionRates"))] - place_production_rates: HashMap>, - - /// Map from object IDs to initial values (nonnegative reals). - #[cfg_attr(feature = "serde", serde(rename = "initialValues"))] - pub initial_values: HashMap, - - /// Duration of simulation. - pub duration: f32, -} +impl ODEParameterType for MassActionParameter {} /// Mass-action ODE analysis for Petri nets. /// @@ -195,6 +171,8 @@ pub struct PetriNetMassActionAnalysis { pub place_ob_type: ModalObType, /// Morphism type for transitions. pub transition_mor_type: ModalMorType, + /// Mass-conservation type. + pub mass_conservation_type: MassConservationType, } impl Default for PetriNetMassActionAnalysis { @@ -203,104 +181,118 @@ impl Default for PetriNetMassActionAnalysis { Self { place_ob_type: ob_type.clone(), transition_mor_type: ModalMorType::Zero(ob_type), + mass_conservation_type: MassConservationType::Balanced, } } } -impl PetriNetMassActionAnalysis { - /// Creates a mass-action system with symbolic rate coefficients. - pub fn build_system( +impl + ODESemanticsAnalysis< + ::ModelType, + ::ParameterType, + > for PetriNetMassActionAnalysis +{ + fn build_system_builder( &self, - model: &ModalDblModel, - mass_conservation_type: MassConservationType, - ) -> PolynomialSystem, i8> { - let mut sys = PolynomialSystem::new(); - for ob in model.ob_generators_with_type(&self.place_ob_type) { - sys.add_term(ob, Polynomial::zero()); + model: &::ModelType, + ) -> PolynomialODESystemBuilder<::ParameterType> + { + let mut builder = PolynomialODESystemBuilder::new(); + + for place in model.ob_generators_with_type(&self.place_ob_type) { + // For each place, we create a variable. + builder.add_variable(place.clone()); } - for mor in model.mor_generators_with_type(&self.transition_mor_type) { - let (inputs, outputs) = transition_interface(model, &mor); - let term: Monomial<_, _> = - inputs.iter().map(|ob| (ob.clone().unwrap_generator(), 1)).collect(); - match mass_conservation_type { - MassConservationType::Balanced => { - let term: Polynomial<_, _, _> = [( - Parameter::generator(FlowParameter::Balanced { transition: mor }), - term.clone(), - )] - .into_iter() - .collect(); - - for input in inputs { - sys.add_term(input.unwrap_generator(), -term.clone()); + for transition in model.mor_generators_with_type(&self.transition_mor_type) { + let interface = transition_interface(model, &transition); + let (inputs, outputs) = + (interface.input_places.clone(), interface.output_places.clone()); + + // Each transition gives a positive contribution to each term corresponding to + // one of its outputs, and a negative contribution to each term corresponding to + // one of its inputs. For example, a single transition T: [a,b] -> [x,y] will give + // four contributions, namely two positive contributions (ab -> x , ab -> y) + // and two negative (ab -> a , ab -> b). + + for output in outputs.clone() { + let id = output.cons(name_seg("ToOutput")).cons(transition.only().unwrap()); + // The transition + // T : [x_1, ..., x_n] -> [y_1, ..., y_n] + // becomes the contributions + // \dot{y_i} += Parameter_! \cdot x_1...x_n + // where Parameter_! depends on `mass_conservation_type`: + // Balanced => Parameter_T + // Unbalanced::PerTransition => Parameter_T^inflow + // Unbalanced::PerPlace => Parameter_{T,y_i}^inflow + let parameter = match self.mass_conservation_type { + MassConservationType::Balanced => { + MassActionParameter::Balanced { flow: transition.clone() } } + MassConservationType::Unbalanced(granularity) => match granularity { + RateGranularity::PerTransition => MassActionParameter::Unbalanced { + direction: Direction::IncomingFlow, + parameter: RateParameter::PerTransition { flow: transition.clone() }, + }, + RateGranularity::PerPlace => MassActionParameter::Unbalanced { + direction: Direction::IncomingFlow, + parameter: RateParameter::PerPlace { + flow: transition.clone(), + stock: output.clone(), + }, + }, + }, + }; + + builder.add_contribution( + id, + output, + ContributionSign::Positive, + parameter, + inputs.clone(), + ); + } - for output in outputs { - sys.add_term(output.unwrap_generator(), term.clone()); - } - } - - MassConservationType::Unbalanced(granularity) => { - for input in inputs { - let input_term: Polynomial<_, _, _> = match granularity { - RateGranularity::PerTransition => [( - Parameter::generator(FlowParameter::Unbalanced { - direction: Direction::OutgoingFlow, - parameter: RateParameter::PerTransition { - transition: mor.clone(), - }, - }), - term.clone(), - )], - RateGranularity::PerPlace => [( - Parameter::generator(FlowParameter::Unbalanced { - direction: Direction::OutgoingFlow, - parameter: RateParameter::PerPlace { - transition: mor.clone(), - place: input.clone().unwrap_generator(), - }, - }), - term.clone(), - )], - } - .into_iter() - .collect(); - - sys.add_term(input.unwrap_generator(), -input_term.clone()); + for input in inputs.clone() { + let id = input.cons(name_seg("ToInput")).cons(transition.only().unwrap()); + // The transition + // T : [x_1, ..., x_n] -> [y_1, ..., y_n] + // becomes the contributions + // \dot{x_i} -= Parameter_! \cdot x_1...x_n + // where Parameter_! depends on `mass_conservation_type`: + // Balanced => Parameter_T + // Unbalanced::PerTransition => Parameter_T^outflow + // Unbalanced::PerPlace => Parameter_{T,x_i}^outflow + let parameter = match self.mass_conservation_type { + MassConservationType::Balanced => { + MassActionParameter::Balanced { flow: transition.clone() } } - for output in outputs { - let output_term: Polynomial<_, _, _> = match granularity { - RateGranularity::PerTransition => [( - Parameter::generator(FlowParameter::Unbalanced { - direction: Direction::IncomingFlow, - parameter: RateParameter::PerTransition { - transition: mor.clone(), - }, - }), - term.clone(), - )], - RateGranularity::PerPlace => [( - Parameter::generator(FlowParameter::Unbalanced { - direction: Direction::IncomingFlow, - parameter: RateParameter::PerPlace { - transition: mor.clone(), - place: output.clone().unwrap_generator(), - }, - }), - term.clone(), - )], - } - .into_iter() - .collect(); - - sys.add_term(output.unwrap_generator(), output_term.clone()); - } - } + MassConservationType::Unbalanced(granularity) => match granularity { + RateGranularity::PerTransition => MassActionParameter::Unbalanced { + direction: Direction::OutgoingFlow, + parameter: RateParameter::PerTransition { flow: transition.clone() }, + }, + RateGranularity::PerPlace => MassActionParameter::Unbalanced { + direction: Direction::OutgoingFlow, + parameter: RateParameter::PerPlace { + flow: transition.clone(), + stock: input.clone(), + }, + }, + }, + }; + + builder.add_contribution( + id, + input, + ContributionSign::Negative, + parameter, + inputs.clone(), + ); } } - sys.normalize() + builder } } @@ -314,156 +306,211 @@ pub struct StockFlowMassActionAnalysis { pub pos_link_mor_type: TabMorType, /// Morphism type for negative links from stocks to flows. pub neg_link_mor_type: TabMorType, + /// Mass-conservation type. + pub mass_conservation_type: MassConservationType, } impl Default for StockFlowMassActionAnalysis { fn default() -> Self { - let stock_ob_type = TabObType::Basic(name("Object")); - let flow_mor_type = TabMorType::Hom(Box::new(stock_ob_type.clone())); + let ob_type = TabObType::Basic(name("Object")); Self { - stock_ob_type, - flow_mor_type, + stock_ob_type: ob_type.clone(), + flow_mor_type: TabMorType::Hom(Box::new(ob_type.clone())), pos_link_mor_type: TabMorType::Basic(name("Link")), neg_link_mor_type: TabMorType::Basic(name("NegativeLink")), + mass_conservation_type: MassConservationType::Balanced, } } } -impl StockFlowMassActionAnalysis { - /// Creates a mass-action system with symbolic rate coefficients. - pub fn build_system( +impl + ODESemanticsAnalysis< + ::ModelType, + ::ParameterType, + > for StockFlowMassActionAnalysis +{ + fn build_system_builder( &self, - model: &DiscreteTabModel, - mass_conservation_type: MassConservationType, - ) -> PolynomialSystem, i8> { - let terms: Vec<_> = self.flow_monomials(model).into_iter().collect(); - - let mut sys = PolynomialSystem::new(); - for ob in model.ob_generators_with_type(&self.stock_ob_type) { - sys.add_term(ob, Polynomial::zero()); + model: &::ModelType, + ) -> PolynomialODESystemBuilder<::ParameterType> + { + let mut builder = PolynomialODESystemBuilder::new(); + + for stock in model.ob_generators_with_type(&self.stock_ob_type) { + // For each stock, we create a variable. + builder.add_variable(stock.clone()); } - for (flow, term) in terms { - let dom = model.mor_generator_dom(&flow).unwrap_basic(); - let cod = model.mor_generator_cod(&flow).unwrap_basic(); - match mass_conservation_type { + + for flow in model.mor_generators_with_type(&self.flow_mor_type) { + let interface = flow_interface(model, &flow); + let (input, output) = (interface.input_stock, interface.output_stock); + + // Each flow gives a positive contribution to the term corresponding to its output, and + // a negative contribution to the term corresponding to its input; the term is given by + // the product of the input with the sources of all incoming links. + let monomial = [interface.input_pos_link_doms, vec![input.clone()]].concat(); + + // The flow + // F : a -> b + // with links + // l_i : x_i -> F + // becomes the contributions + // \dot{b} += Parameter_! \cdot a x_1.. x_n + // \dot{a} -= Parameter_? \cdot a x_1.. x_n + // where Parameter_! and Parameter_? depend on `mass_conservation_type`: + // Balanced => Parameter_! = Parameter_F + // Parameter_? = Parameter_F + // Unbalanced::PerTransition => Parameter_! = Parameter_F^inflow + // Parameter_? = Parameter_F^outflow + + let output_id = output.cons(name_seg("ToOutput")).cons(flow.only().unwrap()); + let output_parameter = match self.mass_conservation_type { MassConservationType::Balanced => { - let param = Parameter::generator(FlowParameter::Balanced { transition: flow }); - let term: Polynomial<_, _, _> = [(param, term.clone())].into_iter().collect(); - sys.add_term(dom, -term.clone()); - sys.add_term(cod, term); + MassActionParameter::Balanced { flow: flow.clone() } } - MassConservationType::Unbalanced(_) => { - let dom_param = Parameter::generator(FlowParameter::Unbalanced { - direction: Direction::OutgoingFlow, - parameter: RateParameter::PerTransition { transition: flow.clone() }, - }); - let cod_param = Parameter::generator(FlowParameter::Unbalanced { - direction: Direction::IncomingFlow, - parameter: RateParameter::PerTransition { transition: flow }, - }); - let dom_term: Polynomial<_, _, _> = - [(dom_param, term.clone())].into_iter().collect(); - let cod_term: Polynomial<_, _, _> = [(cod_param, term)].into_iter().collect(); - sys.add_term(dom, -dom_term); - sys.add_term(cod, cod_term); + MassConservationType::Unbalanced(_) => MassActionParameter::Unbalanced { + direction: Direction::IncomingFlow, + parameter: RateParameter::PerTransition { flow: flow.clone() }, + }, + }; + builder.add_contribution( + output_id, + output.clone(), + ContributionSign::Positive, + output_parameter, + monomial.clone(), + ); + + let input_id = input.cons(name_seg("ToInput")).cons(flow.only().unwrap()); + let input_parameter = match self.mass_conservation_type { + MassConservationType::Balanced => { + MassActionParameter::Balanced { flow: flow.clone() } } - } + MassConservationType::Unbalanced(_) => MassActionParameter::Unbalanced { + direction: Direction::OutgoingFlow, + parameter: RateParameter::PerTransition { flow: flow.clone() }, + }, + }; + builder.add_contribution( + input_id, + input.clone(), + ContributionSign::Negative, + input_parameter, + monomial, + ); } - sys + + builder } +} - /// Constructs a monomial for each flow in the model. - pub(super) fn flow_monomials( - &self, - model: &DiscreteTabModel, - ) -> HashMap> { - let mut terms: HashMap<_, _> = model - .mor_generators_with_type(&self.flow_mor_type) - .map(|flow| { - let dom = model.mor_generator_dom(&flow).unwrap_basic(); - (flow, Monomial::generator(dom)) - }) - .collect(); +/// Data defining an unbalanced mass-action ODE problem for a model. +#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] +#[cfg_attr(feature = "serde-wasm", derive(Tsify))] +#[cfg_attr( + feature = "serde-wasm", + tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object) +)] +pub struct MassActionProblemData { + /// Whether or not mass is conserved. + #[cfg_attr(feature = "serde", serde(rename = "massConservationType"))] + pub mass_conservation_type: MassConservationType, - let mut multiply_for_link = |link: QualifiedName, exponent: i8| { - let dom = model.mor_generator_dom(&link).unwrap_basic(); - let path = model.mor_generator_cod(&link).unwrap_tabulated(); - let Some(TabEdge::Basic(cod)) = path.only() else { - panic!("Codomain of link should be basic morphism"); - }; - if let Some(term) = terms.get_mut(&cod) { - let mon: Monomial<_, i8> = [(dom, exponent)].into_iter().collect(); - *term = std::mem::take(term) * mon; - } else { - panic!("Codomain of link does not belong to model"); - }; - }; + /// Map from morphism IDs to consumption rate coefficients (nonnegative reals), + /// for the balanced per transition case. + /// N.B. This is renamed to "rates" in catlog-wasm for backwards compatibility. + #[cfg_attr(feature = "serde", serde(rename = "rates"))] + transition_rates: HashMap, - for link in model.mor_generators_with_type(&self.pos_link_mor_type) { - multiply_for_link(link, 1); - } - for link in model.mor_generators_with_type(&self.neg_link_mor_type) { - multiply_for_link(link, -1); - } + /// Map from morphism IDs to consumption rate coefficients (nonnegative reals), + /// for the unbalanced per transition case. + #[cfg_attr(feature = "serde", serde(rename = "transitionConsumptionRates"))] + transition_consumption_rates: HashMap, - terms - } -} + /// Map from morphism IDs to production rate coefficients (nonnegative reals), + /// for the unbalanced per transition case. + #[cfg_attr(feature = "serde", serde(rename = "transitionProductionRates"))] + transition_production_rates: HashMap, -/// Substitutes numerical rate coefficients into a symbolic mass-action system. -pub fn extend_mass_action_scalars( - sys: PolynomialSystem, i8>, - data: &MassActionProblemData, -) -> PolynomialSystem { - let sys = sys.extend_scalars(|poly| { - poly.eval(|flow| match flow { - FlowParameter::Balanced { transition } => { - data.transition_rates.get(transition).cloned().unwrap_or_default() - } - FlowParameter::Unbalanced { direction, parameter } => match (direction, parameter) { - (Direction::IncomingFlow, RateParameter::PerTransition { transition }) => { - data.transition_production_rates.get(transition).cloned().unwrap_or_default() - } - (Direction::OutgoingFlow, RateParameter::PerTransition { transition }) => { - data.transition_consumption_rates.get(transition).cloned().unwrap_or_default() - } - (Direction::IncomingFlow, RateParameter::PerPlace { transition, place }) => data - .place_production_rates - .get(transition) - .and_then(|rate| rate.get(place)) - .copied() - .unwrap_or_default(), - (Direction::OutgoingFlow, RateParameter::PerPlace { transition, place }) => data - .place_consumption_rates - .get(transition) - .and_then(|rate| rate.get(place)) - .copied() - .unwrap_or_default(), - }, - }) - }); + /// Map from morphism IDs to (map from input objects to consumption rate coefficients), + /// for the unbalanced per place case (nonnegative reals). + #[cfg_attr(feature = "serde", serde(rename = "placeConsumptionRates"))] + place_consumption_rates: HashMap>, - sys.normalize() + /// Map from morphism IDs to (map from output objects to production rate coefficients), + /// for the unbalanced per place case (nonnegative reals). + #[cfg_attr(feature = "serde", serde(rename = "placeProductionRates"))] + place_production_rates: HashMap>, + + /// Map from object IDs to initial values (nonnegative reals). + #[cfg_attr(feature = "serde", serde(rename = "initialValues"))] + pub initial_values: HashMap, + + /// Duration of simulation. + pub duration: f32, } -/// Builds the numerical ODE analysis for a mass-action system whose scalars have been substituted. -pub fn into_mass_action_analysis( - sys: PolynomialSystem, - data: MassActionProblemData, -) -> ODEAnalysis> { - let ob_index: IndexMap<_, _> = - sys.components.keys().cloned().enumerate().map(|(i, x)| (x, i)).collect(); - let n = ob_index.len(); +impl ODESemanticsProblemData for MassActionProblemData { + fn initial_values(&self) -> HashMap { + self.initial_values.clone() + } - 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); + fn duration(&self) -> f32 { + self.duration + } - let num_sys = sys.to_numerical(); - let problem = ODEProblem::new(num_sys, x0).end_time(data.duration); + fn extend_scalars( + &self, + sys: PolynomialSystem, i8>, + ) -> PolynomialSystem { + let sys = sys.extend_scalars(|poly| { + poly.eval(|flow| match flow { + MassActionParameter::Balanced { flow: transition } => { + self.transition_rates.get(transition).cloned().unwrap_or_default() + } + MassActionParameter::Unbalanced { direction, parameter } => { + match (direction, parameter) { + ( + Direction::IncomingFlow, + RateParameter::PerTransition { flow: transition }, + ) => self + .transition_production_rates + .get(transition) + .cloned() + .unwrap_or_default(), + ( + Direction::OutgoingFlow, + RateParameter::PerTransition { flow: transition }, + ) => self + .transition_consumption_rates + .get(transition) + .cloned() + .unwrap_or_default(), + ( + Direction::IncomingFlow, + RateParameter::PerPlace { flow: transition, stock: place }, + ) => self + .place_production_rates + .get(transition) + .and_then(|rate| rate.get(place)) + .copied() + .unwrap_or_default(), + ( + Direction::OutgoingFlow, + RateParameter::PerPlace { flow: transition, stock: place }, + ) => self + .place_consumption_rates + .get(transition) + .and_then(|rate| rate.get(place)) + .copied() + .unwrap_or_default(), + } + } + }) + }); - ODEAnalysis::new(problem, ob_index) + sys.normalize() + } } #[cfg(test)] @@ -482,8 +529,7 @@ mod tests { fn balanced_stock_flow() { let th = Rc::new(th_category_links()); let model = backward_link(th); - let sys = StockFlowMassActionAnalysis::default() - .build_system(&model, analyses::ode::MassConservationType::Balanced); + let sys = StockFlowMassActionAnalysis::default().build_system(&model); let expected = expect!([r#" dx = -f x y dy = f x y @@ -495,12 +541,13 @@ mod tests { fn unbalanced_stock_flow() { let th = Rc::new(th_category_links()); let model = backward_link(th); - let sys = StockFlowMassActionAnalysis::default().build_system( - &model, - analyses::ode::MassConservationType::Unbalanced( + let sys = StockFlowMassActionAnalysis { + mass_conservation_type: analyses::ode::MassConservationType::Unbalanced( analyses::ode::RateGranularity::PerTransition, ), - ); + ..StockFlowMassActionAnalysis::default() + } + .build_system(&model); let expected = expect!([r#" dx = -Outgoing(f) x y dy = Incoming(f) x y @@ -508,39 +555,6 @@ mod tests { expected.assert_eq(&sys.to_string()); } - // Tests for signed stock-flow diagrams. These all use the negative_backwards_link() - // model, which has a single flow x==f=>y and a single negative link y->f. - - #[test] - fn balanced_signed_stock_flow() { - let th = Rc::new(th_category_signed_links()); - let model = negative_backward_link(th); - let sys = StockFlowMassActionAnalysis::default() - .build_system(&model, analyses::ode::MassConservationType::Balanced); - let expected = expect!([r#" - dx = -f x y^{-1} - dy = f x y^{-1} - "#]); - expected.assert_eq(&sys.to_string()); - } - - #[test] - fn unbalanced_signed_stock_flow() { - let th = Rc::new(th_category_signed_links()); - let model = negative_backward_link(th); - let sys = StockFlowMassActionAnalysis::default().build_system( - &model, - analyses::ode::MassConservationType::Unbalanced( - analyses::ode::RateGranularity::PerTransition, - ), - ); - let expected = expect!([r#" - dx = -Outgoing(f) x y^{-1} - dy = Incoming(f) x y^{-1} - "#]); - expected.assert_eq(&sys.to_string()); - } - // Tests for Petri nets. These all use the catalyzed_reaction() model, which // has a single transition [x,c]-->f-->[y,c]. @@ -548,8 +562,7 @@ mod tests { fn balanced_petri() { let th = Rc::new(th_sym_monoidal_category()); let model = catalyzed_reaction(th); - let sys = PetriNetMassActionAnalysis::default() - .build_system(&model, analyses::ode::MassConservationType::Balanced); + let sys = PetriNetMassActionAnalysis::default().build_system(&model); let expected = expect!([r#" dx = -f c x dy = f c x @@ -562,12 +575,13 @@ mod tests { fn unbalanced_petri_per_transition() { let th = Rc::new(th_sym_monoidal_category()); let model = catalyzed_reaction(th); - let sys = PetriNetMassActionAnalysis::default().build_system( - &model, - analyses::ode::MassConservationType::Unbalanced( + let sys = PetriNetMassActionAnalysis { + mass_conservation_type: analyses::ode::MassConservationType::Unbalanced( analyses::ode::RateGranularity::PerTransition, ), - ); + ..PetriNetMassActionAnalysis::default() + } + .build_system(&model); let expected = expect!([r#" dx = -Outgoing(f) c x dy = Incoming(f) c x @@ -580,12 +594,13 @@ mod tests { fn unbalanced_petri_per_place() { let th = Rc::new(th_sym_monoidal_category()); let model = catalyzed_reaction(th); - let sys = PetriNetMassActionAnalysis::default().build_system( - &model, - analyses::ode::MassConservationType::Unbalanced( + let sys = PetriNetMassActionAnalysis { + mass_conservation_type: analyses::ode::MassConservationType::Unbalanced( analyses::ode::RateGranularity::PerPlace, ), - ); + ..PetriNetMassActionAnalysis::default() + } + .build_system(&model); let expected = expect!([r#" dx = -(x->[f]) c x dy = ([f]->y) c x @@ -600,12 +615,13 @@ mod tests { fn to_latex() { let th = Rc::new(th_category_links()); let model = backward_link(th); - let sys = StockFlowMassActionAnalysis::default().build_system( - &model, - analyses::ode::MassConservationType::Unbalanced( + let sys = StockFlowMassActionAnalysis { + mass_conservation_type: analyses::ode::MassConservationType::Unbalanced( analyses::ode::RateGranularity::PerTransition, ), - ); + ..StockFlowMassActionAnalysis::default() + } + .build_system(&model); let expected = vec![ LatexEquation { lhs: "\\frac{\\mathrm{d}}{\\mathrm{d}t} x".to_string(), diff --git a/packages/catlog/src/stdlib/analyses/ode/mod.rs b/packages/catlog/src/stdlib/analyses/ode/mod.rs index 4c9b2a862..c5ce791d2 100644 --- a/packages/catlog/src/stdlib/analyses/ode/mod.rs +++ b/packages/catlog/src/stdlib/analyses/ode/mod.rs @@ -73,12 +73,12 @@ pub mod kuramoto; pub mod linear_ode; pub mod lotka_volterra; pub mod mass_action; +pub mod ode_semantics; pub mod polynomial_ode; -pub mod signed_coefficients; pub use kuramoto::*; pub use linear_ode::*; pub use lotka_volterra::*; pub use mass_action::*; +pub use ode_semantics::*; pub use polynomial_ode::*; -pub use signed_coefficients::*; diff --git a/packages/catlog/src/stdlib/analyses/ode/ode_semantics.rs b/packages/catlog/src/stdlib/analyses/ode/ode_semantics.rs new file mode 100644 index 000000000..b23a4a67a --- /dev/null +++ b/packages/catlog/src/stdlib/analyses/ode/ode_semantics.rs @@ -0,0 +1,258 @@ +//! Analyses for different ODE semantics on models. +//! +//! Inspired by schema migration, we define the data of an ODE semantics on models in a theory to +//! consist of (in particular) a `PolynomialODESystemBuilder`, which contains all the data needed +//! for [`ode::polynomial_ode::PolynomialODEAnalysis`] to do the following: +//! +//! 1. Build the system as a model of the theory of polynomial ODE systems (i.e. multicategories) +//! with abstract coefficients, using `build_system_custom_parameters()`. +//! 2. Substitute in numerical coefficients, using `extend_polynomial_ode_scalars()`. +//! 3. Build an `ODEAnalysis>` that can be fed into an ODE solver, +//! using `polynomial_ode_analysis()`. +//! +//! In short, this module constructs multicategories from models, and [`ode::polynomial_ode`] then +//! constructs `PolynomialSystem` from multicategories. +//! +//! To implement a new ODE semantics for models in some theory, one essentially needs to create an +//! empty struct and implement `ODESemantics`, and then follow the compiler. For more documentation, +//! see [`ode::polynomial_ode`]; for an example implementation, see [`ode::mass_action`]. +//! +//! [`ode::polynomial_ode`]: crate::stdlib::analyses::ode::polynomial_ode +//! [`ode::polynomial_ode::PolynomialODEAnalysis`]: crate::stdlib::analyses::ode::polynomial_ode::PolynomialODEAnalysis +//! [`ode::mass_action`]: crate::stdlib::analyses::ode::mass_action + +use indexmap::IndexMap; +use nalgebra::DVector; +use std::{collections::HashMap, fmt}; + +use crate::{ + dbl::{ + modal::{List, ModeApp}, + model::{DiscreteDblModel, DiscreteTabModel, ModalDblModel, ModalOb, MutDblModel}, + theory::{NonUnital, Unital}, + }, + one::FgCategory, + simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem}, + stdlib::{ + analyses::ode::{ODEAnalysis, Parameter, PolynomialODEAnalysis}, + th_signed_polynomial_ode_system, + }, + zero::{QualifiedName, name}, +}; + +/// The trait for an ODE semantics on models. +pub trait ODESemantics { + /// The type of the model for which these ODE semantics are intended. + type ModelType: DblModelForODESemantics; + /// The type of the parameters associated to each contribution in the multicategory built from + /// the model. The "default" value for this would be `QualifiedName`, but it can be useful to + /// have a more descriptive type. For example, we might wish for certain parameters to be + /// identified with one another, or to be rendered differently in debug/LaTeX output. For an + /// instructive example, see `MassActionParameter` in `ode::mass_action`. + type ParameterType: ODEParameterType; + /// The data describing the things that the ODE semantics "cares about". (See the documentation + /// for `ODESemanticsAnalysis`). + type AnalysisType: ODESemanticsAnalysis; + /// The data describing how to turn the algebraic system of equations into a simulation, + /// including e.g. which values that appear in the front-end analysis correspond to which + /// parameters within the equations. + type ProblemDataType: ODESemanticsProblemData; +} + +/// The models for which we support ODE semantics need to be sufficiently nice, though +/// these bounds are not particularly restrictive. +pub trait DblModelForODESemantics: + FgCategory + MutDblModel + Clone +{ +} + +impl DblModelForODESemantics for DiscreteDblModel {} +impl DblModelForODESemantics for DiscreteTabModel {} +impl DblModelForODESemantics for ModalDblModel {} +impl DblModelForODESemantics for ModalDblModel {} + +/// The type of the parameters in the ODE system need to be sufficiently nice, though +/// (again) these bounds are not particularly restrictive. +pub trait ODEParameterType: Eq + Ord + Clone + fmt::Display {} + +/// The simplest type for parameters is `QualifiedName`. +impl ODEParameterType for QualifiedName {} + +/// Builder for polynomial ODE systems. +/// +/// This struct is just a convenient interface to construct a model of the theory of polynomial ODE +/// systems. Being an ordinary mutable Rust struct, it does *not* constitute a declarative language +/// to define ODE semantics for models of other theories. However, the idea is that it should be +/// used in a style that can mechanically translated to a future declarative language for model +/// migration. +#[derive(Clone)] +pub struct PolynomialODESystemBuilder { + model: ModalDblModel, + associated_parameters: HashMap, +} + +impl Default for PolynomialODESystemBuilder

{ + fn default() -> Self { + let th = th_signed_polynomial_ode_system(); + Self { + model: ModalDblModel::new(th.into()), + associated_parameters: HashMap::new(), + } + } +} + +impl PolynomialODESystemBuilder

{ + /// Constructs an empty ODE system. + pub fn new() -> Self { + Self::default() + } + + /// Returns a model of the theory of polynomial ODE systems. + pub fn model(self) -> ModalDblModel { + self.model + } + + /// Returns the HashMap of associated parameters, giving the term of type `P: ODEParameterType` + /// corresponding to each monomial. + pub fn associated_parameters(self) -> HashMap { + self.associated_parameters + } + + /// Adds a state variable to the ODE system. + pub fn add_variable(&mut self, var: QualifiedName) { + self.model.add_ob(var, ModeApp::new(name("State"))); + } + + /// Adds a contribution to the ODE system. + pub fn add_contribution( + &mut self, + id: QualifiedName, + target: QualifiedName, + sign: ContributionSign, + parameter: P, + monomial: impl IntoIterator, + ) { + let monomial = monomial.into_iter().map(ModalOb::Generator).collect(); + let sign = match sign { + ContributionSign::Positive => ModeApp::new(name("Contribution")).into(), + ContributionSign::Negative => ModeApp::new(name("NegativeContribution")).into(), + }; + + self.model.add_mor( + id.clone(), + ModalOb::List(List::Symmetric, monomial), + ModalOb::Generator(target), + sign, + ); + + self.associated_parameters.insert(id, parameter); + } +} + +/// This trait is where we define the actual ODE semantics, in the implementation of +/// `build_system_builder()`; `build_system()` will almost certainly always use the default +/// implementation given below. +/// +/// Note that the type that implements this trait is also where you are expected to state everything +/// that your semantics "cares about". For example, the default minimum is to give the values of +/// `ObType` and `MorType` that you want to distinguish between and iterate over. It can also hold +/// any extra data upon which your semantics can depend (see e.g. +/// `ode::mass_action::PetriNetMassActionAnalysis`, which contains the data of some +/// `MassConservationType`, whose value is fundamental in constructing the semantics). However, +/// this is left to the user: the type checker will *not* enforce any of these extras. +pub trait ODESemanticsAnalysis: Default { + /// The implementation of this function is what contains the actual data of the ODE semantics, + /// in the form of a `PolynomialODESystemBuilder`. + fn build_system_builder(&self, model: &T) -> PolynomialODESystemBuilder

; + + /// We simply feed the `PolynomialODESystemBuilder` constructed by the above function into + /// `PolynomialODEAnalysis::build_system_custom_parameters`. + fn build_system(&self, model: &T) -> PolynomialSystem, i8> { + let builder = self.build_system_builder(model); + PolynomialODEAnalysis::default().build_system_custom_parameters( + &builder.clone().model(), + builder.associated_parameters(), + ) + } +} + +/// A contribution to the ODE system consists of all the data that `ModalDblModel::add_mor()` +/// requires to create a multimorphism. +#[derive(Clone)] +pub struct Contribution { + /// The name of the multimorphism. + pub id: QualifiedName, + /// The target of the multimorphism, to be interpreted as the variable whose + /// first derivative is affected by the monomial. + pub target: QualifiedName, + /// The sign of a contribution. + pub sign: ContributionSign, + /// The parameter (coefficient) to be associated with this contribution. + pub parameter: P, + /// The source of the multimorphism (a list of objects), to be interpreted + /// as the monomial given by the product of all the list elements. + pub monomial: Vec, +} + +/// The sign of a contribution, since we work in *signed* multicategories. +#[derive(PartialEq, Eq, PartialOrd, Ord, Clone, Copy)] +pub enum ContributionSign { + /// Positive contribution: (d/dt)y -= x. + Positive, + /// Negative contribution: (d/dt)y += x. + Negative, +} + +/// The trait describing how to turn the formal system of ODEs into a numerical problem, to be +/// solved by an ODE solver and presented to the front-end. At minimum, such data must contain +/// initial values for variables and the intended duration of simulation, as well as the method for +/// converting the parameters (which are of type `ODEParameterType`) into floats. +// REQUEST | If you look at a struct that implements this trait (such as `LotkaVolterraProblemData`), +// FOR | there are a lot of serde statements going on. Should I be able to just move them +// FEEDBACK | (that is, those that come *before* the struct) here and have things all work? I'm still +// _________/ a bit intimidated by all these `crg_attr(feature = "serde")` bits. +// +// #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] +// #[cfg_attr(feature = "serde-wasm", derive(Tsify))] +// #[cfg_attr( +// feature = "serde-wasm", +// tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object) +// )] +pub trait ODESemanticsProblemData { + // REQUEST | The two getters (`initial_values()` and `duration()`) are annoying boilerplate to + // FOR | ask to be implemented. Is there a nice way to get rid of them here? Without them, + // FEEDBACK | the call to `self.initial_values` in `build_analysis()` fails because there is no + // _________/ way of knowing whether a struct implementing this trait actually has those fields. + /// Map from object IDs to initial values (nonnegative reals). + fn initial_values(&self) -> HashMap; + /// Duration of simulation. + fn duration(&self) -> f32; + + /// How to convert the formal parameters of type `ODEParameterType` into floats using values that + /// will eventually be filled in by the user from the front-end. + fn extend_scalars( + &self, + sys: PolynomialSystem, i8>, + ) -> PolynomialSystem; + + /// Converting the polynomial system into a system ready for use in numerical solvers. The default + /// implementation here should essentially always be the desired one. + fn build_analysis( + &self, + sys: PolynomialSystem, + ) -> ODEAnalysis> { + let ob_index: IndexMap<_, _> = + sys.components.keys().cloned().enumerate().map(|(i, x)| (x, i)).collect(); + let n = ob_index.len(); + + let initial_values = ob_index + .keys() + .map(|ob| self.initial_values().get(ob).copied().unwrap_or_default()); + let x0 = DVector::from_iterator(n, initial_values); + + let num_sys = sys.to_numerical(); + let problem = ODEProblem::new(num_sys, x0).end_time(self.duration()); + + ODEAnalysis::new(problem, ob_index) + } +} diff --git a/packages/catlog/src/stdlib/analyses/ode/polynomial_ode.rs b/packages/catlog/src/stdlib/analyses/ode/polynomial_ode.rs index 0b9d49f9f..d1fcfa08e 100644 --- a/packages/catlog/src/stdlib/analyses/ode/polynomial_ode.rs +++ b/packages/catlog/src/stdlib/analyses/ode/polynomial_ode.rs @@ -1,5 +1,17 @@ //! ODE analysis of models of the logic of systems of polynomial ODEs. -use std::collections::HashMap; +//! +//! This is used for the the simulation and equations analyses for models in the theory of +//! systems of polynomial ODEs [`th_polynomial_ode_system()`]. However, *all* ODE analyses +//! now factor through this by implementing [`ode::ode_semantics::ODESemantics`]; for further +//! documentation, see there. +//! +//! The interpretation of multicategories as systems of polynomial ODEs is explained in [RFC-0001]. +//! +//! [`th_polynomial_ode_system()`]: crate::stdlib::theories +//! [`ode::ode_semantics::ODESemantics`]: crate::stdlib::analyses::ode::ode_semantics::ODESemantics +//! [RFC-0001]: https://next.catcolab.org/rfc/0001 + +use std::{collections::HashMap, fmt}; use indexmap::IndexMap; use nalgebra::DVector; @@ -64,11 +76,37 @@ impl Default for PolynomialODEAnalysis { } impl PolynomialODEAnalysis { - /// Creates a system with symbolic coefficients. + /// Creates a `PolynomialSystem` with symbolic coefficients of type `QualifiedName`. pub fn build_system( &self, model: &ModalDblModel, ) -> PolynomialSystem, i8> { + // The default is to build a system whose parameters are in bijective correspondence + // with morphisms, given by using the `QualifiedName` of the morphism as the parameter + // generator. We thus build the graph of the identity function to pass as the HashMap + // of associated parameters. + let mut associated_parameters: HashMap = HashMap::new(); + for mor in model.mor_generators_with_type(&self.positive_contribution_mor_type) { + associated_parameters.insert(mor.clone(), mor.clone()); + } + for mor in model.mor_generators_with_type(&self.negative_contribution_mor_type) { + associated_parameters.insert(mor.clone(), mor.clone()); + } + + self.build_system_custom_parameters::(model, associated_parameters) + } + + /// Creates a `PolynomialSystem` with symbolic coefficients of some generic type. + /// + /// When constructing a system as a derived model from another model (as in e.g. `mass_action`), + /// it is not necessarily the case that each morphism will give rise to a unique parameter. This + /// function allows for the construction of a `PolynomialSystem<_ , Parameter, _>` using some + /// specified `HashMap` that describes how to associate parameters to morphisms. + pub fn build_system_custom_parameters( + &self, + model: &ModalDblModel, + associated_parameters: HashMap, + ) -> PolynomialSystem, i8> { let mut sys = PolynomialSystem::new(); // Create a variable for each object. @@ -76,17 +114,31 @@ impl PolynomialODEAnalysis { sys.add_term(ob, Polynomial::zero()); } + // Every morphism will give a term, i.e. a pair consisting of a monomial and a parameter. + // Although the *monomial* depends only on the input objects to the morphism, the *parameter* + // might be described by external data. For example, multiple morphisms might share the same + // parameter. + // + // This closure builds a term to add to the `PolynomialSystem` given a morphism and the + // hash map `associated_parameters`. let make_term = |mor: QualifiedName| { + // Find the inputs and output of the morphism. let (Some(ModalOb::List(_, inputs)), Some(output)) = (model.get_dom(&mor), model.get_cod(&mor)) else { return None; }; - let term: Monomial<_, _> = + // Construct the monomial given by the product of all of the inputs. + let monomial: Monomial<_, _> = inputs.iter().cloned().map(|ob| (ob.unwrap_generator(), 1)).collect(); - let term: Polynomial<_, _, _> = - [(Parameter::generator(mor), term.clone())].into_iter().collect(); + // Construct the term given by the monomial and the parameter from `associated_parameters`. + let term: Polynomial<_, _, _> = [( + Parameter::generator(associated_parameters.get(&mor).unwrap().clone()), + monomial.clone(), + )] + .into_iter() + .collect(); Some((output.clone().unwrap_generator(), term)) }; @@ -97,7 +149,6 @@ impl PolynomialODEAnalysis { sys.add_term(var, term); } } - // Add a monomial with negative sign for each negative contribution. for mor in model.mor_generators_with_type(&self.negative_contribution_mor_type) { if let Some((var, term)) = make_term(mor) { @@ -153,11 +204,11 @@ mod tests { tt, }; - // (Unsigned) Lotka–Volterra dynamics on a two-level model. + /// (Unsigned) Lotka-Volterra dynamics on a two-level model. #[test] - fn lotka_volterra_equations() { + fn unsigned_lotka_volterra_equations() { let th = Rc::new(th_polynomial_ode_system()); - let model = lotka_volterra_dynamics(th); + let model = unsigned_lotka_volterra_dynamics(th); let sys = PolynomialODEAnalysis::default().build_system(&model); let expected = expect!([r#" dA = A_growth A + BA_interaction A B @@ -167,31 +218,26 @@ mod tests { expected.assert_eq(&sys.to_string()); } - // (Unsigned) Lotka–Volterra dynamics on a two-level model with LaTeX. + /// Lotka-Volterra dynamics on a two-level model with LaTeX. #[test] fn lotka_volterra_equations_latex() { - let th = Rc::new(th_polynomial_ode_system()); - let model = lotka_volterra_dynamics(th); + let th = Rc::new(th_signed_polynomial_ode_system()); + let model = signed_lotka_volterra_dynamics(th); let sys = PolynomialODEAnalysis::default().build_system(&model); let expected = vec![ LatexEquation { lhs: "\\frac{\\mathrm{d}}{\\mathrm{d}t} A".to_string(), - rhs: "A_growth \\cdot A + BA_interaction \\cdot A \\cdot B".to_string(), + rhs: "A_growth \\cdot A - BA_interaction \\cdot A \\cdot B".to_string(), }, LatexEquation { lhs: "\\frac{\\mathrm{d}}{\\mathrm{d}t} B".to_string(), - rhs: "AB_interaction \\cdot A \\cdot B + B_growth \\cdot B + CB_interaction \\cdot B \\cdot C" - .to_string(), - }, - LatexEquation { - lhs: "\\frac{\\mathrm{d}}{\\mathrm{d}t} C".to_string(), - rhs: "BC_interaction \\cdot B \\cdot C + C_growth \\cdot C".to_string(), + rhs: "AB_interaction \\cdot A \\cdot B + B_growth \\cdot B".to_string(), }, ]; assert_eq!(expected, sys.to_latex_equations()); } - // DoubleTT elaboration from text. + /// DoubleTT elaboration from text. #[test] fn ode_system_from_text() { let th = Rc::new(th_polynomial_ode_system()); diff --git a/packages/catlog/src/stdlib/analyses/ode/signed_coefficients.rs b/packages/catlog/src/stdlib/analyses/ode/signed_coefficients.rs deleted file mode 100644 index ce106e909..000000000 --- a/packages/catlog/src/stdlib/analyses/ode/signed_coefficients.rs +++ /dev/null @@ -1,84 +0,0 @@ -//! Helper module to build analyses based on signed coefficient matrices. - -use indexmap::IndexMap; -use nalgebra::DMatrix; -use num_traits::zero; - -use super::Parameter; -use crate::{ - dbl::model::FpDblModel, - zero::{QualifiedName, rig::Monomial}, -}; - -/// Builder for signed coefficient matrices and analyses based on them. -/// -/// Used to construct the [linear](Self::linear_ode_analysis) and -/// [Lotka-Volterra](Self::lotka_volterra_analysis) ODE analyses. -pub struct SignedCoefficientBuilder { - var_ob_type: ObType, - positive_mor_types: Vec, - negative_mor_types: Vec, -} - -impl SignedCoefficientBuilder { - /// Creates a new builder for the given object type. - pub fn new(var_ob_type: ObType) -> Self { - Self { - var_ob_type, - positive_mor_types: Vec::new(), - negative_mor_types: Vec::new(), - } - } - - /// Adds a morphism type defining a positive interaction between objects. - pub fn add_positive(mut self, mor_type: MorType) -> Self { - self.positive_mor_types.push(mor_type); - self - } - - /// Adds a morphism type defining a negative interaction between objects. - pub fn add_negative(mut self, mor_type: MorType) -> Self { - self.negative_mor_types.push(mor_type); - self - } - - /// 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 FpDblModel< - ObType = ObType, - MorType = MorType, - Ob = QualifiedName, - ObGen = QualifiedName, - MorGen = QualifiedName, - >, - ) -> (DMatrix>, IndexMap) { - let ob_index: IndexMap<_, _> = model - .ob_generators_with_type(&self.var_ob_type) - .enumerate() - .map(|(i, x)| (x, i)) - .collect(); - - let n = ob_index.len(); - let mut mat = DMatrix::from_element(n, n, 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)] += (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)] += (-1.0, Monomial::generator(mor)); - } - } - - (mat, ob_index) - } -} diff --git a/packages/catlog/src/stdlib/analyses/petri.rs b/packages/catlog/src/stdlib/analyses/petri.rs index ec28171ca..f00fecaa8 100644 --- a/packages/catlog/src/stdlib/analyses/petri.rs +++ b/packages/catlog/src/stdlib/analyses/petri.rs @@ -1,21 +1,49 @@ //! Helpers for analyses on Petri nets. -use crate::dbl::model::{ModalDblModel, ModalOb, MutDblModel}; +use crate::dbl::model::{ModalDblModel, MutDblModel}; use crate::dbl::theory::Unital; use crate::zero::QualifiedName; +pub struct TransitionInterface { + pub input_places: Vec, + pub output_places: Vec, +} + +// TODO: Unfortunately, in the case of transition_interface, there is a further +// subtlety that isn't addressed by these considerations. The collect_product +// function only collects one level of operation application, as opposed to +// acting recursively. Thus, I'd say it's technically incorrect to unwrap +// generators from the lists returned. This point is a bit academic since in +// the notebook editor you couldn't construct such a model anyway, but it is +// perfectly valid in the text elaborator to write tensor [a, tensor [b, c]] +// and we shouldn't bomb on that. +// +// To do this safely, you should collect recursively rather than at one level; +// however, under the validation assumption, you are allowed (in fact +// encouraged) to panic if you encounter anything that is not an basic object +// or an application of tensor to a list. + /// Gets the inputs and outputs of a transition in a Petri net. pub fn transition_interface( model: &ModalDblModel, id: &QualifiedName, -) -> (Vec, Vec) { +) -> TransitionInterface { let inputs = model .get_dom(id) .and_then(|ob| ob.clone().collect_product(None)) - .unwrap_or_default(); + .unwrap_or_default() + .into_iter() + .map(|ob| ob.unwrap_generator()) + .collect(); let outputs = model .get_cod(id) .and_then(|ob| ob.clone().collect_product(None)) - .unwrap_or_default(); - (inputs, outputs) + .unwrap_or_default() + .into_iter() + .map(|ob| ob.unwrap_generator()) + .collect(); + TransitionInterface { + input_places: inputs, + output_places: outputs, + } } diff --git a/packages/catlog/src/stdlib/analyses/reachability.rs b/packages/catlog/src/stdlib/analyses/reachability.rs index 402b5dac5..8a7c91028 100644 --- a/packages/catlog/src/stdlib/analyses/reachability.rs +++ b/packages/catlog/src/stdlib/analyses/reachability.rs @@ -3,10 +3,10 @@ use itertools::Itertools; use std::collections::HashMap; -use crate::dbl::modal::model::{ModalDblModel, ModalOb}; +use crate::dbl::modal::model::ModalDblModel; use crate::dbl::theory::Unital; use crate::one::category::FgCategory; -use crate::stdlib::analyses::petri::transition_interface; +use crate::stdlib::analyses::petri::{TransitionInterface, transition_interface}; use crate::zero::QualifiedName; #[cfg(feature = "serde")] @@ -57,16 +57,14 @@ pub fn subreachability(m: &ModalDblModel, data: ReachabilityProblemData) for e in m.mor_generators() { let e_idx = *hom_inv.get(&e).unwrap(); - let (inputs, outputs) = transition_interface(m, &e); + let transition_interface: TransitionInterface = transition_interface(m, &e); + let inputs = transition_interface.input_places.clone(); + let outputs = transition_interface.output_places.clone(); for ob in inputs { - if let ModalOb::Generator(u) = ob { - i_mat[*ob_inv.get(&u).unwrap()][e_idx] += 1; - } + i_mat[*ob_inv.get(&ob).unwrap()][e_idx] += 1; } for ob in outputs { - if let ModalOb::Generator(u) = ob { - o_mat[*ob_inv.get(&u).unwrap()][e_idx] += 1; - } + o_mat[*ob_inv.get(&ob).unwrap()][e_idx] += 1; } } let (i_mat_, o_mat_) = (&i_mat, &o_mat); diff --git a/packages/catlog/src/stdlib/analyses/stochastic/mass_action.rs b/packages/catlog/src/stdlib/analyses/stochastic/mass_action.rs index e0edf3094..5aabdb961 100644 --- a/packages/catlog/src/stdlib/analyses/stochastic/mass_action.rs +++ b/packages/catlog/src/stdlib/analyses/stochastic/mass_action.rs @@ -8,7 +8,10 @@ use std::collections::HashMap; use crate::{ dbl::{modal::*, model::FpDblModel, theory::Unital}, - stdlib::analyses::{ode::ODESolution, petri::transition_interface}, + stdlib::analyses::{ + ode::ODESolution, + petri::{TransitionInterface, transition_interface}, + }, zero::{QualifiedName, name}, }; @@ -108,20 +111,16 @@ impl PetriNetStochasticMassActionAnalysis { let mut problem = gillespie::Gillespie::new(initial, false); for mor in model.mor_generators_with_type(&self.transition_mor_type) { - let (inputs, outputs) = transition_interface(model, &mor); + let transition_interface: TransitionInterface = transition_interface(model, &mor); + let inputs = transition_interface.input_places.clone(); + let outputs = transition_interface.output_places.clone(); // 1. convert the inputs/outputs to sequences of counts let input_vec = ob_generators.iter().map(|id| { - inputs - .iter() - .filter(|&ob| matches!(ob, ModalOb::Generator(id2) if id2 == id)) - .count() as u32 + inputs.iter().filter(|&ob| matches!(ob, id2 if id2 == id)).count() as u32 }); let output_vec = ob_generators.iter().map(|id| { - outputs - .iter() - .filter(|&ob| matches!(ob, ModalOb::Generator(id2) if id2 == id)) - .count() as isize + outputs.iter().filter(|&ob| matches!(ob, id2 if id2 == id)).count() as isize }); // 2. output := output - input diff --git a/packages/catlog/src/stdlib/analyses/stock_flow.rs b/packages/catlog/src/stdlib/analyses/stock_flow.rs new file mode 100644 index 000000000..2cb83ee55 --- /dev/null +++ b/packages/catlog/src/stdlib/analyses/stock_flow.rs @@ -0,0 +1,46 @@ +//! Helpers for analyses on stock-flow diagrams. + +use crate::dbl::discrete_tabulator::DiscreteTabModel; +use crate::dbl::discrete_tabulator::TabEdge; +use crate::dbl::discrete_tabulator::TabMorType; +use crate::dbl::model::FpDblModel; +use crate::dbl::model::TabOb; +use crate::one::category::FgCategory; +use crate::zero::QualifiedName; +use crate::zero::name; + +pub struct FlowInterface { + pub input_stock: QualifiedName, + pub input_pos_link_doms: Vec, + pub output_stock: QualifiedName, +} + +/// Gets the inputs (including links) and output of a flow in a stock-flow diagram. +pub fn flow_interface(model: &DiscreteTabModel, flow: &QualifiedName) -> FlowInterface { + let dom = model.mor_generator_dom(flow).unwrap_basic(); + let cod = model.mor_generator_cod(flow).unwrap_basic(); + + let mut input_pos_link_doms: Vec = Vec::new(); + + // Iterate over positive links and add them to the interface if their codomain is the + // link in question. + for link in model.mor_generators_with_type(&TabMorType::Basic(name("Link"))) { + let dom = model.mor_generator_dom(&link); + let path = model.mor_generator_cod(&link).unwrap_tabulated(); + let Some(TabEdge::Basic(cod)) = path.only() else { + panic!("Codomain of link should be basic morphism"); + }; + if cod == *flow { + input_pos_link_doms.push(dom) + }; + } + + FlowInterface { + input_stock: dom, + input_pos_link_doms: input_pos_link_doms + .iter() + .map(|stock| stock.clone().unwrap_basic()) + .collect(), + output_stock: cod, + } +} diff --git a/packages/catlog/src/stdlib/models.rs b/packages/catlog/src/stdlib/models.rs index 6965d47f9..970332d15 100644 --- a/packages/catlog/src/stdlib/models.rs +++ b/packages/catlog/src/stdlib/models.rs @@ -167,14 +167,17 @@ pub fn sir_petri(th: Rc>) -> ModalDblModel { model } -/// An example of Lotka–Volterra dynamics viewed as a non-unital theory for a symmetric multicategory. -pub fn lotka_volterra_dynamics(th: Rc>) -> ModalDblModel { +/// An example of (unsigned) Lotka-Volterra dynamics viewed as a non-unital theory for +/// a symmetric multicategory. +pub fn unsigned_lotka_volterra_dynamics( + th: Rc>, +) -> ModalDblModel { let ob_type = ModalObType::new(name("State")); let mor_type: ModalMorType = ModeApp::new(name("Contribution")).into(); let mut model = ModalDblModel::new(th); - // We're going to build a two-level predator-prey model, but where (in absence of signed - // arrows) all interactions have positive coefficients. + // A two-level predator-prey model, but where (in absence of signed arrows) all + // interactions have positive coefficients. let (a, b, c) = (name("A"), name("B"), name("C")); model.add_ob(a.clone(), ob_type.clone()); @@ -235,6 +238,54 @@ pub fn lotka_volterra_dynamics(th: Rc>) -> ModalDblMod model } +/// An example of Lotka-Volterra dynamics viewed as a non-unital theory for a symmetric multicategory. +pub fn signed_lotka_volterra_dynamics( + th: Rc>, +) -> ModalDblModel { + let ob_type = ModalObType::new(name("State")); + let pos_mor_type: ModalMorType = ModeApp::new(name("Contribution")).into(); + let neg_mor_type: ModalMorType = ModeApp::new(name("NegativeContribution")).into(); + + let mut model = ModalDblModel::new(th); + // We're going to build a simple predator-prey model. + let (a, b) = (name("A"), name("B")); + + model.add_ob(a.clone(), ob_type.clone()); + model.add_ob(b.clone(), ob_type.clone()); + // The growth terms, corresponding to + // dA/dt += g_A A + // dB/dt += g_B B + model.add_mor( + name("A_growth"), + ModalOb::List(List::Symmetric, vec![a.clone().into()]), + a.clone().into(), + pos_mor_type.clone(), + ); + model.add_mor( + name("B_growth"), + ModalOb::List(List::Symmetric, vec![b.clone().into()]), + b.clone().into(), + pos_mor_type.clone(), + ); + // The interaction terms, corresponding to + // dB/dt += k_AB AB + // dA/dt -= k_BA AB + model.add_mor( + name("AB_interaction"), + ModalOb::List(List::Symmetric, vec![a.clone().into(), b.clone().into()]), + b.clone().into(), + pos_mor_type.clone(), + ); + model.add_mor( + name("BA_interaction"), + ModalOb::List(List::Symmetric, vec![a.clone().into(), b.clone().into()]), + a.clone().into(), + neg_mor_type.clone(), + ); + + model +} + #[cfg(test)] mod tests { use super::super::theories::*; @@ -288,6 +339,6 @@ mod tests { #[test] fn polynomial_ode_systems() { let th = Rc::new(th_polynomial_ode_system()); - assert!(lotka_volterra_dynamics(th.clone()).validate().is_ok()); + assert!(unsigned_lotka_volterra_dynamics(th.clone()).validate().is_ok()); } } diff --git a/packages/catlog/src/stdlib/theories.rs b/packages/catlog/src/stdlib/theories.rs index a0a6e3275..0c1a13070 100644 --- a/packages/catlog/src/stdlib/theories.rs +++ b/packages/catlog/src/stdlib/theories.rs @@ -379,6 +379,7 @@ mod tests { assert!(th_sym_multicategory().validate().is_ok()); assert!(modal_th_power_system().validate().is_ok()); assert!(th_polynomial_ode_system().validate().is_ok()); + assert!(th_signed_polynomial_ode_system().validate().is_ok()); } #[test] diff --git a/packages/catlog/src/zero/qualified.rs b/packages/catlog/src/zero/qualified.rs index 908f76d23..60f13c8bf 100644 --- a/packages/catlog/src/zero/qualified.rs +++ b/packages/catlog/src/zero/qualified.rs @@ -294,6 +294,13 @@ impl QualifiedName { } } + /// Prepend a name segment. + pub fn cons(&self, segment: NameSegment) -> Self { + let mut segments = self.0.clone(); + segments.insert(0, segment); + Self(segments) + } + /// Add another segment onto the end. pub fn snoc(&self, segment: NameSegment) -> Self { let mut segments = self.0.clone(); diff --git a/packages/frontend/src/help/analysis/mass-action.mdx b/packages/frontend/src/help/analysis/mass-action.mdx index 68286d712..29c298088 100644 --- a/packages/frontend/src/help/analysis/mass-action.mdx +++ b/packages/frontend/src/help/analysis/mass-action.mdx @@ -7,12 +7,12 @@

Whether or not flows should preserve mass
Rate: $\mathbb{R}_{\geqslant0}$
*(Only if **Mass conservation** = `True`)* The rate coefficient ($r$) of the reaction
-
Rate granularity: `Per transition | Per place`
+
Rate granularity: `Per flow | Per stock`
*(Only if **Mass conservation** = `False`)* If flows can have multiple inputs/outputs (e.g. in the case of Petri nets) then rates can be given per flow or per individual input/output
Consumption: $\mathbb{R}_{\geqslant0}$
-
The consumption rate coefficient ($\kappa$), either per transition (flow) or per place (input/output) depending on **Mass conservation**
+
The consumption rate coefficient ($\kappa$), either per flow or per stock (input/output) depending on **Mass conservation**
Production: $\mathbb{R}_{\geqslant0}$
-
The production rate coefficient ($\rho$), either per transition (flow) or per place (input/output) depending on **Mass conservation**
+
The production rate coefficient ($\rho$), either per flow or per stock (input/output) depending on **Mass conservation**
Duration: $\mathbb{R}_{\geqslant0}$
The total duration of the simulation in units of time
diff --git a/packages/frontend/src/help/logics/petri-net.mdx b/packages/frontend/src/help/logics/petri-net.mdx index 5c652ca8f..3ac06d8ac 100644 --- a/packages/frontend/src/help/logics/petri-net.mdx +++ b/packages/frontend/src/help/logics/petri-net.mdx @@ -53,7 +53,7 @@ The rest of the analysis depends on whether **mass conservation** is checked as - $\dot{X}=r_T AB$ - $\dot{Y}=r_T AB$ -- If **mass conservation** is checked as _false_, and **rate granularity** is set to _per transition_, then a transition $A\xrightarrow{T}B$ between places $A$ and $B$ is interpreted as the equations $\{\dot{A}=-\kappa_T,\dot{B}=\rho_T\}$. Here $\kappa_T$ and $\rho_T$ are the **consumption** and **production** rate coefficients of the transition. +- If **mass conservation** is checked as _false_, and **rate granularity** is set to _per flow_, then a transition $A\xrightarrow{T}B$ between places $A$ and $B$ is interpreted as the equations $\{\dot{A}=-\kappa_T,\dot{B}=\rho_T\}$. Here $\kappa_T$ and $\rho_T$ are the **consumption** and **production** rate coefficients of the transition. A transition $[A,B]\xrightarrow{T}[X,Y]$ gives the equations - $\dot{A}=-\kappa_T AB$ @@ -61,7 +61,7 @@ The rest of the analysis depends on whether **mass conservation** is checked as - $\dot{X}=\rho_T AB$ - $\dot{Y}=\rho_T AB$ -- If **mass conservation** is checked as _false_, and **rate granularity** is set to _per place_, then a transition $A\xrightarrow{T}B$ between places $A$ and $B$ is interpreted as the equations $\{\dot{A}=-\kappa_T^A,\dot{B}=\rho_T^B\}$. Here $\kappa_T^A$ and $\rho_T^B$ are the **consumption** and **production** rate coefficients of the objects $A$ and $B$ with respect to the transition $T$. +- If **mass conservation** is checked as _false_, and **rate granularity** is set to _per stock_, then a transition $A\xrightarrow{T}B$ between places $A$ and $B$ is interpreted as the equations $\{\dot{A}=-\kappa_T^A,\dot{B}=\rho_T^B\}$. Here $\kappa_T^A$ and $\rho_T^B$ are the **consumption** and **production** rate coefficients of the objects $A$ and $B$ with respect to the transition $T$. A transition $[A,B]\xrightarrow{T}[X,Y]$ gives the equations - $\dot{A}=-\kappa_T^A AB$ diff --git a/packages/frontend/src/stdlib/analyses.tsx b/packages/frontend/src/stdlib/analyses.tsx index 2bb74794c..8d4d4d0ab 100644 --- a/packages/frontend/src/stdlib/analyses.tsx +++ b/packages/frontend/src/stdlib/analyses.tsx @@ -1,6 +1,8 @@ import { lazy } from "solid-js"; import type { + LinearODEEquationsData, + LotkaVolterraEquationsData, MassActionEquationsData, MorType, ObType, @@ -133,6 +135,33 @@ export function linearODE( const LinearODE = lazy(() => import("./analyses/linear_ode")); +export function linearODEEquations( + options: Partial & { + getEquations: Simulators.LinearODEEquations; + }, +): ModelAnalysisMeta { + const { + id = "linear-ode-equations", + name = "Linear ODE equations", + description = "Display the symbolic linear ODE dynamics equations", + help = "linear-ode-equations", + ...otherOptions + } = options; + return { + id, + name, + description, + help, + component: (props) => ( + + ), + initialContent: () => ({ + trivialData: true, + }), + }; +} +const LinearODEEquationsDisplay = lazy(() => import("./analyses/linear_ode_equations")); + export function lotkaVolterra( options: Partial & { simulate: Simulators.LotkaVolterraSimulator; @@ -140,8 +169,8 @@ export function lotkaVolterra( ): ModelAnalysisMeta { const { id = "lotka-volterra", - name = "Lotka-Volterra dynamics", - description = "Simulate the system using a Lotka-Volterra ODE", + name = "Lotka–Volterra dynamics", + description = "Simulate the system using a Lotka–Volterra ODE", help = "lotka-volterra", simulate, } = options; @@ -162,6 +191,33 @@ export function lotkaVolterra( const LotkaVolterra = lazy(() => import("./analyses/lotka_volterra")); +export function lotkaVolterraEquations( + options: Partial & { + getEquations: Simulators.LotkaVolterraEquations; + }, +): ModelAnalysisMeta { + const { + id = "lotka-volterra-equations", + name = "Lotka–Volterra equations", + description = "Display the symbolic Lotka–Volterra dynamics equations", + help = "lotka-volterra-equations", + ...otherOptions + } = options; + return { + id, + name, + description, + help, + component: (props) => ( + + ), + initialContent: () => ({ + trivialData: true, + }), + }; +} +const LotkaVolterraEquationsDisplay = lazy(() => import("./analyses/lotka_volterra_equations")); + export function massAction( options: Partial & { ratesHaveGranularity: boolean; diff --git a/packages/frontend/src/stdlib/analyses/linear_ode.tsx b/packages/frontend/src/stdlib/analyses/linear_ode.tsx index 946a156ca..94c42538b 100644 --- a/packages/frontend/src/stdlib/analyses/linear_ode.tsx +++ b/packages/frontend/src/stdlib/analyses/linear_ode.tsx @@ -4,12 +4,14 @@ import { createNumericalColumn, FixedTableEditor, Foldable, + ExpandableTable, + KatexDisplay, } from "catcolab-ui-components"; -import type { DblModel, LinearODEProblemData, QualifiedName } from "catlog-wasm"; +import type { LinearODEProblemData, QualifiedName } from "catlog-wasm"; import type { ModelAnalysisProps } from "../../analysis"; import { morLabelOrDefault } from "../../model"; import { ODEResultPlot } from "../../visualization"; -import { createModelODEPlot } from "./model_ode_plot"; +import { createModelODEPlotWithEquations } from "./model_ode_plot"; import type { LinearODESimulator } from "./simulator_types"; import "./simulation.css"; @@ -70,11 +72,14 @@ export default function LinearODE( }), ]; - const plotResult = createModelODEPlot( + const result = createModelODEPlotWithEquations( () => props.liveModel.validatedModel(), - (model: DblModel) => props.simulate(model, props.content), + (model) => props.simulate(model, props.content), ); + const plotResult = () => result()?.plotData; + const latexEquations = () => result()?.latexEquations ?? []; + return (
@@ -91,7 +96,20 @@ export default function LinearODE(
- + + }, + { cell: () => }, + { cell: (row) => }, + ]} + /> + + + + ); } diff --git a/packages/frontend/src/stdlib/analyses/linear_ode_equations.tsx b/packages/frontend/src/stdlib/analyses/linear_ode_equations.tsx new file mode 100644 index 000000000..d0e65d2db --- /dev/null +++ b/packages/frontend/src/stdlib/analyses/linear_ode_equations.tsx @@ -0,0 +1,36 @@ +import { BlockTitle, ExpandableTable, KatexDisplay } from "catcolab-ui-components"; +import { LinearODEEquationsData } from "catlog-wasm"; +import type { ModelAnalysisProps } from "../../analysis"; +import { createModelODELatex } from "./model_ode_plot"; +import type { LinearODEEquations } from "./simulator_types"; + +import "./simulation.css"; + +/** Display the symbolic mass-action dynamics equations for a model. */ +export default function LinearODEEquationsDisplay( + props: ModelAnalysisProps & { + content: LinearODEEquationsData; + getEquations: LinearODEEquations; + title?: string; + }, +) { + const latexEquations = createModelODELatex( + () => props.liveModel.validatedModel(), + (model) => props.getEquations(model, props.content), + ); + + return ( +
+ + }, + { cell: () => }, + { cell: (row) => }, + ]} + /> +
+ ); +} diff --git a/packages/frontend/src/stdlib/analyses/lotka_volterra.tsx b/packages/frontend/src/stdlib/analyses/lotka_volterra.tsx index 9d6006800..062f28189 100644 --- a/packages/frontend/src/stdlib/analyses/lotka_volterra.tsx +++ b/packages/frontend/src/stdlib/analyses/lotka_volterra.tsx @@ -4,12 +4,14 @@ import { createNumericalColumn, FixedTableEditor, Foldable, + ExpandableTable, + KatexDisplay, } from "catcolab-ui-components"; -import type { DblModel, LotkaVolterraProblemData, QualifiedName } from "catlog-wasm"; +import type { LotkaVolterraProblemData, QualifiedName } from "catlog-wasm"; import type { ModelAnalysisProps } from "../../analysis"; import { morLabelOrDefault } from "../../model"; import { ODEResultPlot } from "../../visualization"; -import { createModelODEPlot } from "./model_ode_plot"; +import { createModelODEPlotWithEquations } from "./model_ode_plot"; import type { LotkaVolterraSimulator } from "./simulator_types"; import "./simulation.css"; @@ -78,11 +80,14 @@ export default function LotkaVolterra( }), ]; - const plotResult = createModelODEPlot( + const result = createModelODEPlotWithEquations( () => props.liveModel.validatedModel(), - (model: DblModel) => props.simulate(model, props.content), + (model) => props.simulate(model, props.content), ); + const plotResult = () => result()?.plotData; + const latexEquations = () => result()?.latexEquations ?? []; + return (
@@ -99,7 +104,20 @@ export default function LotkaVolterra(
- + + }, + { cell: () => }, + { cell: (row) => }, + ]} + /> + + + + ); } diff --git a/packages/frontend/src/stdlib/analyses/lotka_volterra_equations.tsx b/packages/frontend/src/stdlib/analyses/lotka_volterra_equations.tsx new file mode 100644 index 000000000..dcb6271d4 --- /dev/null +++ b/packages/frontend/src/stdlib/analyses/lotka_volterra_equations.tsx @@ -0,0 +1,36 @@ +import { BlockTitle, ExpandableTable, KatexDisplay } from "catcolab-ui-components"; +import { LotkaVolterraEquationsData } from "catlog-wasm"; +import type { ModelAnalysisProps } from "../../analysis"; +import { createModelODELatex } from "./model_ode_plot"; +import type { LotkaVolterraEquations } from "./simulator_types"; + +import "./simulation.css"; + +/** Display the symbolic mass-action dynamics equations for a model. */ +export default function LotkaVolterraEquationsDisplay( + props: ModelAnalysisProps & { + content: LotkaVolterraEquationsData; + getEquations: LotkaVolterraEquations; + title?: string; + }, +) { + const latexEquations = createModelODELatex( + () => props.liveModel.validatedModel(), + (model) => props.getEquations(model, props.content), + ); + + return ( +
+ + }, + { cell: () => }, + { cell: (row) => }, + ]} + /> +
+ ); +} diff --git a/packages/frontend/src/stdlib/analyses/mass_action_config_form.tsx b/packages/frontend/src/stdlib/analyses/mass_action_config_form.tsx index b16365db0..84332f3be 100644 --- a/packages/frontend/src/stdlib/analyses/mass_action_config_form.tsx +++ b/packages/frontend/src/stdlib/analyses/mass_action_config_form.tsx @@ -51,8 +51,8 @@ export function MassActionConfigForm(props: { }); }} > - - + + diff --git a/packages/frontend/src/stdlib/analyses/simulator_types.ts b/packages/frontend/src/stdlib/analyses/simulator_types.ts index e5ac07d98..335b5e7ed 100644 --- a/packages/frontend/src/stdlib/analyses/simulator_types.ts +++ b/packages/frontend/src/stdlib/analyses/simulator_types.ts @@ -3,7 +3,9 @@ import type { KuramotoProblemData, LatexEquations, LinearODEProblemData, + LinearODEEquationsData, LotkaVolterraProblemData, + LotkaVolterraEquationsData, MassActionEquationsData, MassActionProblemData, ODEResult, @@ -22,20 +24,31 @@ export type { }; export type KuramotoSimulator = (model: DblModel, data: KuramotoProblemData) => ODEResult; -export type LinearODESimulator = (model: DblModel, data: LinearODEProblemData) => ODEResult; -export type LotkaVolterraSimulator = (model: DblModel, data: LotkaVolterraProblemData) => ODEResult; +export type LinearODESimulator = ( + model: DblModel, + data: LinearODEProblemData, +) => ODEResultWithEquations; +export type LinearODEEquations = (model: DblModel, data: LinearODEEquationsData) => LatexEquations; +export type LotkaVolterraSimulator = ( + model: DblModel, + data: LotkaVolterraProblemData, +) => ODEResultWithEquations; +export type LotkaVolterraEquations = ( + model: DblModel, + data: LotkaVolterraEquationsData, +) => LatexEquations; export type MassActionSimulator = ( model: DblModel, data: MassActionProblemData, ) => ODEResultWithEquations; -export type StochasticMassActionSimulator = ( - model: DblModel, - data: StochasticMassActionProblemData, -) => ODEResult; export type MassActionEquations = ( model: DblModel, data: MassActionEquationsData, ) => LatexEquations; +export type StochasticMassActionSimulator = ( + model: DblModel, + data: StochasticMassActionProblemData, +) => ODEResult; export type PolynomialODESimulator = ( model: DblModel, data: PolynomialODEProblemData, diff --git a/packages/frontend/src/stdlib/theories/causal-loop.ts b/packages/frontend/src/stdlib/theories/causal-loop.ts index 174d40108..8d81ad97b 100644 --- a/packages/frontend/src/stdlib/theories/causal-loop.ts +++ b/packages/frontend/src/stdlib/theories/causal-loop.ts @@ -76,9 +76,19 @@ export default function createCausalLoopTheory(theoryMeta: TheoryMeta): Theory { analyses.linearODE({ simulate: (model, data) => thSignedCategory.linearODE(model, data), }), + analyses.linearODEEquations({ + getEquations(model) { + return thSignedCategory.linearODEEquations(model); + }, + }), analyses.lotkaVolterra({ simulate: (model, data) => thSignedCategory.lotkaVolterra(model, data), }), + analyses.lotkaVolterraEquations({ + getEquations(model) { + return thSignedCategory.lotkaVolterraEquations(model); + }, + }), ], }); } diff --git a/packages/frontend/src/stdlib/theories/primitive-signed-stock-flow.ts b/packages/frontend/src/stdlib/theories/primitive-signed-stock-flow.ts index b99784d6f..a9774c5cb 100644 --- a/packages/frontend/src/stdlib/theories/primitive-signed-stock-flow.ts +++ b/packages/frontend/src/stdlib/theories/primitive-signed-stock-flow.ts @@ -68,22 +68,6 @@ export default function createPrimitiveSignedStockFlowTheory(theoryMeta: TheoryM description: "Visualize the stock and flow diagram", help: "visualization", }), - analyses.massAction({ - ratesHaveGranularity: false, - simulate(model, data) { - return thCategorySignedLinks.massAction(model, data); - }, - transitionType: { - tag: "Hom", - content: { tag: "Basic", content: "Object" }, - }, - }), - analyses.massActionEquations({ - ratesHaveGranularity: false, - getEquations(model, data) { - return thCategorySignedLinks.massActionEquations(model, data); - }, - }), ], }); } diff --git a/packages/frontend/src/stdlib/theories/reg-net.ts b/packages/frontend/src/stdlib/theories/reg-net.ts index ff6751faf..29f627c9c 100644 --- a/packages/frontend/src/stdlib/theories/reg-net.ts +++ b/packages/frontend/src/stdlib/theories/reg-net.ts @@ -75,9 +75,17 @@ export default function createRegulatoryNetworkTheory(theoryMeta: TheoryMeta): T analyses.linearODE({ simulate: (model, data) => thSignedCategory.linearODE(model, data), }), + analyses.linearODEEquations({ + getEquations(model) { + return thSignedCategory.linearODEEquations(model); + }, + }), analyses.lotkaVolterra({ - simulate(model, data) { - return thSignedCategory.lotkaVolterra(model, data); + simulate: (model, data) => thSignedCategory.lotkaVolterra(model, data), + }), + analyses.lotkaVolterraEquations({ + getEquations(model) { + return thSignedCategory.lotkaVolterraEquations(model); }, }), ],