From bfef835aedcc8f44f9d5aa318274d28cf4d6f21a Mon Sep 17 00:00:00 2001 From: Matt Cuffaro Date: Wed, 23 Jul 2025 15:27:20 -0700 Subject: [PATCH 1/6] ENH: added theory and model --- packages/catlog/src/stdlib/models.rs | 37 +++++++++++++++++++++++++- packages/catlog/src/stdlib/theories.rs | 24 ++++++++++++++++- 2 files changed, 59 insertions(+), 2 deletions(-) diff --git a/packages/catlog/src/stdlib/models.rs b/packages/catlog/src/stdlib/models.rs index e0fc91bac..e596c3ade 100644 --- a/packages/catlog/src/stdlib/models.rs +++ b/packages/catlog/src/stdlib/models.rs @@ -1,7 +1,7 @@ //! Standard library of models of double theories. use std::rc::Rc; -use ustr::{Ustr, ustr}; +use ustr::{ustr, Ustr}; use crate::dbl::{model::*, theory::*}; use crate::one::Path; @@ -141,6 +141,41 @@ pub fn catalyzed_reaction(th: Rc) -> UstrModalDblModel { model } +/** + */ +pub fn water(th: Rc) -> UstrModalDblModel { + let (state_type, aux_type) = + (ModalObType::new(ustr("State")), ModalObType::new(ustr("Auxiliary"))); + let mut model = UstrModalDblModel::new(th); + let (water, container) = (ustr("water"), ustr("container")); + model.add_ob(water, state_type.clone()); + model.add_ob(container, state_type.clone()); + let (bwater, bcontainer) = (ustr("&water"), ustr("&container")); + model.add_ob(bwater, aux_type.clone()); + model.add_ob(bcontainer, aux_type); + let borrow = ModalMorType::One(ModeApp::new(ustr("borrow"))); + model.add_mor( + ustr("borrow1"), + ModalOb::Generator(water), + ModalOb::Generator(bwater), + borrow.clone(), + ); + model.add_mor( + ustr("borrow2"), + ModalOb::Generator(container), + ModalOb::Generator(bcontainer), + borrow, + ); + let (comparator, comparator_out) = (ustr("comparator"), ustr("comparator_out")); + model.add_mor( + comparator, + ModalOb::List(List::Plain, vec![bwater.into(), bcontainer.into()]).into(), + ModalOb::Generator(comparator_out), + ModalMorType::One(ModeApp::new(ustr("function"))), + ); + model +} + #[cfg(test)] mod tests { use super::super::theories::*; diff --git a/packages/catlog/src/stdlib/theories.rs b/packages/catlog/src/stdlib/theories.rs index 948281e1b..061e29d34 100644 --- a/packages/catlog/src/stdlib/theories.rs +++ b/packages/catlog/src/stdlib/theories.rs @@ -3,7 +3,7 @@ use ustr::ustr; use crate::dbl::theory::*; -use crate::one::{Path, fp_category::UstrFpCategory}; +use crate::one::{fp_category::UstrFpCategory, Path}; /** The empty theory, which has a single model, the empty model. @@ -151,6 +151,11 @@ pub fn th_sym_monoidal_category() -> UstrModalDblTheory { th_list_algebra(List::Symmetric) } +/// +pub fn th_modal_state_aux() -> UstrModalDblTheory { + th_modal_sf(List::Plain) +} + /** The theory of a strict algebra of a list monad. This is a modal double theory, parametric over which variant of the double list @@ -209,6 +214,22 @@ fn th_generalized_multicategory(list: List) -> UstrModalDblTheory { th } +/// +fn th_modal_sf(list: List) -> UstrModalDblTheory { + let m = Modality::List(list); + let mut th: UstrModalDblTheory = Default::default(); + let (state, aux) = (ustr("State"), ustr("Auxiliary")); + th.add_ob_type(state); + th.add_ob_type(aux); + let function = ustr("function"); + th.add_mor_type(function, ModeApp::new(aux).apply(m), ModeApp::new(aux)); + let (borrow, outpos, outneg) = (ustr("borrow"), ustr("out-pos"), ustr("out-neg")); + th.add_mor_type(borrow, ModeApp::new(state), ModeApp::new(aux)); + th.add_mor_type(outpos, ModeApp::new(aux), ModeApp::new(state)); + th.add_mor_type(outneg, ModeApp::new(aux), ModeApp::new(state)); + th +} + #[cfg(test)] mod tests { use super::*; @@ -237,6 +258,7 @@ mod tests { assert!(th_monoidal_category().validate().is_ok()); assert!(th_lax_monoidal_category().validate().is_ok()); assert!(th_multicategory().validate().is_ok()); + assert!(th_modal_state_aux().validate().is_ok()); } #[test] From 44e629c4c7ce1b32d0b602fe7de69c2d800ae5ae Mon Sep 17 00:00:00 2001 From: Matt Cuffaro Date: Sun, 27 Jul 2025 20:00:01 -0700 Subject: [PATCH 2/6] WIP: adding analysis --- packages/catlog-wasm/src/analyses.rs | 4 + packages/catlog-wasm/src/theories.rs | 26 +++ packages/catlog/src/dbl/modal/model.rs | 2 +- packages/catlog/src/dbl/modal/theory.rs | 2 +- .../catlog/src/simulate/ode/linear_ode.rs | 43 ++++ .../src/stdlib/analyses/ode/mass_action.rs | 194 +++++++++++++++++- packages/catlog/src/stdlib/models.rs | 69 +++++-- 7 files changed, 315 insertions(+), 25 deletions(-) diff --git a/packages/catlog-wasm/src/analyses.rs b/packages/catlog-wasm/src/analyses.rs index db5c1dfc8..d7d398a3e 100644 --- a/packages/catlog-wasm/src/analyses.rs +++ b/packages/catlog-wasm/src/analyses.rs @@ -20,3 +20,7 @@ pub struct LinearODEModelData(pub analyses::ode::LinearODEProblemData); #[derive(Serialize, Deserialize, Tsify)] #[tsify(into_wasm_abi, from_wasm_abi)] pub struct MassActionModelData(pub analyses::ode::MassActionProblemData); + +#[derive(Serialize, Deserialize, Tsify)] +#[tsify(into_wasm_abi, from_wasm_abi)] +pub struct AnotherMassActionModelData(pub analyses::ode::AnotherMassActionProblemData); diff --git a/packages/catlog-wasm/src/theories.rs b/packages/catlog-wasm/src/theories.rs index 12835c8ac..c60e82bd4 100644 --- a/packages/catlog-wasm/src/theories.rs +++ b/packages/catlog-wasm/src/theories.rs @@ -337,6 +337,32 @@ impl ThSymMonoidalCategory { } } +/// The theory of state/aux interactions +pub struct ThModalStateAuxCategory(Rc); + +#[wasm_bindgen] +impl ThModalStateAuxCategory { + #[wasm_bindgen(constructor)] + pub fn new() -> Self { + Self(Rc::new(theories::th_modal_sf())) + } + + #[wasm_bindgen] + pub fn theory(&self) -> DblTheory { + DblTheory(self.0.clone().into()) + } + + /// Simulates a mass-action ODE system with additional configurations. + #[wasm_bindgen(js_name = "stateAuxMassAction")] + pub fn state_aux_mass_action( + &self, + model: &DblModel, + data: AnotherMassActionModelData, + ) -> Result<_, String> { + let model: &model::ModalDblModel<_, > = (&model.0).try_into().map_err(|_| "Model should be of a modal theory")?; + Ok(()) + } + #[cfg(test)] mod tests { use super::*; diff --git a/packages/catlog/src/dbl/modal/model.rs b/packages/catlog/src/dbl/modal/model.rs index d43cba7be..fbf7592eb 100644 --- a/packages/catlog/src/dbl/modal/model.rs +++ b/packages/catlog/src/dbl/modal/model.rs @@ -17,7 +17,7 @@ use crate::validate::{self, Validate}; use crate::{one::computad::*, one::*, zero::*}; /// Object in a model of a modal double theory. -#[derive(Clone, Debug, PartialEq, Eq, From)] +#[derive(Clone, Debug, PartialEq, Eq, From, Hash)] pub enum ModalOb { /// Generating object. #[from] diff --git a/packages/catlog/src/dbl/modal/theory.rs b/packages/catlog/src/dbl/modal/theory.rs index 6878c1575..113fdcaaa 100644 --- a/packages/catlog/src/dbl/modal/theory.rs +++ b/packages/catlog/src/dbl/modal/theory.rs @@ -54,7 +54,7 @@ double category of sets admits, besides the [plain](Self::Plain) list double monad, a number of variations decorating the spans of lists with extra combinatorial data. */ -#[derive(Clone, Copy, Debug, PartialEq, Eq)] +#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash)] pub enum List { /// Lists of objects and morphisms (of same length). Plain, diff --git a/packages/catlog/src/simulate/ode/linear_ode.rs b/packages/catlog/src/simulate/ode/linear_ode.rs index 6376e11ff..f786e53dc 100644 --- a/packages/catlog/src/simulate/ode/linear_ode.rs +++ b/packages/catlog/src/simulate/ode/linear_ode.rs @@ -1,6 +1,7 @@ //! Constant-coefficient linear first-order differential equations. use nalgebra::{DMatrix, DVector}; +use std::collections::HashMap; #[cfg(test)] use super::ODEProblem; @@ -30,6 +31,48 @@ impl ODESystem for LinearODESystem { } } +// -------------------------------- + +// TODO move into a context where Id is defined +#[derive(Clone, Debug)] +pub enum Condition { + Geq(usize, usize), +} + +impl Condition { + pub fn call(&self, x: &DVector) -> bool { + match self { + Condition::Geq(left, right) => x.row(*left) >= x.row(*right), + } + } +} + +/** + */ +#[derive(Clone)] +pub struct HybridLinearODESystem { + /// Linear ODEs + pub subsystems: HashMap, +} + +impl HybridLinearODESystem { + // TODO I dont like checking for the first condition to be met. + pub fn get_subsystem(&self, x: &DVector) -> Option<&LinearODESystem> { + self.subsystems + .iter() + .find(|(cond, _)| cond.call(&x.clone())) + .map(|(_, sys)| sys) + } +} + +impl ODESystem for HybridLinearODESystem { + fn vector_field(&self, dx: &mut DVector, x: &DVector, _t: f32) { + if let Some(system) = self.get_subsystem(x) { + system.vector_field(dx, x, _t) + } + } +} + #[cfg(test)] pub(crate) fn create_neg_loops_pos_connector() -> ODEProblem { use nalgebra::{dmatrix, dvector}; diff --git a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs index bc58f28d1..137515022 100644 --- a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs +++ b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs @@ -4,13 +4,15 @@ Such ODEs are based on the *law of mass action* familiar from chemistry and mathematical epidemiology. */ -use std::collections::{BTreeMap, HashMap}; +use std::collections::{BTreeMap, HashMap, HashSet}; use std::fmt::Debug; use std::hash::Hash; +use itertools::Itertools; + use nalgebra::DVector; use num_traits::Zero; -use ustr::{Ustr, ustr}; +use ustr::{ustr, Ustr}; #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; @@ -19,8 +21,8 @@ use tsify::Tsify; use super::ODEAnalysis; use crate::dbl::{ - model::{DiscreteTabModel, FgDblModel, ModalDblModel, MutDblModel, TabEdge}, - theory::{ModalMorType, ModalObType, TabMorType, TabObType}, + model::{DiscreteTabModel, FgDblModel, ModalDblModel, ModalMor, ModalOb, MutDblModel, TabEdge}, + theory::{ModalMorType, ModalObType, ModeApp, TabMorType, TabObType}, }; use crate::one::FgCategory; use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem}; @@ -148,6 +150,7 @@ impl StockFlowMassActionAnalysis { &self, model: &DiscreteTabModel, ) -> PolynomialSystem, u8> { + // associate each flow to its domain. let mut terms: HashMap> = model .mor_generators_with_type(&self.flow_mor_type) .map(|flow| { @@ -223,6 +226,180 @@ fn into_numerical_system( ODEAnalysis::new(problem, ob_index) } +// ------------------------------------------------------------------------- // + +/// Data defining a mass-action ODE problem for a model with function parameters +#[derive(Clone)] +#[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 AnotherMassActionProblemData +where + Id: Eq + Hash, +{ + /// mass action problem data + pub mass: MassActionProblemData, + + /// functions associated to T(aux) --> aux + pub functions: HashMap, +} + +// TODO rename +/** + */ +pub struct PetriNetMassActionFunctionAnalysis { + /// Object type for states + pub state_ob_type: ModalObType, + /// Object type for auxiliary variables + pub aux_ob_type: ModalObType, + /// Morphism type for ... + pub fun_mor_type: ModalMorType, + /// + pub borrow_mor_type: ModalMorType, + /// + pub outpos_mor_type: ModalMorType, + /// + pub outneg_mor_type: ModalMorType, +} + +// TODO +impl Default for PetriNetMassActionFunctionAnalysis { + fn default() -> Self { + Self { + state_ob_type: ModalObType::new(ustr("State")), + aux_ob_type: ModalObType::new(ustr("Auxiliary")), + fun_mor_type: ModalMorType::One(ModeApp::new(ustr("function"))), + borrow_mor_type: ModalMorType::One(ModeApp::new(ustr("borrowing"))), + outpos_mor_type: ModalMorType::One(ModeApp::new(ustr("out-pos"))), + outneg_mor_type: ModalMorType::One(ModeApp::new(ustr("out-neg"))), + } + } +} + +impl PetriNetMassActionFunctionAnalysis { + fn comps(&self, model: &ModalDblModel) -> HashMap> + where + Id: Clone + Debug + Hash + Eq + Ord, + { + let arrows = model.mor_generators_with_type(&self.fun_mor_type); + fn search( + mut component: Vec, + input: ModalOb, + model: &ModalDblModel, + fun_mor_type: &ModalMorType, + ) -> Vec + where + Id: Clone + Debug + Hash + Eq + Ord, + { + let mut out = Vec::new(); + if let Some(arrow) = model + .mor_generators_with_type(fun_mor_type) + .find(|arr| model.mor_generator_cod(arr) == input.clone()) + { + component.push(arrow.clone()); + for el in component.clone() { + out.push(el) + } + for input in vec![model.mor_generator_dom(&arrow)].iter() { + search(component.clone(), input.clone(), model, fun_mor_type); + } + } + out + } + let outpos = model.mor_generators_with_type(&self.outpos_mor_type); + let outneg = model.mor_generators_with_type(&self.outneg_mor_type); + let outpos_dom = outpos.map(|p| model.mor_generator_dom(&p)); + let mut outneg_dom = outneg.map(|p| model.mor_generator_dom(&p)); + + let mediators: Vec<_> = outpos_dom.filter(|el| outneg_dom.contains(el)).collect(); + let programs: HashMap<_, _> = HashMap::from_iter(mediators.into_iter().map(|m| { + (m.clone().unwrap_generator(), search(vec![], m, &model, &self.fun_mor_type)) + })); + programs + } + + pub fn build_system( + &self, + model: &ModalDblModel, + filter_id: Option, + ) -> PolynomialSystem, u8> { + let mut flowneg = HashMap::new(); + let mut flowpos = HashMap::new(); + let outpos = model.mor_generators_with_type(&self.outpos_mor_type); + let outneg = model.mor_generators_with_type(&self.outneg_mor_type); + let outpos_dom: HashSet<_> = HashSet::from_iter::<_>( + outpos + .map(|p| { + let dom = model.mor_generator_dom(&p).unwrap_generator(); + flowpos.insert(dom.clone(), model.mor_generator_cod(&p).unwrap_generator()); + dom + }) + .collect::>(), + ); + let outneg_dom: HashSet<_> = HashSet::from_iter::<_>( + outneg + .map(|p| { + let dom = model.mor_generator_dom(&p).unwrap_generator(); + flowneg.insert(dom.clone(), model.mor_generator_cod(&p).unwrap_generator()); + dom + }) + .collect::>(), + ); + + let mediators: Vec<_> = outpos_dom.intersection(&outneg_dom).cloned().collect(); + let terms: Vec<(Id, Polynomial, u8>)> = mediators + .iter() + .filter(|m| match &filter_id { + Some(k) => *k != **m, + None => true, + }) + .map(|mediator| { + let param = Parameter::generator(mediator.clone()); + let term = Monomial::generator(flowneg[mediator].clone()); + (mediator.clone(), [(param, term)].into_iter().collect()) + }) + .collect(); + + let mut sys = PolynomialSystem::new(); + for ob in model.ob_generators_with_type(&self.state_ob_type) { + sys.add_term(ob, Polynomial::zero()); + } + for (mediator, term) in terms.iter() { + sys.add_term(flowneg[&mediator.clone()].clone(), -term.clone()); + sys.add_term(flowpos[&mediator].clone(), term.clone()); + } + + for (var, component) in sys.clone().components.iter() { + println!("d{var} = {component}") + } + println!("--------------------------"); + + sys + } + + pub fn build_switching_system( + &self, + model: &ModalDblModel, + ) -> HashMap, PolynomialSystem, u8>> { + let programs = self.comps(model); + if programs.is_empty() { + let sys = self.build_system(model, None); + HashMap::from([(None, sys)]) + } else { + HashMap::from_iter(programs.into_iter().map(|(k, v)| { + if v.is_empty() { + (None, self.build_system(model, None)) + } else { + (Some(v[0].clone()), self.build_system(model, Some(k))) + } + })) + } + } +} + #[cfg(test)] mod tests { use expect_test::expect; @@ -255,4 +432,13 @@ mod tests { "#]); expected.assert_eq(&sys.to_string()); } + + #[test] + fn water_bath() { + let th = Rc::new(th_modal_state_aux()); + let model = water(th); + let sys = PetriNetMassActionFunctionAnalysis::default().build_switching_system(&model); + // dbg!(&sys); + assert!(true); + } } diff --git a/packages/catlog/src/stdlib/models.rs b/packages/catlog/src/stdlib/models.rs index e596c3ade..3aa5e4db3 100644 --- a/packages/catlog/src/stdlib/models.rs +++ b/packages/catlog/src/stdlib/models.rs @@ -147,32 +147,56 @@ pub fn water(th: Rc) -> UstrModalDblModel { let (state_type, aux_type) = (ModalObType::new(ustr("State")), ModalObType::new(ustr("Auxiliary"))); let mut model = UstrModalDblModel::new(th); - let (water, container) = (ustr("water"), ustr("container")); - model.add_ob(water, state_type.clone()); - model.add_ob(container, state_type.clone()); - let (bwater, bcontainer) = (ustr("&water"), ustr("&container")); - model.add_ob(bwater, aux_type.clone()); - model.add_ob(bcontainer, aux_type); - let borrow = ModalMorType::One(ModeApp::new(ustr("borrow"))); + let (watershed, lake, sediment) = (ustr("watershed"), ustr("lake"), ustr("sediment")); + model.add_ob(watershed, state_type.clone()); + model.add_ob(lake, state_type.clone()); + model.add_ob(sediment, state_type.clone()); + let blake = ustr("&lake"); + model.add_ob(blake, aux_type.clone()); + let borrow = ModalMorType::One(ModeApp::new(ustr("borrowing"))); model.add_mor( - ustr("borrow1"), - ModalOb::Generator(water), - ModalOb::Generator(bwater), + ustr("borrow"), + ModalOb::Generator(lake.clone()), + ModalOb::Generator(blake.clone()), borrow.clone(), ); - model.add_mor( - ustr("borrow2"), - ModalOb::Generator(container), - ModalOb::Generator(bcontainer), - borrow, - ); - let (comparator, comparator_out) = (ustr("comparator"), ustr("comparator_out")); + let container = ustr("container"); + model.add_ob(container, state_type.clone()); + let (comparator, comparator_out) = (ustr("comparator"), ustr("lake_sediment")); model.add_mor( comparator, - ModalOb::List(List::Plain, vec![bwater.into(), bcontainer.into()]).into(), - ModalOb::Generator(comparator_out), + ModalOb::List(List::Plain, vec![blake.into(), sediment.into()]).into(), + ModalOb::Generator(comparator_out.clone()), ModalMorType::One(ModeApp::new(ustr("function"))), ); + // watershed outflow + let watershed_lake = ustr("watershed_lake"); + model.add_ob(watershed_lake.clone(), aux_type.clone()); + model.add_mor( + ustr("watershed-lake-outflow"), + ModalOb::Generator(watershed_lake.clone()), + ModalOb::Generator(watershed.clone()), + ModalMorType::One(ModeApp::new(ustr("out-neg"))), + ); + model.add_mor( + ustr("watershed-lake-inflow"), + ModalOb::Generator(watershed_lake.clone()), + ModalOb::Generator(lake.clone()), + ModalMorType::One(ModeApp::new(ustr("out-pos"))), + ); + // + model.add_mor( + ustr("lake-sediment-outflow"), + ModalOb::Generator(comparator_out.clone()), + ModalOb::Generator(lake.clone()), + ModalMorType::One(ModeApp::new(ustr("out-neg"))), + ); + model.add_mor( + ustr("lake-sediment-infow"), + ModalOb::Generator(comparator_out.clone()), + ModalOb::Generator(sediment.clone()), + ModalMorType::One(ModeApp::new(ustr("out-pos"))), + ); model } @@ -217,4 +241,11 @@ mod tests { let th = Rc::new(th_sym_monoidal_category()); assert!(catalyzed_reaction(th).validate().is_ok()); } + + #[test] + fn modal_state_aux() { + let th = Rc::new(th_modal_state_aux()); + // dbg!(water(th)); + assert!(true) + } } From 2d7898d18e657accbf404731597d874ec391c731 Mon Sep 17 00:00:00 2001 From: Matt Cuffaro Date: Tue, 29 Jul 2025 18:53:47 -0700 Subject: [PATCH 3/6] WIP: simulating analysis --- packages/catlog-wasm/src/theories.rs | 21 +++- .../catlog/src/simulate/ode/polynomial.rs | 97 ++++++++++++++++++- .../src/stdlib/analyses/ode/mass_action.rs | 58 ++++++++++- .../catlog/src/stdlib/analyses/ode/mod.rs | 14 +++ packages/catlog/src/stdlib/models.rs | 10 +- 5 files changed, 185 insertions(+), 15 deletions(-) diff --git a/packages/catlog-wasm/src/theories.rs b/packages/catlog-wasm/src/theories.rs index c60e82bd4..eaf45d652 100644 --- a/packages/catlog-wasm/src/theories.rs +++ b/packages/catlog-wasm/src/theories.rs @@ -13,7 +13,7 @@ use catlog::dbl::theory; use catlog::one::Path; use catlog::stdlib::{analyses, models, theories, theory_morphisms}; -use super::model_morphism::{MotifsOptions, motifs}; +use super::model_morphism::{motifs, MotifsOptions}; use super::{analyses::*, model::DblModel, theory::DblTheory}; /// The empty or initial theory. @@ -338,13 +338,14 @@ impl ThSymMonoidalCategory { } /// The theory of state/aux interactions +#[wasm_bindgen] pub struct ThModalStateAuxCategory(Rc); #[wasm_bindgen] impl ThModalStateAuxCategory { #[wasm_bindgen(constructor)] pub fn new() -> Self { - Self(Rc::new(theories::th_modal_sf())) + Self(Rc::new(theories::th_modal_state_aux())) } #[wasm_bindgen] @@ -358,10 +359,20 @@ impl ThModalStateAuxCategory { &self, model: &DblModel, data: AnotherMassActionModelData, - ) -> Result<_, String> { - let model: &model::ModalDblModel<_, > = (&model.0).try_into().map_err(|_| "Model should be of a modal theory")?; - Ok(()) + ) -> Result { + let model: &model::ModalDblModel<_, _> = + (&model.0).try_into().map_err(|_| "Model should be of a modal theory")?; + // Ok(()) + Ok(ODEResult( + analyses::ode::PetriNetMassActionFunctionAnalysis::default() + .build_numerical_switching_system(model, data.0) + // TODO Switching System + .solve_with_defaults() + .map_err(|err| format!("{err:?}")) + .into(), + )) } +} #[cfg(test)] mod tests { diff --git a/packages/catlog/src/simulate/ode/polynomial.rs b/packages/catlog/src/simulate/ode/polynomial.rs index cd6bf65e1..6e80fb5cb 100644 --- a/packages/catlog/src/simulate/ode/polynomial.rs +++ b/packages/catlog/src/simulate/ode/polynomial.rs @@ -1,7 +1,9 @@ //! Polynomial differential equations. -use std::collections::BTreeMap; +use std::cmp::{Eq, Ord}; +use std::collections::{BTreeMap, HashMap}; use std::fmt::Display; +use std::hash::Hash; use std::ops::Add; use derivative::Derivative; @@ -127,6 +129,7 @@ where Such a system is ready for use in numerical solvers: the coefficients are floating point numbers and the variables are consecutive integer indices. */ +#[derive(Clone)] pub struct NumericalPolynomialSystem { /// Components of the vector field. pub components: Vec>, @@ -144,6 +147,45 @@ where } } +/** + */ +#[derive(Clone, Derivative)] +#[derivative(Default(bound = ""))] +pub struct NumericalPolynomialSwitchingSystem { + /// Components of a switching system. + pub subsystems: BTreeMap, NumericalPolynomialSystem>, +} + +impl NumericalPolynomialSwitchingSystem { + // TODO + fn get_current(&self, x: DVector) -> &NumericalPolynomialSystem { + &self.subsystems[&None] + } +} + +impl From> + for NumericalPolynomialSwitchingSystem, Exp> +{ + fn from(p: NumericalPolynomialSystem) -> Self { + let subsystems = BTreeMap::from([(None, p)]); + NumericalPolynomialSwitchingSystem { subsystems } + } +} + +impl ODESystem for NumericalPolynomialSwitchingSystem +where + Exp: Clone + Ord, + f32: Pow, + Id: Eq + Hash + Ord, +{ + fn vector_field(&self, dx: &mut DVector, x: &DVector, _t: f32) { + let subsystem = self.get_current(x.clone()); + for i in 0..dx.len() { + dx[i] = subsystem.components[i].eval(|var| x[*var]) + } + } +} + #[cfg(test)] mod tests { use expect_test::expect; @@ -208,4 +250,57 @@ mod tests { "#]]; expected.assert_eq(&textplot_ode_result(&problem, &result)); } + + #[test] + fn water_sim() { + let param = |c: char| Parameter::<_>::generator(c); + let var = |c: char| Polynomial::<_, Parameter<_>, u8>::generator(c); + + let terms = [ + ('l', -var('l') * param('L') + var('w') * param('W')), + ('s', var('l') * param('L')), + ('w', -var('w') * param('W')), + ]; + let starting_sys: PolynomialSystem<_, _, _> = terms.into_iter().collect(); + let starting_sys = starting_sys.extend_scalars(|p| p.eval(|_| 1.0)); + + let terms = [ + ('l', -var('l') * param('L') + var('w') * param('W')), + ('s', var('l') * param('L')), + ('w', -var('w') * param('W')), + ]; + let sys: PolynomialSystem<_, _, _> = terms.into_iter().collect(); + let sys = sys.extend_scalars(|p| p.eval(|_| 1.0)); + + let initial = DVector::from_column_slice(&[1.0, 0.0, 4.0]); + let sys: NumericalPolynomialSwitchingSystem<_, _> = + NumericalPolynomialSwitchingSystem::from(sys.to_numerical()); + let problem = ODEProblem::new(sys, initial).end_time(5.0); + let result = problem.solve_rk4(0.1).unwrap(); + let expected = expect![[r#" + ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⠤⠤⠤⠒⠒⠒⠒⠒⠉⠉⠉⠉⠁ 4.9 + ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⠤⠒⠒⠉⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠤⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠚⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠘⡄⠀⢀⠤⠒⠤⡀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⢣⡔⠁⠀⠀⠀⠈⢦⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⡜⡄⠀⠀⠀⠀⢠⠃⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⡸⠀⢣⠀⠀⠀⢠⠃⠀⠀⠀⠣⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⢠⠃⠀⠘⡄⠀⢠⠃⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⢂⠇⠀⠀⠀⠱⣠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡝⠀⠀⠀⠀⢠⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠅⠀⠀⠀⢠⠃⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⢠⠃⠀⠀⠀⠣⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⡠⠃⠀⠀⠀⠀⠀⠈⠒⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⠒⠤⠤⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⢄⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⠒⠤⠤⠤⠤⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣉⣉⣒⣒⣒⣒⣤⣤⣤⣤⠤⣀⣀⣀⣀⡀ + ⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠉⠉⠉⠉⠉⠁ 0.0 + 0.0 5.0 + "#]]; + expected.assert_eq(&textplot_ode_result(&problem, &result)); + } } diff --git a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs index 137515022..33b971fd2 100644 --- a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs +++ b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs @@ -25,7 +25,9 @@ use crate::dbl::{ theory::{ModalMorType, ModalObType, ModeApp, TabMorType, TabObType}, }; use crate::one::FgCategory; -use crate::simulate::ode::{NumericalPolynomialSystem, ODEProblem, PolynomialSystem}; +use crate::simulate::ode::{ + NumericalPolynomialSwitchingSystem, NumericalPolynomialSystem, ODEProblem, PolynomialSystem, +}; use crate::zero::{alg::Polynomial, rig::Monomial}; /// Data defining a mass-action ODE problem for a model. @@ -295,6 +297,9 @@ impl PetriNetMassActionFunctionAnalysis { Id: Clone + Debug + Hash + Eq + Ord, { let mut out = Vec::new(); + dbg!(model + .mor_generators_with_type(fun_mor_type) + .find(|arr| model.mor_generator_cod(arr) == input.clone())); if let Some(arrow) = model .mor_generators_with_type(fun_mor_type) .find(|arr| model.mor_generator_cod(arr) == input.clone()) @@ -304,9 +309,12 @@ impl PetriNetMassActionFunctionAnalysis { out.push(el) } for input in vec![model.mor_generator_dom(&arrow)].iter() { + dbg!(&input); search(component.clone(), input.clone(), model, fun_mor_type); } } + // TODO return inputs here + dbg!(&input, &out); out } let outpos = model.mor_generators_with_type(&self.outpos_mor_type); @@ -383,13 +391,13 @@ impl PetriNetMassActionFunctionAnalysis { pub fn build_switching_system( &self, model: &ModalDblModel, - ) -> HashMap, PolynomialSystem, u8>> { + ) -> BTreeMap, PolynomialSystem, u8>> { let programs = self.comps(model); if programs.is_empty() { let sys = self.build_system(model, None); - HashMap::from([(None, sys)]) + BTreeMap::from([(None, sys)]) } else { - HashMap::from_iter(programs.into_iter().map(|(k, v)| { + BTreeMap::from_iter(programs.into_iter().map(|(k, v)| { if v.is_empty() { (None, self.build_system(model, None)) } else { @@ -398,6 +406,48 @@ impl PetriNetMassActionFunctionAnalysis { })) } } + + pub fn build_numerical_system( + &self, + model: &ModalDblModel, + data: AnotherMassActionProblemData, + ) -> ODEAnalysis> { + // XXX we forget some data from `data` here + into_numerical_switching_system(self.build_switching_system(model), data.mass) + } +} + +fn into_numerical_switching_system( + switching_system: BTreeMap, PolynomialSystem, u8>>, + data: MassActionProblemData, +) -> ODEAnalysis> { + let subsystems = BTreeMap::from_iter(switching_system.into_iter().map(|(id, sys)| { + let ob_index: BTreeMap<_, _> = + sys.components.keys().cloned().enumerate().map(|(i, x)| (x, i)).collect(); + 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 sys = sys + .extend_scalars(|poly| { + poly.eval(|flow| data.rates.get(flow).copied().unwrap_or_default()) + }) + .to_numerical(); + (id, sys) + })); + let sys = NumericalPolynomialSwitchingSystem { subsystems }; + + let ob_index = BTreeMap::new(); + let x0 = DVector::default(); + // ODE Problem expects polynomial + let problem = ODEProblem::new(sys, x0).end_time(data.duration); + ODEAnalysis { + problem, + variable_index: ob_index, + } } #[cfg(test)] diff --git a/packages/catlog/src/stdlib/analyses/ode/mod.rs b/packages/catlog/src/stdlib/analyses/ode/mod.rs index 63adba85a..708ffb3a4 100644 --- a/packages/catlog/src/stdlib/analyses/ode/mod.rs +++ b/packages/catlog/src/stdlib/analyses/ode/mod.rs @@ -7,6 +7,8 @@ use derivative::Derivative; use derive_more::Constructor; use ode_solvers::dop_shared::IntegrationError; +use nalgebra::DVector; + #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; #[cfg(feature = "serde-wasm")] @@ -67,6 +69,18 @@ impl ODEAnalysis { .collect(), }) } + + /// Solves the ODE construed as a switching system + pub fn solve_switching_with_defaults(self) -> Result, IntegrationError> + where + Id: Eq + Hash, + Sys: ODESystem, + { + // if self.variable_index.is_empty() { + // return Ok(Default::default()); + // } + Ok(Default::default()) + } } pub mod linear_ode; diff --git a/packages/catlog/src/stdlib/models.rs b/packages/catlog/src/stdlib/models.rs index 3aa5e4db3..808e4cc0c 100644 --- a/packages/catlog/src/stdlib/models.rs +++ b/packages/catlog/src/stdlib/models.rs @@ -162,11 +162,11 @@ pub fn water(th: Rc) -> UstrModalDblModel { ); let container = ustr("container"); model.add_ob(container, state_type.clone()); - let (comparator, comparator_out) = (ustr("comparator"), ustr("lake_sediment")); + let (comparator, lake_sediment) = (ustr("comparator"), ustr("lake_sediment")); model.add_mor( comparator, - ModalOb::List(List::Plain, vec![blake.into(), sediment.into()]).into(), - ModalOb::Generator(comparator_out.clone()), + ModalOb::List(List::Plain, vec![blake.into(), container.into()]).into(), + ModalOb::Generator(lake_sediment.clone()), ModalMorType::One(ModeApp::new(ustr("function"))), ); // watershed outflow @@ -187,13 +187,13 @@ pub fn water(th: Rc) -> UstrModalDblModel { // model.add_mor( ustr("lake-sediment-outflow"), - ModalOb::Generator(comparator_out.clone()), + ModalOb::Generator(lake_sediment.clone()), ModalOb::Generator(lake.clone()), ModalMorType::One(ModeApp::new(ustr("out-neg"))), ); model.add_mor( ustr("lake-sediment-infow"), - ModalOb::Generator(comparator_out.clone()), + ModalOb::Generator(lake_sediment.clone()), ModalOb::Generator(sediment.clone()), ModalMorType::One(ModeApp::new(ustr("out-pos"))), ); From b7520f47b0071df970b94f741089ea9788b71be6 Mon Sep 17 00:00:00 2001 From: Matt Cuffaro Date: Thu, 31 Jul 2025 15:35:31 -0700 Subject: [PATCH 4/6] ENH: toposort + folding over subsystems --- packages/catlog/src/one/graph.rs | 16 +- packages/catlog/src/one/graph_algorithms.rs | 174 +++++- .../catlog/src/simulate/ode/polynomial.rs | 240 ++++++-- .../src/stdlib/analyses/ode/mass_action.rs | 558 +++++++++++++----- packages/catlog/src/stdlib/models.rs | 2 +- 5 files changed, 781 insertions(+), 209 deletions(-) diff --git a/packages/catlog/src/one/graph.rs b/packages/catlog/src/one/graph.rs index b162afde9..89d38061e 100644 --- a/packages/catlog/src/one/graph.rs +++ b/packages/catlog/src/one/graph.rs @@ -194,20 +194,20 @@ of [`FinGraph`]. */ pub trait ColumnarFinGraph: ColumnarGraph< - Vertices: FinSet, - Edges: FinSet, - Src: Column, - Tgt: Column, - > + Vertices: FinSet, + Edges: FinSet, + Src: Column, + Tgt: Column, +> { } /// A columnar graph with mutable columns. pub trait MutColumnarGraph: ColumnarGraph< - Src: MutMapping, - Tgt: MutMapping, - > + Src: MutMapping, + Tgt: MutMapping, +> { /// Variant of [`src_map`](ColumnarGraph::src_map) that returns a mutable /// reference. diff --git a/packages/catlog/src/one/graph_algorithms.rs b/packages/catlog/src/one/graph_algorithms.rs index 6395b16c1..7fd15fc64 100644 --- a/packages/catlog/src/one/graph_algorithms.rs +++ b/packages/catlog/src/one/graph_algorithms.rs @@ -1,6 +1,6 @@ //! Algorithms on graphs. -use std::collections::{HashSet, VecDeque}; +use std::collections::{HashMap, HashSet, VecDeque}; use std::hash::Hash; use super::graph::*; @@ -168,6 +168,136 @@ where result } +fn out_neighbors(graph: &G, v: &G::V) -> impl Iterator +where + G: FinGraph, + G::V: Hash, +{ + graph.out_edges(v).map(|e| graph.tgt(&e)) +} + +fn in_neighbors(graph: &G, v: &G::V) -> impl Iterator +where + G: FinGraph, + G::V: Hash, +{ + graph.in_edges(v).map(|e| graph.src(&e)) +} + +#[derive(Clone, Debug)] +struct VisitMap { + visited: Vec, +} + +impl VisitMap { + fn visit(&mut self, idx: usize) -> bool { + let previous = self.visited[idx]; + self.visited[idx] = true; + !previous + } + + fn is_visited(&self, idx: usize) -> bool { + self.visited[idx] + } +} + +#[derive(Clone, Debug)] +pub struct Dfs +where + G: FinGraph, + G::V: Hash, +{ + // + pub stack: Vec, + // + pub discovered: VisitMap, +} + +// TODO +// 1. replace String with Cycle error +/** Computes a topological sorting for a given graph. + +This algorithm was borrowed from `petgraph`. + */ +pub fn toposort(graph: &G) -> Result, String> +where + G: FinGraph, + G::V: Hash + std::fmt::Debug, +{ + // XXX dont clone + let n = graph.vertices().collect::>().len(); + let mut discovered = VisitMap { + visited: vec![false; n.clone()], + }; + let mut finished = VisitMap { + visited: vec![false; n.clone()], + }; + let mut finish_stack: Vec = Vec::new(); + let mut stack = Vec::new(); + + // we shouldn't need to do this + let gmap: HashMap<_, _> = HashMap::from_iter(graph.vertices().enumerate().map(|(k, v)| (v, k))); + + for (idx, v) in graph.vertices().enumerate() { + if discovered.is_visited(idx) { + continue; + } + stack.push(v); + while let Some(nx) = stack.clone().last() { + if discovered.visit(gmap[&nx]) { + for succ in out_neighbors(graph, &nx) { + if succ == *nx { + return Err("self cycle".to_owned()); + } + if !discovered.is_visited(gmap[&succ]) { + stack.push(succ); + } + } + } else { + stack.pop(); + if finished.visit(gmap[&nx]) { + finish_stack.push(nx.clone()); + } + } + } + } + finish_stack.reverse(); + + // dfs.reset(g); + let mut discovered = VisitMap { + visited: vec![false; n.clone()], + }; + for i in &finish_stack { + // dfs.move_to(i); + stack.clear(); + stack.push(i.clone()); + // + let mut cycle = false; + while let Some(j) = { + let mut out = None; + while let Some(node) = stack.pop() { + if discovered.visit(gmap[&node]) { + for succ in in_neighbors(graph, &node) { + if !discovered.is_visited(gmap[&succ]) { + stack.push(succ); + } + } + out = Some(node); + break; + } + } + out + } { + if cycle { + return Err(format!("cycle detected involving node {:#?}", j).to_owned()); + } + cycle = true; + } + } + + Ok(finish_stack) +} + #[cfg(test)] mod tests { use super::GraphElem::*; @@ -231,4 +361,46 @@ mod tests { let g = SkelGraph::cycle(1); assert_eq!(spec_order_all(&g), vec![Vertex(0), Edge(0)]); } + + #[test] + fn toposorting() { + let mut g = SkelGraph::path(5); + assert_eq!(toposort(&g), Ok(vec![0, 1, 2, 3, 4])); + + let mut g = SkelGraph::path(3); + g.add_vertices(1); + let _ = g.add_edge(2, 3); + let _ = g.add_edge(3, 0); + assert_eq!(toposort(&g), Err("cycle detected involving node 3".to_owned())); + + let g = SkelGraph::triangle(); + assert_eq!(toposort(&g), Ok(vec![0, 1, 2])); + + let mut g = SkelGraph::path(4); + g.add_vertices(2); + let _ = g.add_edge(1, 4); + let _ = g.add_edge(4, 3); + let _ = g.add_edge(5, 2); + assert_eq!(toposort(&g), Ok(vec![5, 0, 1, 2, 4, 3])); + + let mut g: HashGraph = Default::default(); + g.add_vertices(vec![0, 1, 2, 3, 4, 5]); + g.add_edge("0-1", 0, 1); + g.add_edge("1-2", 1, 2); + g.add_edge("2-3", 2, 3); + g.add_edge("1-4", 1, 4); + g.add_edge("4-3", 4, 3); + g.add_edge("5-2", 5, 2); + // TODO non-deterministic + // assert_eq!(toposort(&g), Ok(vec![5, 0, 1, 2, 4, 3])); + } + + #[test] + #[test] + fn neighbors() { + let g = SkelGraph::triangle(); + let out_neighbors = &g.out_edges(&1).map(|e| g.tgt(&e)).collect::>(); + // dbg!(&out_neighbors); + assert!(true); + } } diff --git a/packages/catlog/src/simulate/ode/polynomial.rs b/packages/catlog/src/simulate/ode/polynomial.rs index 6e80fb5cb..e771b0258 100644 --- a/packages/catlog/src/simulate/ode/polynomial.rs +++ b/packages/catlog/src/simulate/ode/polynomial.rs @@ -2,14 +2,18 @@ use std::cmp::{Eq, Ord}; use std::collections::{BTreeMap, HashMap}; -use std::fmt::Display; +use std::fmt::{Debug, Display}; use std::hash::Hash; use std::ops::Add; +use itertools::{EitherOrBoth::*, Itertools}; + use derivative::Derivative; use nalgebra::DVector; use num_traits::{One, Pow, Zero}; +use crate::stdlib::analyses::ode::{ComputeGraph, EligibleFunctions, Output}; + #[cfg(test)] use super::ODEProblem; use super::ODESystem; @@ -135,6 +139,26 @@ pub struct NumericalPolynomialSystem { pub components: Vec>, } +impl Add for NumericalPolynomialSystem +where + Exp: Ord, +{ + type Output = Self; + + fn add(self, rhs: Self) -> Self::Output { + let components = self + .components + .into_iter() + .zip_longest(rhs.components.into_iter()) + .map(|pair| match pair { + Both(l, r) => l + r, + Left(x) | Right(x) => x, + }) + .collect::>(); + Self { components } + } +} + impl ODESystem for NumericalPolynomialSystem where Exp: Clone + Ord, @@ -151,37 +175,148 @@ where */ #[derive(Clone, Derivative)] #[derivative(Default(bound = ""))] -pub struct NumericalPolynomialSwitchingSystem { +pub struct NumericalPolynomialSwitchingSystem +where + Id: Hash + Clone + Eq + Debug, + Exp: Clone, +{ + /// + pub ob_index: BTreeMap, + + // TODO public? + /// Null model + pub null_model: Option>, + /// Components of a switching system. - pub subsystems: BTreeMap, NumericalPolynomialSystem>, + pub subsystems: Vec<(ComputeGraph, NumericalPolynomialSystem)>, + + /// Analysis data + pub functions: HashMap, } -impl NumericalPolynomialSwitchingSystem { - // TODO - fn get_current(&self, x: DVector) -> &NumericalPolynomialSystem { - &self.subsystems[&None] +impl NumericalPolynomialSwitchingSystem +where + Id: Hash + Clone + Eq + Debug + Ord, + Exp: Ord + Clone, +{ + // TODO accumulate over first arguments + fn get_current(self, x: DVector) -> NumericalPolynomialSystem { + self.subsystems + .into_iter() + .filter(|(graph, _)| graph.compute(&x, self.functions.clone(), self.ob_index.clone())) + .fold(self.null_model.unwrap(), |acc, (_, sys)| acc + sys) + // if x[0] >= 4.5f32 { + // let (g, sys0) = self.subsystems[0].clone(); + // let (_, sys1) = self.subsystems[1].clone(); + // let _ = &g.compute(&x, self.functions); + // sys0 + sys1 + // } else { + // let (graph, sys) = self.subsystems[1].clone(); + // sys + // } } } -impl From> - for NumericalPolynomialSwitchingSystem, Exp> +// impl From> +// for NumericalPolynomialSwitchingSystem +// where +// Exp: Clone, +// { +// fn from(nps: NumericalPolynomialSystem) -> Self { +// let subsystems = vec![(ComputeGraph::::new(), nps)]; +// let functions = HashMap::from([(ustr::ustr("comparator"), EligibleFunctions::Geq())]); +// NumericalPolynomialSwitchingSystem { +// subsystems, +// functions, +// } +// } +// } + +impl ComputeGraph +where + Id: Eq + Hash + Debug + Clone + Ord, { - fn from(p: NumericalPolynomialSystem) -> Self { - let subsystems = BTreeMap::from([(None, p)]); - NumericalPolynomialSwitchingSystem { subsystems } + // + fn compute( + &self, + x: &DVector, + functions: HashMap, + ob_index: BTreeMap, + ) -> bool { + // binding name + let mut bindings: Vec<(_, Output)> = Default::default(); + // dbg!(&self.borrows, &self.obs); + for var in self.toposort.iter() { + if self.obs.contains(&var) { + let index = ob_index[var]; + bindings.push((var, Output::Float(x[index]))); + } else if self.borrows.keys().contains(&var) { + let index = ob_index[&self.borrows[var]]; + bindings.push((var, Output::Float(x[index]))); + } else if let Some(val) = self.auxes.get(&var) { + bindings.push((var, Output::Float(*val))); + } else if self.funcs.keys().contains(&var) { + let args = &self.funcs[&var] + .clone() + .into_iter() + .map(|arg| self.fetch(x, arg, ob_index.clone())) + .collect::>(); + let out = functions[&var].clone().eval(args); + bindings.push((var, out)); + } + } + if let Some((_, value)) = bindings.last() { + match value { + Output::Bool(x) => *x, + _ => false, + } + } else { + true + } + } + + // TODO parameterize + fn fetch(&self, x: &DVector, arg: Id, ob_index: BTreeMap) -> f32 { + if self.obs.contains(&arg) { + x[ob_index[&arg]] + } else if self.borrows.keys().contains(&arg) { + x[ob_index[&self.borrows[&arg]]] + } else { + self.auxes[&arg] + } } } +// impl From, PolynomialSystem)>> +// for NumericalPolynomialSwitchingSystem +// where +// Var: Ord, +// Exp: Ord, +// Id: Hash + Clone + Eq + Debug, +// { +// fn from(subsystems: Vec<(ComputeGraph, PolynomialSystem)>) -> Self { +// let subsystems = subsystems +// .into_iter() +// .map(|(id, poly)| (id, poly.to_numerical())) +// .collect::>(); +// NumericalPolynomialSwitchingSystem { subsystems } +// } +// } + impl ODESystem for NumericalPolynomialSwitchingSystem where - Exp: Clone + Ord, + Exp: Clone + Ord + Display + num_traits::One, f32: Pow, - Id: Eq + Hash + Ord, + Id: Eq + Hash + Ord + Clone + Debug, { fn vector_field(&self, dx: &mut DVector, x: &DVector, _t: f32) { - let subsystem = self.get_current(x.clone()); + let subsystem = self.clone().get_current(x.clone()); for i in 0..dx.len() { - dx[i] = subsystem.components[i].eval(|var| x[*var]) + // TODO do we still need to .get() ? + dx[i] = match subsystem.components.get(i) { + Some(p) => p.eval(|var| x[*var]), + None => 0f32, + } } } } @@ -252,55 +387,46 @@ mod tests { } #[test] - fn water_sim() { - let param = |c: char| Parameter::<_>::generator(c); - let var = |c: char| Polynomial::<_, Parameter<_>, u8>::generator(c); + fn water_simulation() { + let th = std::rc::Rc::new(crate::stdlib::th_modal_state_aux()); + let model = crate::stdlib::water(th); + let system = crate::stdlib::analyses::ode::PetriNetMassActionFunctionAnalysis::default() + .build_switching_system(&model); - let terms = [ - ('l', -var('l') * param('L') + var('w') * param('W')), - ('s', var('l') * param('L')), - ('w', -var('w') * param('W')), - ]; - let starting_sys: PolynomialSystem<_, _, _> = terms.into_iter().collect(); - let starting_sys = starting_sys.extend_scalars(|p| p.eval(|_| 1.0)); + // TODO this should be obtained from the petri net analysis - let terms = [ - ('l', -var('l') * param('L') + var('w') * param('W')), - ('s', var('l') * param('L')), - ('w', -var('w') * param('W')), - ]; - let sys: PolynomialSystem<_, _, _> = terms.into_iter().collect(); - let sys = sys.extend_scalars(|p| p.eval(|_| 1.0)); + let functions = HashMap::from([(ustr::ustr("comparator"), EligibleFunctions::Geq())]); + let sys = system.to_numerical(functions); - let initial = DVector::from_column_slice(&[1.0, 0.0, 4.0]); - let sys: NumericalPolynomialSwitchingSystem<_, _> = - NumericalPolynomialSwitchingSystem::from(sys.to_numerical()); + // sediment, watershed, lake + let initial = DVector::from_column_slice(&[1.0, 2.0, 4.0]); let problem = ODEProblem::new(sys, initial).end_time(5.0); let result = problem.solve_rk4(0.1).unwrap(); let expected = expect![[r#" - ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⠤⠤⠤⠒⠒⠒⠒⠒⠉⠉⠉⠉⠁ 4.9 - ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⠤⠒⠒⠉⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠤⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠒⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠚⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡁⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠄⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠂⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡁⠀⠘⡄⠀⢀⠤⠒⠤⡀⠀⢠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠄⠀⠀⢣⡔⠁⠀⠀⠀⠈⢦⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠂⠀⠀⡜⡄⠀⠀⠀⠀⢠⠃⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡁⠀⡸⠀⢣⠀⠀⠀⢠⠃⠀⠀⠀⠣⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠄⢠⠃⠀⠘⡄⠀⢠⠃⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⢂⠇⠀⠀⠀⠱⣠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡝⠀⠀⠀⠀⢠⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠅⠀⠀⠀⢠⠃⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⠂⠀⠀⢠⠃⠀⠀⠀⠣⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⡁⠀⡠⠃⠀⠀⠀⠀⠀⠈⠒⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⠒⠤⠤⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ - ⢄⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⠒⠤⠤⠤⠤⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⣉⣉⣒⣒⣒⣒⣤⣤⣤⣤⠤⣀⣀⣀⣀⡀ - ⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠁⠈⠀⠉⠉⠉⠉⠉⠁ 0.0 + ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠒⢣⠀⣀⠤⠤⠒⠒⠉⠒⠒⠒⠉⠉⠒⠒⠒⠒⠉⠉⠉⠉⠉⠒⠒⠒⠒⠒⠒⠂ 4.5 + ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠒⠉⠀⠀⠀⠀⠀⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠤⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⣱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⢇⠀⠀⠀⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠘⡄⠀⠀⠀⠀⠀⡔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠱⡀⠀⠀⠀⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⢣⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⠀⢇⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⣀⣀⣀⡀ + ⡁⠀⠀⠀⡎⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⣀⣀⣀⣀⠤⠤⠤⠤⠤⠒⠒⠒⠒⠒⠒⠒⠒⠊⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⡜⠀⠈⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠒⠒⢲⠓⠒⠒⠚⢖⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⢀⠇⠀⠀⠀⠀⠀⠣⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⡜⠀⠀⠀⠀⠀⠀⠀⠑⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⢲⠁⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠒⠤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠒⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⠒⠤⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠒⠒⠒⠤⠤⠤⠤⣀⣀⣀⣀⣀⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠉⠉⠉⠉⠁ 0.0 0.0 5.0 "#]]; + // println!("{}", &textplot_ode_result(&problem, &result)); expected.assert_eq(&textplot_ode_result(&problem, &result)); } } diff --git a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs index 33b971fd2..0b5520e86 100644 --- a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs +++ b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs @@ -5,7 +5,7 @@ mathematical epidemiology. */ use std::collections::{BTreeMap, HashMap, HashSet}; -use std::fmt::Debug; +use std::fmt::{Debug, Display}; use std::hash::Hash; use itertools::Itertools; @@ -21,10 +21,10 @@ use tsify::Tsify; use super::ODEAnalysis; use crate::dbl::{ - model::{DiscreteTabModel, FgDblModel, ModalDblModel, ModalMor, ModalOb, MutDblModel, TabEdge}, + model::{DiscreteTabModel, FgDblModel, ModalDblModel, ModalOb, MutDblModel, TabEdge}, theory::{ModalMorType, ModalObType, ModeApp, TabMorType, TabObType}, }; -use crate::one::FgCategory; +use crate::one::{FgCategory, FinGraph, HashGraph}; use crate::simulate::ode::{ NumericalPolynomialSwitchingSystem, NumericalPolynomialSystem, ODEProblem, PolynomialSystem, }; @@ -230,6 +230,68 @@ fn into_numerical_system( // ------------------------------------------------------------------------- // +// trait Evaluable { +// type Output; +// fn eval(&self) -> Self::Output; +// } + +// pub struct GeqFunction(T, T); +// pub struct IdentityFunction(T); + +// impl Evaluable for GeqFunction { +// type Output = bool; + +// fn eval(&self) -> Self::Output { +// self.0 >= self.1 +// } +// } + +// impl Evaluable for IdentityFunction +// where +// T: Clone, +// { +// type Output = T; + +// fn eval(&self) -> Self::Output { +// self.0.clone() +// } +// } + +// pub enum EligibleFunctions { +// Geq(Y), +// Identity(Y), +// } + +// impl EligibleFunctions where T: PartialOrd { +// pub fn eval(&self, args: Vec) -> Y { + +#[derive(Debug)] +pub enum Output { + Bool(bool), + Float(f32), +} + +#[derive(Debug, Clone)] +#[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 enum EligibleFunctions { + Geq(), + Identity(), +} + +impl EligibleFunctions { + pub fn eval(self, args: &Vec) -> Output { + match self { + EligibleFunctions::Geq() => Output::Bool(args[0] >= args[1]), + EligibleFunctions::Identity() => Output::Float(args[0]), + } + } +} + /// Data defining a mass-action ODE problem for a model with function parameters #[derive(Clone)] #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] @@ -240,13 +302,194 @@ fn into_numerical_system( )] pub struct AnotherMassActionProblemData where - Id: Eq + Hash, + Id: Eq + Hash + Clone + Debug, { /// mass action problem data pub mass: MassActionProblemData, /// functions associated to T(aux) --> aux - pub functions: HashMap, + pub functions: HashMap, +} + +#[derive(Debug, Clone)] +pub struct ComputeGraph +where + Id: Clone + Eq + Hash + Debug, +{ + /// compute graph + pub graph: HashGraph, // TODO parameterize edge type + + pub obs: Vec, + + pub funcs: HashMap>, + + pub auxes: HashMap, + + pub borrows: HashMap, + + // /// in-ports + + // /// out-ports + // inp + /// + pub toposort: Vec, +} + +impl ComputeGraph { + pub fn new() -> Self { + Self { + graph: HashGraph::::default(), + obs: Vec::::new(), + auxes: HashMap::::new(), + funcs: HashMap::>::new(), // store inputs + borrows: HashMap::::new(), + toposort: vec![], + } + } + + // TODO need to verify that arrow is an arrow. Otherwise, return an identity function + // XXX update this + /// This extracts a compute graph associated to an arrow in a model. + pub fn complete(model: &ModalDblModel, arrow: Id) -> Self { + let mut cg = ComputeGraph::::new(); + match model.mor_generator_dom(&arrow.clone()) { + ModalOb::Generator(id) => { + cg.add_vertex(arrow.clone()); + cg.add_vertex(id.clone()); + cg.connect(id, arrow.clone()); + } + ModalOb::List(_, xs) => { + cg.add_vertex(arrow.clone()); + for x in xs { + // TODO `x` may not be Generator + cg.add_vertex(x.clone().unwrap_generator()); + cg.connect(x.clone().unwrap_generator(), arrow.clone()); + } + } + _ => todo!(), + } + match model.mor_generator_cod(&arrow) { + ModalOb::Generator(id) => { + cg.add_vertex(arrow.clone()); + cg.add_vertex(id.clone()); + cg.connect(arrow.clone(), id.clone()); + } + ModalOb::List(_, xs) => { + cg.add_vertex(arrow.clone()); + for x in xs { + // TODO `x` may not be Generator + cg.add_vertex(x.clone().unwrap_generator()); + cg.connect(arrow.clone(), x.unwrap_generator()); + } + } + _ => todo!(), + } + + cg.toposort = crate::one::graph_algorithms::toposort(&cg.graph).expect("!"); + cg + } + + pub fn add_vertex(&mut self, v: Id) -> bool { + self.graph.add_vertex(v) + } + + pub fn add_edge(&mut self, e: String, dom: Id, cod: Id) -> bool { + self.graph.add_edge(e, dom, cod) + } + + // TODO check that the name exists? + /// Connects a src and tgt by an edge whose name is generated automatically. + pub fn connect(&mut self, dom: Id, cod: Id) -> bool { + self.graph.add_edge(format!("{:?}=>{:?}", dom.clone(), cod.clone()), dom, cod) + } +} + +#[derive(Clone)] +pub struct SwitchingSystem +where + Id: Eq + Clone + Hash + Debug, +{ + pub null_model: PolynomialSystem, u8>, + pub subsystems: Vec<(ComputeGraph, PolynomialSystem, u8>)>, + pub ob_index: BTreeMap, +} + +impl SwitchingSystem +where + Id: Hash + Clone + Eq + Ord + Debug, +{ + pub fn to_numerical( + &self, + functions: HashMap, + ) -> NumericalPolynomialSwitchingSystem { + let subsystems: Vec<(ComputeGraph, _)> = self + .subsystems + .clone() + .into_iter() + .map(|(id, poly)| (id, poly.extend_scalars(|p| p.eval(|_| 1.0)).to_numerical())) + .collect(); + // TODO option is used here to avoid implementing Default. + let null_model = + Some(self.null_model.clone().extend_scalars(|p| p.eval(|_| 1.0)).to_numerical()); + NumericalPolynomialSwitchingSystem { + ob_index: self.ob_index.clone(), + functions, + null_model, + subsystems, + } + } +} + +/// Convenience struct +#[derive(Clone, Debug)] +pub struct PetriNetSystemData +where + Id: Eq + Clone + Hash + Debug, +{ + pub flowneg: HashMap, + pub flowpos: HashMap, + pub outpos_dom: HashSet, + pub outneg_dom: HashSet, + pub mediators: Vec, +} + +impl PetriNetSystemData +where + Id: Hash + Eq + Clone + Debug, +{ + fn make(model: &ModalDblModel, pna: &PetriNetMassActionFunctionAnalysis) -> Self { + let mut flowneg = HashMap::new(); + let mut flowpos = HashMap::new(); + let outpos_dom: HashSet<_> = HashSet::from_iter::<_>( + model + .mor_generators_with_type(&pna.outpos_mor_type) + .map(|p| { + let dom = model.mor_generator_dom(&p).unwrap_generator(); + flowpos.insert(dom.clone(), model.mor_generator_cod(&p).unwrap_generator()); + dom + }) + .collect::>(), + ); + let outneg_dom: HashSet<_> = HashSet::from_iter::<_>( + model + .mor_generators_with_type(&pna.outneg_mor_type) + .map(|p| { + let dom = model.mor_generator_dom(&p).unwrap_generator(); + flowneg.insert(dom.clone(), model.mor_generator_cod(&p).unwrap_generator()); + dom + }) + .collect::>(), + ); + let mediators: Vec<_> = outpos_dom.intersection(&outneg_dom).cloned().collect(); + + Self { + flowneg, + flowpos, + outpos_dom, + outneg_dom, + mediators, + } + } } // TODO rename @@ -257,17 +500,16 @@ pub struct PetriNetMassActionFunctionAnalysis { pub state_ob_type: ModalObType, /// Object type for auxiliary variables pub aux_ob_type: ModalObType, - /// Morphism type for ... + /// Morphism type for functions between auxiliary variables pub fun_mor_type: ModalMorType, - /// + /// Morphism type for "borrowing," or maps from states to auxiliary variables pub borrow_mor_type: ModalMorType, - /// + /// Morphism type for positively-signed flows pub outpos_mor_type: ModalMorType, - /// + /// Morphism type for negatively-signed flows pub outneg_mor_type: ModalMorType, } -// TODO impl Default for PetriNetMassActionFunctionAnalysis { fn default() -> Self { Self { @@ -282,91 +524,86 @@ impl Default for PetriNetMassActionFunctionAnalysis { } impl PetriNetMassActionFunctionAnalysis { - fn comps(&self, model: &ModalDblModel) -> HashMap> - where - Id: Clone + Debug + Hash + Eq + Ord, - { - let arrows = model.mor_generators_with_type(&self.fun_mor_type); - fn search( - mut component: Vec, - input: ModalOb, - model: &ModalDblModel, - fun_mor_type: &ModalMorType, - ) -> Vec - where - Id: Clone + Debug + Hash + Eq + Ord, - { - let mut out = Vec::new(); - dbg!(model - .mor_generators_with_type(fun_mor_type) - .find(|arr| model.mor_generator_cod(arr) == input.clone())); - if let Some(arrow) = model - .mor_generators_with_type(fun_mor_type) - .find(|arr| model.mor_generator_cod(arr) == input.clone()) - { - component.push(arrow.clone()); - for el in component.clone() { - out.push(el) + fn build_graph( + &self, + model: &ModalDblModel, + ) -> ComputeGraph { + let mut cg = ComputeGraph::::new(); + cg.obs = model.ob_generators_with_type(&self.state_ob_type).collect::>(); + // TODO these should just be auxes entering the graph + cg.auxes = HashMap::from_iter( + model.ob_generators_with_type(&self.aux_ob_type).map(|aux| (aux, 4.5f32)), + ); + cg.borrows = + HashMap::from_iter(model.mor_generators_with_type(&self.borrow_mor_type).map(|f| { + ( + model.mor_generator_cod(&f.clone()).unwrap_generator(), + model.mor_generator_dom(&f.clone()).unwrap_generator(), + ) + })); + let arrows = model.mor_generators_with_type(&self.fun_mor_type).collect::>(); + for arrow in arrows { + let mut inputs: Vec = vec![]; + match model.mor_generator_dom(&arrow) { + ModalOb::Generator(id) => { + inputs.push(id.clone()); + cg.add_vertex(arrow.clone()); + cg.add_vertex(id.clone()); + cg.connect(id, arrow.clone()); + } + ModalOb::List(_, xs) => { + cg.add_vertex(arrow.clone()); + for x in xs { + inputs.push(x.clone().unwrap_generator()); + // TODO `x` may not be Generator + cg.add_vertex(x.clone().unwrap_generator()); + cg.connect(x.clone().unwrap_generator(), arrow.clone()); + } + } + _ => todo!(), + } + cg.funcs.insert(arrow.clone(), inputs); + // + match model.mor_generator_cod(&arrow) { + ModalOb::Generator(id) => { + cg.add_vertex(arrow.clone()); + cg.add_vertex(id.clone()); + cg.connect(arrow.clone(), id.clone()); } - for input in vec![model.mor_generator_dom(&arrow)].iter() { - dbg!(&input); - search(component.clone(), input.clone(), model, fun_mor_type); + ModalOb::List(_, xs) => { + cg.add_vertex(arrow.clone()); + for x in xs { + // TODO `x` may not be Generator + cg.add_vertex(x.clone().unwrap_generator()); + cg.connect(arrow.clone(), x.unwrap_generator()); + } } + _ => todo!(), } - // TODO return inputs here - dbg!(&input, &out); - out } - let outpos = model.mor_generators_with_type(&self.outpos_mor_type); - let outneg = model.mor_generators_with_type(&self.outneg_mor_type); - let outpos_dom = outpos.map(|p| model.mor_generator_dom(&p)); - let mut outneg_dom = outneg.map(|p| model.mor_generator_dom(&p)); - - let mediators: Vec<_> = outpos_dom.filter(|el| outneg_dom.contains(el)).collect(); - let programs: HashMap<_, _> = HashMap::from_iter(mediators.into_iter().map(|m| { - (m.clone().unwrap_generator(), search(vec![], m, &model, &self.fun_mor_type)) - })); - programs + + cg.toposort = crate::one::graph_algorithms::toposort(&cg.graph).expect("!"); + cg } - pub fn build_system( + pub fn build_system( &self, model: &ModalDblModel, - filter_id: Option, - ) -> PolynomialSystem, u8> { - let mut flowneg = HashMap::new(); - let mut flowpos = HashMap::new(); - let outpos = model.mor_generators_with_type(&self.outpos_mor_type); - let outneg = model.mor_generators_with_type(&self.outneg_mor_type); - let outpos_dom: HashSet<_> = HashSet::from_iter::<_>( - outpos - .map(|p| { - let dom = model.mor_generator_dom(&p).unwrap_generator(); - flowpos.insert(dom.clone(), model.mor_generator_cod(&p).unwrap_generator()); - dom - }) - .collect::>(), - ); - let outneg_dom: HashSet<_> = HashSet::from_iter::<_>( - outneg - .map(|p| { - let dom = model.mor_generator_dom(&p).unwrap_generator(); - flowneg.insert(dom.clone(), model.mor_generator_cod(&p).unwrap_generator()); - dom - }) - .collect::>(), - ); - - let mediators: Vec<_> = outpos_dom.intersection(&outneg_dom).cloned().collect(); - let terms: Vec<(Id, Polynomial, u8>)> = mediators + model_data: PetriNetSystemData, // TODO rename + filter_ids: Vec, + ) -> PolynomialSystem, u8> + where + Id: Eq + Clone + Hash + Debug + Ord + Display, + { + // TODO what are these + let terms: Vec<(Id, Polynomial, u8>)> = model_data + .mediators .iter() - .filter(|m| match &filter_id { - Some(k) => *k != **m, - None => true, - }) + // keep the mediator if it is in the list + .filter(|m| filter_ids.clone().contains(m)) .map(|mediator| { let param = Parameter::generator(mediator.clone()); - let term = Monomial::generator(flowneg[mediator].clone()); + let term = Monomial::generator(model_data.flowneg[mediator].clone()); (mediator.clone(), [(param, term)].into_iter().collect()) }) .collect(); @@ -375,9 +612,10 @@ impl PetriNetMassActionFunctionAnalysis { for ob in model.ob_generators_with_type(&self.state_ob_type) { sys.add_term(ob, Polynomial::zero()); } + // TODO for (mediator, term) in terms.iter() { - sys.add_term(flowneg[&mediator.clone()].clone(), -term.clone()); - sys.add_term(flowpos[&mediator].clone(), term.clone()); + sys.add_term(model_data.flowneg[&mediator.clone()].clone(), -term.clone()); + sys.add_term(model_data.flowpos[&mediator].clone(), term.clone()); } for (var, component) in sys.clone().components.iter() { @@ -388,68 +626,95 @@ impl PetriNetMassActionFunctionAnalysis { sys } - pub fn build_switching_system( - &self, - model: &ModalDblModel, - ) -> BTreeMap, PolynomialSystem, u8>> { - let programs = self.comps(model); - if programs.is_empty() { - let sys = self.build_system(model, None); - BTreeMap::from([(None, sys)]) - } else { - BTreeMap::from_iter(programs.into_iter().map(|(k, v)| { - if v.is_empty() { - (None, self.build_system(model, None)) - } else { - (Some(v[0].clone()), self.build_system(model, Some(k))) - } - })) - } - } + pub fn build_switching_system(&self, model: &ModalDblModel) -> SwitchingSystem + where + Id: Clone + Eq + Hash + Debug + Ord + Display, + { + let modeldata = PetriNetSystemData::make(model.into(), &self); - pub fn build_numerical_system( - &self, - model: &ModalDblModel, - data: AnotherMassActionProblemData, - ) -> ODEAnalysis> { - // XXX we forget some data from `data` here - into_numerical_switching_system(self.build_switching_system(model), data.mass) - } -} + // TODO ensure they are connected components + let programs = self.build_graph(model); -fn into_numerical_switching_system( - switching_system: BTreeMap, PolynomialSystem, u8>>, - data: MassActionProblemData, -) -> ODEAnalysis> { - let subsystems = BTreeMap::from_iter(switching_system.into_iter().map(|(id, sys)| { - let ob_index: BTreeMap<_, _> = - sys.components.keys().cloned().enumerate().map(|(i, x)| (x, i)).collect(); - 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 sys = sys - .extend_scalars(|poly| { - poly.eval(|flow| data.rates.get(flow).copied().unwrap_or_default()) + // build subsystem with no programs. + // This means we exclude all flows are are mediated by + // any program + let affected_mediators = modeldata + .clone() + .mediators + .into_iter() + .filter(|m| programs.graph.vertices().contains(m)); + let init_model = self.build_system( + &model, + modeldata.clone(), + modeldata + .clone() + .mediators + .into_iter() + .filter(|m| !programs.graph.vertices().contains(m)) + .collect::>(), + ); + let null_model = self.build_system(&model, modeldata.clone(), vec![]); + let mut subsystems = affected_mediators + .map(|m| { + let sys = self.build_system(&model, modeldata.clone(), vec![m]); + (programs.clone(), sys) }) - .to_numerical(); - (id, sys) - })); - let sys = NumericalPolynomialSwitchingSystem { subsystems }; - - let ob_index = BTreeMap::new(); - let x0 = DVector::default(); - // ODE Problem expects polynomial - let problem = ODEProblem::new(sys, x0).end_time(data.duration); - ODEAnalysis { - problem, - variable_index: ob_index, + .collect::, PolynomialSystem, u8>)>>(); + subsystems.push((ComputeGraph::::new(), init_model)); + let ob_index: BTreeMap = + null_model.components.keys().cloned().enumerate().map(|(i, x)| (x, i)).collect(); + SwitchingSystem { + null_model, + subsystems, + ob_index, + } } + + // pub fn build_numerical_system( + // &self, + // model: &ModalDblModel, + // data: AnotherMassActionProblemData, + // ) -> ODEAnalysis> { + // // XXX we forget some data from `data` here + // into_numerical_switching_system(self.build_switching_system(model), data.mass) + // } } +// fn into_numerical_switching_system( +// system: SwitchingSystem, +// data: MassActionProblemData, +// ) -> ODEAnalysis> { +// // +// let subsystems = system.subsystems.into_iter().map(|(graph, sys)| { +// let ob_index: BTreeMap<_, _> = +// sys.components.keys().cloned().enumerate().map(|(i, x)| (x, i)).collect(); +// 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 sys = sys +// .extend_scalars(|poly| { +// poly.eval(|flow| data.rates.get(flow).copied().unwrap_or_default()) +// }) +// .to_numerical(); +// (graph, sys) +// }); +// let sys = NumericalPolynomialSwitchingSystem { subsystems: vec![] }; + +// // TODO +// let variable_index = BTreeMap::new(); +// let x0 = DVector::default(); +// // ODE Problem expects polynomial +// let problem = ODEProblem::new(sys, x0).end_time(data.duration); +// ODEAnalysis { +// problem, +// variable_index, +// } +// } + #[cfg(test)] mod tests { use expect_test::expect; @@ -491,4 +756,13 @@ mod tests { // dbg!(&sys); assert!(true); } + + #[test] + fn graph() { + let th = Rc::new(th_modal_state_aux()); + let model = water(th); + let sys = PetriNetMassActionFunctionAnalysis::default().build_switching_system(&model); + // dbg!(&sys); + assert!(true); + } } diff --git a/packages/catlog/src/stdlib/models.rs b/packages/catlog/src/stdlib/models.rs index 808e4cc0c..3e1745332 100644 --- a/packages/catlog/src/stdlib/models.rs +++ b/packages/catlog/src/stdlib/models.rs @@ -161,7 +161,7 @@ pub fn water(th: Rc) -> UstrModalDblModel { borrow.clone(), ); let container = ustr("container"); - model.add_ob(container, state_type.clone()); + model.add_ob(container, aux_type.clone()); let (comparator, lake_sediment) = (ustr("comparator"), ustr("lake_sediment")); model.add_mor( comparator, From 67ee0e85b8903547d8e0d3f99582e87eb6e093aa Mon Sep 17 00:00:00 2001 From: Matt Cuffaro Date: Mon, 4 Aug 2025 20:33:49 -0700 Subject: [PATCH 5/6] WIP: applying changes to frontend --- packages/catlog-wasm/src/analyses.rs | 2 +- packages/catlog-wasm/src/theories.rs | 10 +- .../catlog/src/simulate/ode/polynomial.rs | 11 +- .../src/stdlib/analyses/ode/mass_action.rs | 88 ++++++----- .../catlog/src/stdlib/analyses/ode/mod.rs | 12 -- .../frontend/src/stdlib/analyses/index.ts | 1 + .../stdlib/analyses/switching_mass_action.tsx | 147 ++++++++++++++++++ packages/frontend/src/stdlib/theories.ts | 11 ++ .../src/stdlib/theories/stock-flow-aux.ts | 78 ++++++++++ 9 files changed, 287 insertions(+), 73 deletions(-) create mode 100644 packages/frontend/src/stdlib/analyses/switching_mass_action.tsx create mode 100644 packages/frontend/src/stdlib/theories/stock-flow-aux.ts diff --git a/packages/catlog-wasm/src/analyses.rs b/packages/catlog-wasm/src/analyses.rs index d7d398a3e..e942e846a 100644 --- a/packages/catlog-wasm/src/analyses.rs +++ b/packages/catlog-wasm/src/analyses.rs @@ -23,4 +23,4 @@ pub struct MassActionModelData(pub analyses::ode::MassActionProblemData); #[derive(Serialize, Deserialize, Tsify)] #[tsify(into_wasm_abi, from_wasm_abi)] -pub struct AnotherMassActionModelData(pub analyses::ode::AnotherMassActionProblemData); +pub struct SwitchingMassActionModelData(pub analyses::ode::SwitchingMassActionProblemData); diff --git a/packages/catlog-wasm/src/theories.rs b/packages/catlog-wasm/src/theories.rs index eaf45d652..0ad893bb8 100644 --- a/packages/catlog-wasm/src/theories.rs +++ b/packages/catlog-wasm/src/theories.rs @@ -354,19 +354,15 @@ impl ThModalStateAuxCategory { } /// Simulates a mass-action ODE system with additional configurations. - #[wasm_bindgen(js_name = "stateAuxMassAction")] + #[wasm_bindgen(js_name = "massAction")] pub fn state_aux_mass_action( &self, model: &DblModel, - data: AnotherMassActionModelData, + data: SwitchingMassActionModelData, ) -> Result { - let model: &model::ModalDblModel<_, _> = - (&model.0).try_into().map_err(|_| "Model should be of a modal theory")?; - // Ok(()) Ok(ODEResult( analyses::ode::PetriNetMassActionFunctionAnalysis::default() - .build_numerical_switching_system(model, data.0) - // TODO Switching System + .build_numerical_system(model.modal()?, data.0) .solve_with_defaults() .map_err(|err| format!("{err:?}")) .into(), diff --git a/packages/catlog/src/simulate/ode/polynomial.rs b/packages/catlog/src/simulate/ode/polynomial.rs index e771b0258..85779dd8d 100644 --- a/packages/catlog/src/simulate/ode/polynomial.rs +++ b/packages/catlog/src/simulate/ode/polynomial.rs @@ -199,21 +199,12 @@ where Id: Hash + Clone + Eq + Debug + Ord, Exp: Ord + Clone, { - // TODO accumulate over first arguments + // TODO do we need a `null_model`? fn get_current(self, x: DVector) -> NumericalPolynomialSystem { self.subsystems .into_iter() .filter(|(graph, _)| graph.compute(&x, self.functions.clone(), self.ob_index.clone())) .fold(self.null_model.unwrap(), |acc, (_, sys)| acc + sys) - // if x[0] >= 4.5f32 { - // let (g, sys0) = self.subsystems[0].clone(); - // let (_, sys1) = self.subsystems[1].clone(); - // let _ = &g.compute(&x, self.functions); - // sys0 + sys1 - // } else { - // let (graph, sys) = self.subsystems[1].clone(); - // sys - // } } } diff --git a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs index 0b5520e86..1f69b34bd 100644 --- a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs +++ b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs @@ -300,7 +300,7 @@ impl EligibleFunctions { feature = "serde-wasm", tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object) )] -pub struct AnotherMassActionProblemData +pub struct SwitchingMassActionProblemData where Id: Eq + Hash + Clone + Debug, { @@ -670,50 +670,52 @@ impl PetriNetMassActionFunctionAnalysis { } } - // pub fn build_numerical_system( - // &self, - // model: &ModalDblModel, - // data: AnotherMassActionProblemData, - // ) -> ODEAnalysis> { - // // XXX we forget some data from `data` here - // into_numerical_switching_system(self.build_switching_system(model), data.mass) - // } + pub fn build_numerical_system( + &self, + model: &ModalDblModel, + data: SwitchingMassActionProblemData, + ) -> ODEAnalysis> { + // XXX we forget some data from `data` here + into_numerical_switching_system(self.build_switching_system(model), data) + } } -// fn into_numerical_switching_system( -// system: SwitchingSystem, -// data: MassActionProblemData, -// ) -> ODEAnalysis> { -// // -// let subsystems = system.subsystems.into_iter().map(|(graph, sys)| { -// let ob_index: BTreeMap<_, _> = -// sys.components.keys().cloned().enumerate().map(|(i, x)| (x, i)).collect(); -// 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 sys = sys -// .extend_scalars(|poly| { -// poly.eval(|flow| data.rates.get(flow).copied().unwrap_or_default()) -// }) -// .to_numerical(); -// (graph, sys) -// }); -// let sys = NumericalPolynomialSwitchingSystem { subsystems: vec![] }; - -// // TODO -// let variable_index = BTreeMap::new(); -// let x0 = DVector::default(); -// // ODE Problem expects polynomial -// let problem = ODEProblem::new(sys, x0).end_time(data.duration); -// ODEAnalysis { -// problem, -// variable_index, -// } -// } +fn into_numerical_switching_system( + system: SwitchingSystem, + data: SwitchingMassActionProblemData, +) -> ODEAnalysis> { + let ob_index: BTreeMap<_, _> = system.ob_index.clone(); + let n = ob_index.len(); + let null_model = + Some(system.null_model.clone().extend_scalars(|p| p.eval(|_| 1.0)).to_numerical()); + let initial_values = ob_index + .keys() + .map(|ob| data.mass.initial_values.get(ob).copied().unwrap_or_default()); + let x0 = DVector::from_iterator(n, initial_values); + let subsystems = system + .subsystems + .into_iter() + .map(|(graph, sys)| { + let sys = sys + .extend_scalars(|poly| { + poly.eval(|flow| data.mass.rates.get(flow).copied().unwrap_or_default()) + }) + .to_numerical(); + (graph, sys) + }) + .collect::>(); + let sys = NumericalPolynomialSwitchingSystem { + ob_index: ob_index.clone(), + subsystems, + null_model, + functions: data.functions, + }; + let problem = ODEProblem::new(sys, x0).end_time(data.mass.duration); + ODEAnalysis { + problem, + variable_index: ob_index, + } +} #[cfg(test)] mod tests { diff --git a/packages/catlog/src/stdlib/analyses/ode/mod.rs b/packages/catlog/src/stdlib/analyses/ode/mod.rs index 708ffb3a4..e1ef038bd 100644 --- a/packages/catlog/src/stdlib/analyses/ode/mod.rs +++ b/packages/catlog/src/stdlib/analyses/ode/mod.rs @@ -69,18 +69,6 @@ impl ODEAnalysis { .collect(), }) } - - /// Solves the ODE construed as a switching system - pub fn solve_switching_with_defaults(self) -> Result, IntegrationError> - where - Id: Eq + Hash, - Sys: ODESystem, - { - // if self.variable_index.is_empty() { - // return Ok(Default::default()); - // } - Ok(Default::default()) - } } pub mod linear_ode; diff --git a/packages/frontend/src/stdlib/analyses/index.ts b/packages/frontend/src/stdlib/analyses/index.ts index f1f066d95..03e7729fa 100644 --- a/packages/frontend/src/stdlib/analyses/index.ts +++ b/packages/frontend/src/stdlib/analyses/index.ts @@ -7,3 +7,4 @@ export * from "./model_graph"; export * from "./petri_net_visualization"; export * from "./stock_flow_diagram"; export * from "./submodel_graphs"; +export * from "./switching_mass_action"; diff --git a/packages/frontend/src/stdlib/analyses/switching_mass_action.tsx b/packages/frontend/src/stdlib/analyses/switching_mass_action.tsx new file mode 100644 index 000000000..3ad10968d --- /dev/null +++ b/packages/frontend/src/stdlib/analyses/switching_mass_action.tsx @@ -0,0 +1,147 @@ +import { createMemo } from "solid-js"; + +import type { + DblModel, + ODEResult, + SwitchingMassActionModelData, + SwitchingMassActionProblemData, +} from "catlog-wasm"; +import type { ModelAnalysisProps } from "../../analysis"; +import { + type ColumnSchema, + FixedTableEditor, + Foldable, + createNumericalColumn, +} from "../../components"; +import type { MorphismDecl, ObjectDecl } from "../../model"; +import type { ModelAnalysisMeta } from "../../theory"; +import { ODEResultPlot } from "../../visualization"; +import { createModelODEPlot } from "./simulation"; + +import "./simulation.css"; + +/** Configuration for a mass-action ODE analysis of a model. */ +export type SwitchingMassActionContent = SwitchingMassActionProblemData; + +type Simulator = (model: DblModel, data: SwitchingMassActionModelData) => ODEResult; + +/** Configure a mass-action ODE analysis for use with models of a theory. */ +export function configureSwitchingMassAction(options: { + id?: string; + name?: string; + description?: string; + help?: string; + simulate: Simulator; + isState?: (ob: ObjectDecl) => boolean; + isTransition?: (mor: MorphismDecl) => boolean; +}): ModelAnalysisMeta { + const { + id = "switching-mass-action", + name = "Switching Mass-action dynamics", + description = "Simulate a switching system using the law of mass action", + help = "mass-action", + ...otherOptions + } = options; + return { + id, + name, + description, + help, + component: (props) => , + initialContent: () => ({ + mass: { rates: {}, initialValues: {}, duration: 10 }, + functions: {}, + }), + }; +} + +/** Analyze a model using mass-action dynamics. */ +export function SwitchingMassAction( + props: ModelAnalysisProps & { + simulate: Simulator; + isState?: (ob: ObjectDecl) => boolean; + isTransition?: (mor: MorphismDecl) => boolean; + title?: string; + }, +) { + const obDecls = createMemo(() => { + return props.liveModel + .formalJudgments() + .filter((jgmt) => jgmt.tag === "object") + .filter((ob) => props.isState?.(ob) ?? true); + }, []); + + const morDecls = createMemo(() => { + return props.liveModel + .formalJudgments() + .filter((jgmt) => jgmt.tag === "morphism"); + // .filter((mor) => props.isTransition?.(mor) ?? true); + }, []); + console.log(morDecls()); + + const obSchema: ColumnSchema[] = [ + { + contentType: "string", + header: true, + content: (ob) => ob.name, + }, + createNumericalColumn({ + name: "Initial value", + data: (ob) => props.content.mass.initialValues[ob.id], + validate: (_, data) => data >= 0, + setData: (ob, data) => + props.changeContent((content) => { + content.mass.initialValues[ob.id] = data; + }), + }), + ]; + + const morSchema: ColumnSchema[] = [ + { + contentType: "string", + header: true, + content: (mor) => mor.name, + }, + createNumericalColumn({ + name: "Rate", + data: (mor) => props.content.mass.rates[mor.id], + default: 1, + validate: (_, data) => data >= 0, + setData: (mor, data) => + props.changeContent((content) => { + content.mass.rates[mor.id] = data; + }), + }), + ]; + + // TODO Duration is a simulation-level quantity and should be moved outside of "mass" + const toplevelSchema: ColumnSchema[] = [ + createNumericalColumn({ + name: "Duration", + data: (_) => props.content.mass.duration, + validate: (_, data) => data >= 0, + setData: (_, data) => + props.changeContent((content) => { + content.mass.duration = data; + }), + }), + ]; + + const plotResult = createModelODEPlot( + () => props.liveModel, + (model: DblModel) => props.simulate(model, props.content), + ); + + return ( +
+ +
+ + + +
+ + +
+ ); +} diff --git a/packages/frontend/src/stdlib/theories.ts b/packages/frontend/src/stdlib/theories.ts index 216246196..a495d7de2 100644 --- a/packages/frontend/src/stdlib/theories.ts +++ b/packages/frontend/src/stdlib/theories.ts @@ -84,6 +84,17 @@ stdTheories.add( async () => (await import("./theories/primitive-stock-flow")).default, ); +stdTheories.add( + { + id: "stock-flow-with-auxes", + name: "Stock and flow with logic", + description: "Model accumulation (stocks) and change (flows)", + group: "System Dynamics", + }, + async () => (await import("./theories/stock-flow-aux")).default, +); + + stdTheories.add( { id: "reg-net", diff --git a/packages/frontend/src/stdlib/theories/stock-flow-aux.ts b/packages/frontend/src/stdlib/theories/stock-flow-aux.ts new file mode 100644 index 000000000..6e4b434b7 --- /dev/null +++ b/packages/frontend/src/stdlib/theories/stock-flow-aux.ts @@ -0,0 +1,78 @@ +import { ThModalStateAuxCategory } from "catlog-wasm"; + +import { Theory } from "../../theory"; +import * as analyses from "../analyses"; +import type { TheoryMeta } from "../types"; + +import styles from "../styles.module.css"; +import svgStyles from "../svg_styles.module.css"; + +export default function createStockFlowAuxTheory(theoryMeta: TheoryMeta): Theory { + const thCategoryStockFlowAux = new ThModalStateAuxCategory(); + + return new Theory({ + ...theoryMeta, + theory: thCategoryStockFlowAux.theory(), + onlyFreeModels: true, + modelTypes: [ + { + tag: "ObType", + obType: { tag: "Basic", content: "State" }, + name: "Stock", + description: "Thing with an amount", + shortcut: ["S"], + cssClasses: [styles.box], + svgClasses: [svgStyles.box], + }, + { + tag: "MorType", + morType: { tag: "Basic", content: "out-neg" }, + name: "Outflow", + description: "Flow from one stock to another", + shortcut: ["+"], + arrowStyle: "double", + }, + { + tag: "MorType", + morType: { tag: "Basic", content: "out-pos" }, + name: "Inflow", + description: "Flow from one stock to another", + shortcut: ["-"], + arrowStyle: "double", + }, + { + tag: "ObType", + obType: { tag: "Basic", content: "Auxiliary" }, + name: "Auxiliary", + description: "Variable used in computation", + shortcut: ["A"], + cssClasses: [styles.box], + svgClasses: [svgStyles.box], + }, + { + tag: "MorType", + morType: { tag: "Basic", content: "Borrow" }, + name: "Borrow", + description: "Influence of a stock on a flow", + preferUnnamed: true, + shortcut: ["&"], + }, + ], + modelAnalyses: [ + analyses.configureStockFlowDiagram({ + id: "diagram", + name: "Visualization", + description: "Visualize the stock and flow diagram", + help: "visualization", + }), + analyses.configureSwitchingMassAction({ + simulate(model, data) { + return thCategoryStockFlowAux.massAction(model, data); + }, + isTransition(mor) { + return mor.morType.tag === "Hom"; + }, + }), + ], + }); +} From 6dba66727eb4b48ccbf661339ecf6558b28782fd Mon Sep 17 00:00:00 2001 From: Matt Cuffaro Date: Tue, 5 Aug 2025 18:24:54 -0700 Subject: [PATCH 6/6] ENH: added to frontend --- Cargo.lock | 1 + packages/catlog/Cargo.toml | 1 + .../catlog/src/simulate/ode/polynomial.rs | 66 ++++++++++++++++--- .../src/stdlib/analyses/ode/mass_action.rs | 41 ++++++------ packages/catlog/src/stdlib/models.rs | 58 +++++++++++++++- .../frontend/src/notebook/notebook_editor.tsx | 1 + .../stdlib/analyses/switching_mass_action.tsx | 54 +++++++++++++-- .../src/stdlib/theories/stock-flow-aux.ts | 18 +++-- 8 files changed, 196 insertions(+), 44 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index c57089af0..736ab9f3c 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -380,6 +380,7 @@ dependencies = [ "ustr", "uuid", "wasm-bindgen", + "web-sys", ] [[package]] diff --git a/packages/catlog/Cargo.toml b/packages/catlog/Cargo.toml index 938839524..716afac63 100644 --- a/packages/catlog/Cargo.toml +++ b/packages/catlog/Cargo.toml @@ -28,6 +28,7 @@ tsify = { version = "0.5", features = ["js"], optional = true } ustr = "1" uuid = { version = "=1.11", optional = true } wasm-bindgen = { version = "0.2.100", optional = true } +web-sys = { version = "0.3.77", features = ["console"] } [dev-dependencies] expect-test = "1.5" diff --git a/packages/catlog/src/simulate/ode/polynomial.rs b/packages/catlog/src/simulate/ode/polynomial.rs index 85779dd8d..931197340 100644 --- a/packages/catlog/src/simulate/ode/polynomial.rs +++ b/packages/catlog/src/simulate/ode/polynomial.rs @@ -191,7 +191,7 @@ where pub subsystems: Vec<(ComputeGraph, NumericalPolynomialSystem)>, /// Analysis data - pub functions: HashMap, + pub functions: HashMap, } impl NumericalPolynomialSwitchingSystem @@ -231,12 +231,11 @@ where fn compute( &self, x: &DVector, - functions: HashMap, + functions: HashMap, ob_index: BTreeMap, ) -> bool { // binding name let mut bindings: Vec<(_, Output)> = Default::default(); - // dbg!(&self.borrows, &self.obs); for var in self.toposort.iter() { if self.obs.contains(&var) { let index = ob_index[var]; @@ -252,7 +251,16 @@ where .into_iter() .map(|arg| self.fetch(x, arg, ob_index.clone())) .collect::>(); - let out = functions[&var].clone().eval(args); + let out = if let Some(function) = functions.get(&var) { + let res = &args[0] >= &args[1]; + match function.as_str() { + "Identity" => Output::Bool(true), // do nothing + "Geq" => Output::Bool(args[0] >= args[1]), + _ => Output::Float(args[0]), + } + } else { + Output::Bool(true) + }; bindings.push((var, out)); } } @@ -273,7 +281,7 @@ where } else if self.borrows.keys().contains(&arg) { x[ob_index[&self.borrows[&arg]]] } else { - self.auxes[&arg] + 7f32 // TODO remove } } } @@ -303,7 +311,6 @@ where fn vector_field(&self, dx: &mut DVector, x: &DVector, _t: f32) { let subsystem = self.clone().get_current(x.clone()); for i in 0..dx.len() { - // TODO do we still need to .get() ? dx[i] = match subsystem.components.get(i) { Some(p) => p.eval(|var| x[*var]), None => 0f32, @@ -384,14 +391,54 @@ mod tests { let system = crate::stdlib::analyses::ode::PetriNetMassActionFunctionAnalysis::default() .build_switching_system(&model); - // TODO this should be obtained from the petri net analysis + let functions = HashMap::from([(ustr::ustr("comparator"), String::from("Geq"))]); + let sys = system.to_numerical(functions); + + // container, lake, sediment, watershed + let initial = DVector::from_column_slice(&[6.0, 4.0, 1.0, 4.0]); + // let initial = DVector::from_column_slice(&[1.0, 2.0, 4.0]); + let problem = ODEProblem::new(sys, initial).end_time(10.0); + let result = problem.solve_rk4(0.1).unwrap(); + let expected = expect![[r#" + ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠒⢣⠀⣀⠤⠤⠒⠒⠉⠒⠒⠒⠉⠉⠒⠒⠒⠒⠉⠉⠉⠉⠉⠒⠒⠒⠒⠒⠒⠂ 4.5 + ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠒⠉⠀⠀⠀⠀⠀⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠤⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⣱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡠⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⢇⠀⠀⠀⠀⠀⠀⠀⡠⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠘⡄⠀⠀⠀⠀⠀⡔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠱⡀⠀⠀⠀⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⢣⠀⢀⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⠀⢇⠎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⣀⣀⣀⡀ + ⡁⠀⠀⠀⡎⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⣀⣀⣀⣀⣀⣀⠤⠤⠤⠤⠤⠒⠒⠒⠒⠒⠒⠒⠒⠊⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⡜⠀⠈⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠒⠒⢲⠓⠒⠒⠚⢖⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠒⠊⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⢀⠇⠀⠀⠀⠀⠀⠣⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⡜⠀⠀⠀⠀⠀⠀⠀⠑⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⢲⠁⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠒⠤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠒⠤⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠒⠒⠤⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠒⠒⠒⠤⠤⠤⠤⣀⣀⣀⣀⣀⣀⣀⠀⠀⠀⠀⠀⠀⠀⠀ + ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠉⠉⠉⠉⠁ 0.0 + 0.0 5.0 + "#]]; + expected.assert_eq(&textplot_ode_result(&problem, &result)); + } + + #[test] + fn water_no_logic_simulation() { + let th = std::rc::Rc::new(crate::stdlib::th_modal_state_aux()); + let model = crate::stdlib::water_no_logic(th); + let system = crate::stdlib::analyses::ode::PetriNetMassActionFunctionAnalysis::default() + .build_switching_system(&model); - let functions = HashMap::from([(ustr::ustr("comparator"), EligibleFunctions::Geq())]); + let functions = HashMap::from([(ustr::ustr("comparator"), String::from("Geq"))]); let sys = system.to_numerical(functions); // sediment, watershed, lake let initial = DVector::from_column_slice(&[1.0, 2.0, 4.0]); - let problem = ODEProblem::new(sys, initial).end_time(5.0); + let problem = ODEProblem::new(sys, initial).end_time(10.0); let result = problem.solve_rk4(0.1).unwrap(); let expected = expect![[r#" ⡁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣀⠤⠒⠒⢣⠀⣀⠤⠤⠒⠒⠉⠒⠒⠒⠉⠉⠒⠒⠒⠒⠉⠉⠉⠉⠉⠒⠒⠒⠒⠒⠒⠂ 4.5 @@ -417,7 +464,6 @@ mod tests { ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠉⠉⠉⠉⠁ 0.0 0.0 5.0 "#]]; - // println!("{}", &textplot_ode_result(&problem, &result)); expected.assert_eq(&textplot_ode_result(&problem, &result)); } } diff --git a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs index 1f69b34bd..a95869773 100644 --- a/packages/catlog/src/stdlib/analyses/ode/mass_action.rs +++ b/packages/catlog/src/stdlib/analyses/ode/mass_action.rs @@ -308,7 +308,7 @@ where pub mass: MassActionProblemData, /// functions associated to T(aux) --> aux - pub functions: HashMap, + pub functions: HashMap, } #[derive(Debug, Clone)] @@ -420,7 +420,7 @@ where { pub fn to_numerical( &self, - functions: HashMap, + functions: HashMap, ) -> NumericalPolynomialSwitchingSystem { let subsystems: Vec<(ComputeGraph, _)> = self .subsystems @@ -516,7 +516,7 @@ impl Default for PetriNetMassActionFunctionAnalysis { state_ob_type: ModalObType::new(ustr("State")), aux_ob_type: ModalObType::new(ustr("Auxiliary")), fun_mor_type: ModalMorType::One(ModeApp::new(ustr("function"))), - borrow_mor_type: ModalMorType::One(ModeApp::new(ustr("borrowing"))), + borrow_mor_type: ModalMorType::One(ModeApp::new(ustr("borrow"))), outpos_mor_type: ModalMorType::One(ModeApp::new(ustr("out-pos"))), outneg_mor_type: ModalMorType::One(ModeApp::new(ustr("out-neg"))), } @@ -531,9 +531,9 @@ impl PetriNetMassActionFunctionAnalysis { let mut cg = ComputeGraph::::new(); cg.obs = model.ob_generators_with_type(&self.state_ob_type).collect::>(); // TODO these should just be auxes entering the graph - cg.auxes = HashMap::from_iter( - model.ob_generators_with_type(&self.aux_ob_type).map(|aux| (aux, 4.5f32)), - ); + // cg.auxes = HashMap::from_iter( + // model.ob_generators_with_type(&self.aux_ob_type).map(|aux| (aux, 4.5f32)), + // ); cg.borrows = HashMap::from_iter(model.mor_generators_with_type(&self.borrow_mor_type).map(|f| { ( @@ -643,6 +643,7 @@ impl PetriNetMassActionFunctionAnalysis { .mediators .into_iter() .filter(|m| programs.graph.vertices().contains(m)); + println!("INIT MODEL"); let init_model = self.build_system( &model, modeldata.clone(), @@ -653,16 +654,19 @@ impl PetriNetMassActionFunctionAnalysis { .filter(|m| !programs.graph.vertices().contains(m)) .collect::>(), ); + // println!("DISCRETE MODEL"); let null_model = self.build_system(&model, modeldata.clone(), vec![]); let mut subsystems = affected_mediators .map(|m| { + println!("MEDIATOR: {:?}", &m); let sys = self.build_system(&model, modeldata.clone(), vec![m]); (programs.clone(), sys) }) .collect::, PolynomialSystem, u8>)>>(); - subsystems.push((ComputeGraph::::new(), init_model)); + // subsystems.push((ComputeGraph::::new(), init_model)); let ob_index: BTreeMap = null_model.components.keys().cloned().enumerate().map(|(i, x)| (x, i)).collect(); + let null_model = init_model.clone(); SwitchingSystem { null_model, subsystems, @@ -675,7 +679,6 @@ impl PetriNetMassActionFunctionAnalysis { model: &ModalDblModel, data: SwitchingMassActionProblemData, ) -> ODEAnalysis> { - // XXX we forget some data from `data` here into_numerical_switching_system(self.build_switching_system(model), data) } } @@ -686,8 +689,15 @@ fn into_numerical_switching_system( ) -> ODEAnalysis> { let ob_index: BTreeMap<_, _> = system.ob_index.clone(); let n = ob_index.len(); - let null_model = - Some(system.null_model.clone().extend_scalars(|p| p.eval(|_| 1.0)).to_numerical()); + let null_model = Some( + system + .null_model + .clone() + .extend_scalars(|p| { + p.eval(|flow| data.mass.rates.get(flow).copied().unwrap_or_default()) + }) + .to_numerical(), + ); let initial_values = ob_index .keys() .map(|ob| data.mass.initial_values.get(ob).copied().unwrap_or_default()); @@ -755,16 +765,7 @@ mod tests { let th = Rc::new(th_modal_state_aux()); let model = water(th); let sys = PetriNetMassActionFunctionAnalysis::default().build_switching_system(&model); - // dbg!(&sys); - assert!(true); - } - - #[test] - fn graph() { - let th = Rc::new(th_modal_state_aux()); - let model = water(th); - let sys = PetriNetMassActionFunctionAnalysis::default().build_switching_system(&model); - // dbg!(&sys); + // TODO assert!(true); } } diff --git a/packages/catlog/src/stdlib/models.rs b/packages/catlog/src/stdlib/models.rs index 3e1745332..3daa4abed 100644 --- a/packages/catlog/src/stdlib/models.rs +++ b/packages/catlog/src/stdlib/models.rs @@ -155,17 +155,25 @@ pub fn water(th: Rc) -> UstrModalDblModel { model.add_ob(blake, aux_type.clone()); let borrow = ModalMorType::One(ModeApp::new(ustr("borrowing"))); model.add_mor( - ustr("borrow"), + ustr("borrow-lake"), ModalOb::Generator(lake.clone()), ModalOb::Generator(blake.clone()), borrow.clone(), ); let container = ustr("container"); - model.add_ob(container, aux_type.clone()); + model.add_ob(container, state_type.clone()); + let bcontainer = ustr("&container"); + model.add_ob(bcontainer, aux_type.clone()); + model.add_mor( + ustr("borrow-container"), + ModalOb::Generator(container.clone()), + ModalOb::Generator(bcontainer.clone()), + borrow.clone(), + ); let (comparator, lake_sediment) = (ustr("comparator"), ustr("lake_sediment")); model.add_mor( comparator, - ModalOb::List(List::Plain, vec![blake.into(), container.into()]).into(), + ModalOb::List(List::Plain, vec![blake.into(), bcontainer.into()]).into(), ModalOb::Generator(lake_sediment.clone()), ModalMorType::One(ModeApp::new(ustr("function"))), ); @@ -200,6 +208,50 @@ pub fn water(th: Rc) -> UstrModalDblModel { model } +/** + */ +pub fn water_no_logic(th: Rc) -> UstrModalDblModel { + let (state_type, aux_type) = + (ModalObType::new(ustr("State")), ModalObType::new(ustr("Auxiliary"))); + let mut model = UstrModalDblModel::new(th); + let (watershed, lake, sediment) = (ustr("watershed"), ustr("lake"), ustr("sediment")); + model.add_ob(watershed, state_type.clone()); + model.add_ob(lake, state_type.clone()); + model.add_ob(sediment, state_type.clone()); + let container = ustr("container"); + model.add_ob(container, aux_type.clone()); + let lake_sediment = ustr("lake_sediment"); + // watershed outflow + let watershed_lake = ustr("watershed_lake"); + model.add_ob(watershed_lake.clone(), aux_type.clone()); + model.add_mor( + ustr("watershed-lake-outflow"), + ModalOb::Generator(watershed_lake.clone()), + ModalOb::Generator(watershed.clone()), + ModalMorType::One(ModeApp::new(ustr("out-neg"))), + ); + model.add_mor( + ustr("watershed-lake-inflow"), + ModalOb::Generator(watershed_lake.clone()), + ModalOb::Generator(lake.clone()), + ModalMorType::One(ModeApp::new(ustr("out-pos"))), + ); + // + model.add_mor( + ustr("lake-sediment-outflow"), + ModalOb::Generator(lake_sediment.clone()), + ModalOb::Generator(lake.clone()), + ModalMorType::One(ModeApp::new(ustr("out-neg"))), + ); + model.add_mor( + ustr("lake-sediment-infow"), + ModalOb::Generator(lake_sediment.clone()), + ModalOb::Generator(sediment.clone()), + ModalMorType::One(ModeApp::new(ustr("out-pos"))), + ); + model +} + #[cfg(test)] mod tests { use super::super::theories::*; diff --git a/packages/frontend/src/notebook/notebook_editor.tsx b/packages/frontend/src/notebook/notebook_editor.tsx index 25bc7577f..c278e5007 100644 --- a/packages/frontend/src/notebook/notebook_editor.tsx +++ b/packages/frontend/src/notebook/notebook_editor.tsx @@ -301,6 +301,7 @@ export function NotebookEditor(props: { }; const cell = props.notebook.cellContents[cellId]; + console.log(props.notebook.cellContents); invariant(cell, `Failed to find contents for cell '${cellId}'`); if (cell.tag !== "rich-text") { diff --git a/packages/frontend/src/stdlib/analyses/switching_mass_action.tsx b/packages/frontend/src/stdlib/analyses/switching_mass_action.tsx index 3ad10968d..1c6279311 100644 --- a/packages/frontend/src/stdlib/analyses/switching_mass_action.tsx +++ b/packages/frontend/src/stdlib/analyses/switching_mass_action.tsx @@ -64,20 +64,33 @@ export function SwitchingMassAction( title?: string; }, ) { + console.log(props.liveModel.formalJudgments()); const obDecls = createMemo(() => { return props.liveModel .formalJudgments() .filter((jgmt) => jgmt.tag === "object") + .filter((ob) => ob.obType.content === "State") .filter((ob) => props.isState?.(ob) ?? true); }, []); - const morDecls = createMemo(() => { + const auxDecls = createMemo(() => { return props.liveModel .formalJudgments() - .filter((jgmt) => jgmt.tag === "morphism"); - // .filter((mor) => props.isTransition?.(mor) ?? true); + .filter((jgmt) => jgmt.tag === "object") + .filter((ob) => ob.obType.content === "Auxiliary"); }, []); - console.log(morDecls()); + + // const morDecls = createMemo(() => { + // return props.liveModel + // .formalJudgments() + // .filter((jgmt) => jgmt.tag === "morphism") + // // .filter((mor) => mor.morType.content === "out-pos") + // .filter((mor) => props.isTransition?.(mor) ?? true); + // }, []); + + const functionDecls = createMemo(() => { + return props.liveModel.formalJudgments().filter((jgmt) => jgmt.tag === "morphism").filter((mor) => mor.morType.content === "function").filter((mor) => props.isTransition?.(mor) ?? true); + }, []); const obSchema: ColumnSchema[] = [ { @@ -96,7 +109,7 @@ export function SwitchingMassAction( }), ]; - const morSchema: ColumnSchema[] = [ + const auxSchema: ColumnSchema[] = [ { contentType: "string", header: true, @@ -109,6 +122,7 @@ export function SwitchingMassAction( validate: (_, data) => data >= 0, setData: (mor, data) => props.changeContent((content) => { + console.log(content.mass.rates); content.mass.rates[mor.id] = data; }), }), @@ -127,6 +141,33 @@ export function SwitchingMassAction( }), ]; + + + const functionSchema: ColumnSchema[] = [ + { + contentType: "string", + header: true, + content: (ob) => ob.name, + }, + { + contentType: "enum", + name: "Function", + variants() { + return ["Identity", "Geq"]; + }, + content: () => "Geq", + setContent: (ob, value) => + props.changeContent((content) => { + if (value === null) { + delete content.functions[ob.id]; + } else { + content.functions[ob.id] = value; + } + }), + }, + // TODO createEnumColumn + ]; + const plotResult = createModelODEPlot( () => props.liveModel, (model: DblModel) => props.simulate(model, props.content), @@ -137,7 +178,8 @@ export function SwitchingMassAction(
- + +
diff --git a/packages/frontend/src/stdlib/theories/stock-flow-aux.ts b/packages/frontend/src/stdlib/theories/stock-flow-aux.ts index 6e4b434b7..89846ba1d 100644 --- a/packages/frontend/src/stdlib/theories/stock-flow-aux.ts +++ b/packages/frontend/src/stdlib/theories/stock-flow-aux.ts @@ -46,17 +46,25 @@ export default function createStockFlowAuxTheory(theoryMeta: TheoryMeta): Theory name: "Auxiliary", description: "Variable used in computation", shortcut: ["A"], - cssClasses: [styles.box], - svgClasses: [svgStyles.box], + // cssClasses: [styles.box], + // svgClasses: [svgStyles.box], }, { tag: "MorType", - morType: { tag: "Basic", content: "Borrow" }, + morType: { tag: "Basic", content: "borrow" }, name: "Borrow", description: "Influence of a stock on a flow", preferUnnamed: true, shortcut: ["&"], }, + { + tag: "MorType", + morType: { tag: "Basic", content: "function" }, + name: "Function", + description: "Function consuming auxiliary variables", + preferUnnamed: true, + shortcut: ["f"], + } ], modelAnalyses: [ analyses.configureStockFlowDiagram({ @@ -67,10 +75,10 @@ export default function createStockFlowAuxTheory(theoryMeta: TheoryMeta): Theory }), analyses.configureSwitchingMassAction({ simulate(model, data) { - return thCategoryStockFlowAux.massAction(model, data); + return thCategoryStockFlowAux.massAction(model, data); }, isTransition(mor) { - return mor.morType.tag === "Hom"; + return mor.morType.tag === "Basic"; }, }), ],