catlog/stdlib/analyses/ode/
kuramoto.rs

1//! Kuramoto ODE analsysis of models.
2
3use std::collections::HashMap;
4
5use indexmap::IndexMap;
6use nalgebra::{DMatrix, DVector};
7
8#[cfg(feature = "serde")]
9use serde::{Deserialize, Serialize};
10#[cfg(feature = "serde-wasm")]
11use tsify::Tsify;
12
13use super::{ODEAnalysis, ODEProblem};
14use crate::dbl::model::{DiscreteDblModel, FgDblModel};
15use crate::one::{FgCategory, QualifiedPath};
16use crate::simulate::ode::{KuramotoOrder, KuramotoSystem};
17use crate::zero::QualifiedName;
18
19/// Data defining a Kuramoto ODE problem for a model.
20#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
21#[cfg_attr(feature = "serde", serde(tag = "order"))]
22#[cfg_attr(feature = "serde-wasm", derive(Tsify))]
23#[cfg_attr(
24    feature = "serde-wasm",
25    tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object)
26)]
27pub enum KuramotoProblemData {
28    /// Data defining a first-order Kuramoto ODE problem.
29    #[serde(rename = "first")]
30    FirstOrder(CommonKuramotoProblemData),
31
32    /// Data defining a second-order Kuramoto ODE problem.
33    #[serde(rename = "second")]
34    SecondOrder {
35        /// Data common to any Kuramoto problem.
36        #[cfg_attr(feature = "serde", serde(flatten))]
37        common: CommonKuramotoProblemData,
38
39        /// Map from object IDs to initial values of angular frequencies.
40        #[cfg_attr(feature = "serde", serde(rename = "initialFrequencies"))]
41        initial_frequencies: HashMap<QualifiedName, f32>,
42    },
43}
44
45impl KuramotoProblemData {
46    fn common(&self) -> &CommonKuramotoProblemData {
47        match self {
48            Self::FirstOrder(common) => common,
49            Self::SecondOrder { common, .. } => common,
50        }
51    }
52
53    fn order(&self) -> KuramotoOrder {
54        match self {
55            Self::FirstOrder(_) => KuramotoOrder::First,
56            Self::SecondOrder { .. } => KuramotoOrder::Second,
57        }
58    }
59}
60
61/// Data common to both first- and second-order Kuramoto problems.
62#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
63#[cfg_attr(feature = "serde-wasm", derive(Tsify))]
64#[cfg_attr(
65    feature = "serde-wasm",
66    tsify(into_wasm_abi, from_wasm_abi, hashmap_as_object)
67)]
68pub struct CommonKuramotoProblemData {
69    /// Map from morphism IDs to coupling coefficients (nonnegative reals).
70    #[cfg_attr(feature = "serde", serde(rename = "couplingCoefficients"))]
71    coupling_coeffs: HashMap<QualifiedName, f32>,
72
73    /// Map from object IDs to damping coefficients (nonnegative reals).
74    #[cfg_attr(feature = "serde", serde(rename = "dampingCoefficients"))]
75    damping_coeffs: HashMap<QualifiedName, f32>,
76
77    /// Map from object IDs to forcing parameters (reals).
78    #[cfg_attr(feature = "serde", serde(rename = "forcingParameters"))]
79    forcing_params: HashMap<QualifiedName, f32>,
80
81    /// Map from object IDs to initial values of phases.
82    #[cfg_attr(feature = "serde", serde(rename = "initialPhases"))]
83    initial_phases: HashMap<QualifiedName, f32>,
84
85    /// Duration of simulation.
86    duration: f32,
87}
88
89/// Kuramoto ODE analysis of a model.
90pub struct KuramotoAnalysis {
91    node_ob_type: QualifiedName,
92    link_mor_types: Vec<QualifiedPath>,
93}
94
95impl KuramotoAnalysis {
96    /// Constructs a Kuramoto analysis with nodes the objects of given type.
97    pub fn new(ob_type: QualifiedName) -> Self {
98        Self {
99            node_ob_type: ob_type,
100            link_mor_types: Default::default(),
101        }
102    }
103
104    /// Adds a type of morphism to be treated as links between nodes.
105    pub fn add_link_type(mut self, mor_type: QualifiedPath) -> Self {
106        self.link_mor_types.push(mor_type);
107        self
108    }
109
110    /// Creates a Kuramoto system from a model plus numerical data.
111    pub fn build_system(
112        &self,
113        model: &DiscreteDblModel,
114        data: &KuramotoProblemData,
115    ) -> ODEAnalysis<KuramotoSystem> {
116        let common = data.common();
117
118        let ob_index: IndexMap<_, _> = model
119            .ob_generators_with_type(&self.node_ob_type)
120            .enumerate()
121            .map(|(i, x)| (x, i))
122            .collect();
123        let ob_ids = || ob_index.keys();
124        let n = ob_index.len();
125
126        let mut coupling_coeffs = DMatrix::from_element(n, n, 0.0f32);
127        for mor_type in self.link_mor_types.iter() {
128            for mor in model.mor_generators_with_type(mor_type) {
129                let coef = common.coupling_coeffs.get(&mor).copied().unwrap_or_default();
130                let (dom, cod) = (model.mor_generator_dom(&mor), model.mor_generator_cod(&mor));
131                let (i, j) = (*ob_index.get(&dom).unwrap(), *ob_index.get(&cod).unwrap());
132                // Coupling matrix is assumed symmetric in the literature.
133                coupling_coeffs[(i, j)] += coef;
134                coupling_coeffs[(j, i)] += coef;
135            }
136        }
137
138        let damping_coeffs =
139            ob_ids().map(|id| common.damping_coeffs.get(id).copied().unwrap_or_default());
140        let forcing_params =
141            ob_ids().map(|id| common.forcing_params.get(id).copied().unwrap_or_default());
142
143        let system = KuramotoSystem {
144            order: data.order(),
145            coupling_coeffs,
146            damping_coeffs: DVector::from_iterator(n, damping_coeffs),
147            forcing_params: DVector::from_iterator(n, forcing_params),
148        };
149        let initial_phases =
150            ob_ids().map(|id| common.initial_phases.get(id).copied().unwrap_or_default());
151        let initial_values = match data {
152            KuramotoProblemData::FirstOrder(_) => DVector::from_iterator(n, initial_phases),
153            KuramotoProblemData::SecondOrder {
154                initial_frequencies,
155                ..
156            } => {
157                let initial_frequencies =
158                    ob_ids().map(|id| initial_frequencies.get(id).copied().unwrap_or_default());
159                DVector::from_iterator(2 * n, initial_phases.chain(initial_frequencies))
160            }
161        };
162        let problem = ODEProblem::new(system, initial_values).end_time(common.duration);
163        ODEAnalysis::new(problem, ob_index)
164    }
165}