Modal Double Theory with State and Auxiliary Variables

packages/catlog-wasm/src/analyses.rs

@@ -20,3 +20,7 @@ pub struct LinearODEModelData(pub analyses::ode::LinearODEProblemData<Uuid>);
 #[derive(Serialize, Deserialize, Tsify)]
 #[tsify(into_wasm_abi, from_wasm_abi)]
 pub struct MassActionModelData(pub analyses::ode::MassActionProblemData<Uuid>);
+
+#[derive(Serialize, Deserialize, Tsify)]
+#[tsify(into_wasm_abi, from_wasm_abi)]
+pub struct SwitchingMassActionModelData(pub analyses::ode::SwitchingMassActionProblemData<Uuid>);

packages/catlog-wasm/src/theories.rs

@@ -10,10 +10,10 @@
 use wasm_bindgen::prelude::*;
 
 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.
@@ -337,6 +337,39 @@
     }
 }
 
+/// The theory of state/aux interactions
+#[wasm_bindgen]
+pub struct ThModalStateAuxCategory(Rc<theory::UstrModalDblTheory>);
+
+#[wasm_bindgen]
+impl ThModalStateAuxCategory {
+    #[wasm_bindgen(constructor)]
+    pub fn new() -> Self {
+        Self(Rc::new(theories::th_modal_state_aux()))
+    }
+
+    #[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 = "massAction")]
+    pub fn state_aux_mass_action(
+        &self,
+        model: &DblModel,
+        data: SwitchingMassActionModelData,
+    ) -> Result<ODEResult, String> {
+        Ok(ODEResult(
+            analyses::ode::PetriNetMassActionFunctionAnalysis::default()
+                .build_numerical_system(model.modal()?, data.0)
+                .solve_with_defaults()
+                .map_err(|err| format!("{err:?}"))
+                .into(),
+        ))
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;

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"

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<Id, ThId> {
     /// Generating object.
     #[from]

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,

packages/catlog/src/one/graph.rs

@@ -191,23 +191,23 @@
 Such a graph is defined in copresheaf style by two [finite sets](FinSet) and two
 [columns](Column). Implementing this trait provides a *blanket implementation*
 of [`FinGraph`].
  */
 pub trait ColumnarFinGraph:
     ColumnarGraph<
-        Vertices: FinSet<Elem = Self::V>,
-        Edges: FinSet<Elem = Self::E>,
-        Src: Column<Dom = Self::E, Cod = Self::V>,
-        Tgt: Column<Dom = Self::E, Cod = Self::V>,
-    >
+    Vertices: FinSet<Elem = Self::V>,
+    Edges: FinSet<Elem = Self::E>,
+    Src: Column<Dom = Self::E, Cod = Self::V>,
+    Tgt: Column<Dom = Self::E, Cod = Self::V>,
+>
 {
 }
 
 /// A columnar graph with mutable columns.
 pub trait MutColumnarGraph:
     ColumnarGraph<
-        Src: MutMapping<Dom = Self::E, Cod = Self::V>,
-        Tgt: MutMapping<Dom = Self::E, Cod = Self::V>,
-    >
+    Src: MutMapping<Dom = Self::E, Cod = Self::V>,
+    Tgt: MutMapping<Dom = Self::E, Cod = Self::V>,
+>
 {
     /// Variant of [`src_map`](ColumnarGraph::src_map) that returns a mutable
     /// reference.

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 @@
     result
 }
 
+fn out_neighbors<G>(graph: &G, v: &G::V) -> impl Iterator<Item = G::V>
+where
+    G: FinGraph,
+    G::V: Hash,
+{
+    graph.out_edges(v).map(|e| graph.tgt(&e))
+}
+
+fn in_neighbors<G>(graph: &G, v: &G::V) -> impl Iterator<Item = G::V>
+where
+    G: FinGraph,
+    G::V: Hash,
+{
+    graph.in_edges(v).map(|e| graph.src(&e))
+}
+
+#[derive(Clone, Debug)]
+struct VisitMap {
+    visited: Vec<bool>,
+}
+
+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<G>
+where
+    G: FinGraph,
+    G::V: Hash,
+{
+    //
+    pub stack: Vec<G::V>,
+    //
+    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<G>(graph: &G) -> Result<Vec<G::V>, String>
+where
+    G: FinGraph,
+    G::V: Hash + std::fmt::Debug,
+{
+    // XXX dont clone
+    let n = graph.vertices().collect::<Vec<_>>().len();
+    let mut discovered = VisitMap {
+        visited: vec![false; n.clone()],
+    };
+    let mut finished = VisitMap {
+        visited: vec![false; n.clone()],
+    };
+    let mut finish_stack: Vec<G::V> = 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 @@
         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<u8, &str> = 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::<Vec<_>>();
+        // dbg!(&out_neighbors);
+        assert!(true);
+    }
 }

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 @@
     }
 }
 
+// --------------------------------
+
+// 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<f32>) -> bool {
+        match self {
+            Condition::Geq(left, right) => x.row(*left) >= x.row(*right),
+        }
+    }
+}
+
+/**
+ */
+#[derive(Clone)]
+pub struct HybridLinearODESystem {
+    /// Linear ODEs
+    pub subsystems: HashMap<Condition, LinearODESystem>,
+}
+
+impl HybridLinearODESystem {
+    // TODO I dont like checking for the first condition to be met.
+    pub fn get_subsystem(&self, x: &DVector<f32>) -> 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<f32>, x: &DVector<f32>, _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<LinearODESystem> {
     use nalgebra::{dmatrix, dvector};