oldies_xppaut/
lib.rs

1//! # XPPAUT-RS
2//!
3//! Revival of XPPAUT (XPP + AUTO) bifurcation analysis in Rust.
4//!
5//! ## History
6//!
7//! XPPAUT was created by Bard Ermentrout at University of Pittsburgh.
8//! It's described by the author as "sort of a hobby of mine" with a
9//! "dated user interface" that "sometimes crashes" and has "no scripting
10//! interface."
11//!
12//! This crate provides a modern, safe, high-performance implementation
13//! of XPPAUT's bifurcation analysis capabilities.
14//!
15//! ## Capabilities
16//!
17//! 1. **ODE Integration**: Multiple solvers (Euler, RK4, adaptive)
18//! 2. **Bifurcation Analysis**: Fixed points, limit cycles, Hopf bifurcations
19//! 3. **Continuation**: Parameter continuation via AUTO algorithm
20//! 4. **Phase Portraits**: Nullclines, vector fields
21//! 5. **Stability Analysis**: Eigenvalues, Floquet multipliers
22
23use oldies_core::{OdeSystem, Result, StateVector, Time, OldiesError};
24use nalgebra::{DMatrix, DVector};
25use ndarray::Array1;
26use num_complex::Complex64;
27use serde::{Deserialize, Serialize};
28
29/// Bifurcation types
30#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
31pub enum BifurcationType {
32    /// Saddle-node (fold)
33    SaddleNode,
34    /// Transcritical
35    Transcritical,
36    /// Pitchfork
37    Pitchfork,
38    /// Hopf (supercritical or subcritical)
39    Hopf { supercritical: bool },
40    /// Period-doubling
41    PeriodDoubling,
42    /// Limit point of cycles
43    LimitPointCycles,
44    /// Torus (Neimark-Sacker)
45    Torus,
46}
47
48/// Fixed point with stability info
49#[derive(Debug, Clone, Serialize, Deserialize)]
50pub struct FixedPoint {
51    /// State at fixed point
52    pub state: Vec<f64>,
53    /// Parameter value
54    pub parameter: f64,
55    /// Eigenvalues
56    pub eigenvalues: Vec<Complex64>,
57    /// Is stable?
58    pub stable: bool,
59    /// Type (node, focus, saddle, etc.)
60    pub point_type: FixedPointType,
61}
62
63/// Fixed point classification
64#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
65pub enum FixedPointType {
66    StableNode,
67    UnstableNode,
68    StableFocus,
69    UnstableFocus,
70    Saddle,
71    Center,
72    Unknown,
73}
74
75/// Limit cycle with stability info
76#[derive(Debug, Clone, Serialize, Deserialize)]
77pub struct LimitCycle {
78    /// Period
79    pub period: f64,
80    /// Amplitude (max deviation from mean)
81    pub amplitude: f64,
82    /// Parameter value
83    pub parameter: f64,
84    /// Floquet multipliers
85    pub floquet_multipliers: Vec<Complex64>,
86    /// Is stable?
87    pub stable: bool,
88}
89
90/// Bifurcation point
91#[derive(Debug, Clone, Serialize, Deserialize)]
92pub struct BifurcationPoint {
93    /// Type of bifurcation
94    pub bifurcation_type: BifurcationType,
95    /// Parameter value at bifurcation
96    pub parameter: f64,
97    /// State at bifurcation
98    pub state: Vec<f64>,
99    /// Additional info (e.g., Hopf frequency)
100    pub info: Option<String>,
101}
102
103/// Bifurcation diagram
104#[derive(Debug, Clone, Serialize, Deserialize)]
105pub struct BifurcationDiagram {
106    /// Parameter name
107    pub parameter_name: String,
108    /// Parameter range
109    pub parameter_range: (f64, f64),
110    /// State variable index being tracked
111    pub state_index: usize,
112    /// Fixed point branches
113    pub fixed_points: Vec<FixedPoint>,
114    /// Limit cycle branches
115    pub limit_cycles: Vec<LimitCycle>,
116    /// Bifurcation points
117    pub bifurcations: Vec<BifurcationPoint>,
118}
119
120/// XPPAUT-style ODE model
121#[derive(Debug, Clone)]
122pub struct XppModel {
123    /// Model name
124    pub name: String,
125    /// Variable names
126    pub variables: Vec<String>,
127    /// Parameter names and values
128    pub parameters: Vec<(String, f64)>,
129    /// ODE right-hand sides (as Rust closures would require Box<dyn Fn>)
130    dimension: usize,
131}
132
133impl XppModel {
134    /// Create a new model
135    pub fn new(name: &str, variables: Vec<String>) -> Self {
136        let dimension = variables.len();
137        Self {
138            name: name.to_string(),
139            variables,
140            parameters: Vec::new(),
141            dimension,
142        }
143    }
144
145    /// Add a parameter
146    pub fn add_parameter(&mut self, name: &str, value: f64) {
147        self.parameters.push((name.to_string(), value));
148    }
149
150    /// Get parameter value
151    pub fn get_parameter(&self, name: &str) -> Option<f64> {
152        self.parameters.iter()
153            .find(|(n, _)| n == name)
154            .map(|(_, v)| *v)
155    }
156
157    /// Set parameter value
158    pub fn set_parameter(&mut self, name: &str, value: f64) -> Result<()> {
159        for (n, v) in &mut self.parameters {
160            if n == name {
161                *v = value;
162                return Ok(());
163            }
164        }
165        Err(OldiesError::ModelNotFound(format!("Parameter {} not found", name)))
166    }
167}
168
169/// Bifurcation analyzer
170pub struct BifurcationAnalyzer {
171    /// Model
172    model: XppModel,
173    /// Numerical tolerance
174    tolerance: f64,
175    /// Maximum iterations for Newton's method
176    max_iterations: usize,
177}
178
179impl BifurcationAnalyzer {
180    /// Create a new analyzer
181    pub fn new(model: XppModel) -> Self {
182        Self {
183            model,
184            tolerance: 1e-10,
185            max_iterations: 100,
186        }
187    }
188
189    /// Find fixed points at current parameter values
190    pub fn find_fixed_points<F>(&self, rhs: F, initial_guesses: &[Vec<f64>]) -> Vec<FixedPoint>
191    where
192        F: Fn(&[f64], &[(String, f64)]) -> Vec<f64>,
193    {
194        let mut fixed_points = Vec::new();
195
196        for guess in initial_guesses {
197            if let Some(fp) = self.newton_raphson(&rhs, guess) {
198                // Check if we already found this one
199                let is_new = fixed_points.iter().all(|existing: &FixedPoint| {
200                    let dist: f64 = existing.state.iter()
201                        .zip(&fp)
202                        .map(|(a, b)| (a - b).powi(2))
203                        .sum::<f64>()
204                        .sqrt();
205                    dist > self.tolerance * 100.0
206                });
207
208                if is_new {
209                    // Compute Jacobian and eigenvalues
210                    let jacobian = self.numerical_jacobian(&rhs, &fp);
211                    let eigenvalues = self.compute_eigenvalues(&jacobian);
212                    let stable = eigenvalues.iter().all(|e| e.re < 0.0);
213                    let point_type = classify_fixed_point(&eigenvalues);
214
215                    fixed_points.push(FixedPoint {
216                        state: fp,
217                        parameter: 0.0, // Would need parameter tracking
218                        eigenvalues,
219                        stable,
220                        point_type,
221                    });
222                }
223            }
224        }
225
226        fixed_points
227    }
228
229    /// Newton-Raphson for finding zeros
230    fn newton_raphson<F>(&self, rhs: F, initial: &[f64]) -> Option<Vec<f64>>
231    where
232        F: Fn(&[f64], &[(String, f64)]) -> Vec<f64>,
233    {
234        let mut x = initial.to_vec();
235        let n = x.len();
236
237        for _ in 0..self.max_iterations {
238            let f = rhs(&x, &self.model.parameters);
239
240            // Check convergence
241            let norm: f64 = f.iter().map(|v| v.powi(2)).sum::<f64>().sqrt();
242            if norm < self.tolerance {
243                return Some(x);
244            }
245
246            // Compute Jacobian
247            let jacobian = self.numerical_jacobian(&rhs, &x);
248
249            // Solve J * delta = -f
250            let j = DMatrix::from_row_slice(n, n, &jacobian);
251            let f_vec = DVector::from_vec(f.iter().map(|v| -*v).collect());
252
253            if let Some(lu) = j.lu().solve(&f_vec) {
254                for i in 0..n {
255                    x[i] += lu[i];
256                }
257            } else {
258                return None; // Singular Jacobian
259            }
260        }
261
262        None
263    }
264
265    /// Numerical Jacobian via finite differences
266    fn numerical_jacobian<F>(&self, rhs: F, x: &[f64]) -> Vec<f64>
267    where
268        F: Fn(&[f64], &[(String, f64)]) -> Vec<f64>,
269    {
270        let n = x.len();
271        let h = 1e-8;
272        let mut jacobian = vec![0.0; n * n];
273
274        let f0 = rhs(x, &self.model.parameters);
275
276        for j in 0..n {
277            let mut x_plus = x.to_vec();
278            x_plus[j] += h;
279            let f_plus = rhs(&x_plus, &self.model.parameters);
280
281            for i in 0..n {
282                jacobian[i * n + j] = (f_plus[i] - f0[i]) / h;
283            }
284        }
285
286        jacobian
287    }
288
289    /// Compute eigenvalues of a matrix
290    fn compute_eigenvalues(&self, matrix: &[f64]) -> Vec<Complex64> {
291        let n = (matrix.len() as f64).sqrt() as usize;
292        let m = DMatrix::from_row_slice(n, n, matrix);
293
294        // Use nalgebra's eigenvalue computation
295        if let Some(eigen) = m.clone().try_symmetric_eigen(1e-10, 1000) {
296            eigen.eigenvalues.iter()
297                .map(|&v| Complex64::new(v, 0.0))
298                .collect()
299        } else {
300            // Fall back to Schur decomposition for non-symmetric
301            vec![Complex64::new(0.0, 0.0); n]
302        }
303    }
304}
305
306/// Classify fixed point based on eigenvalues
307fn classify_fixed_point(eigenvalues: &[Complex64]) -> FixedPointType {
308    let all_real = eigenvalues.iter().all(|e| e.im.abs() < 1e-10);
309    let all_negative = eigenvalues.iter().all(|e| e.re < 0.0);
310    let all_positive = eigenvalues.iter().all(|e| e.re > 0.0);
311    let any_zero = eigenvalues.iter().any(|e| e.re.abs() < 1e-10);
312
313    if any_zero {
314        FixedPointType::Center
315    } else if all_real {
316        if all_negative {
317            FixedPointType::StableNode
318        } else if all_positive {
319            FixedPointType::UnstableNode
320        } else {
321            FixedPointType::Saddle
322        }
323    } else {
324        if all_negative {
325            FixedPointType::StableFocus
326        } else if all_positive {
327            FixedPointType::UnstableFocus
328        } else {
329            FixedPointType::Saddle
330        }
331    }
332}
333
334/// Common dynamical systems
335pub mod examples {
336    use super::*;
337
338    /// Lorenz system
339    pub fn lorenz(sigma: f64, rho: f64, beta: f64) -> XppModel {
340        let mut model = XppModel::new("Lorenz", vec!["x".into(), "y".into(), "z".into()]);
341        model.add_parameter("sigma", sigma);
342        model.add_parameter("rho", rho);
343        model.add_parameter("beta", beta);
344        model
345    }
346
347    /// Lorenz RHS
348    pub fn lorenz_rhs(state: &[f64], params: &[(String, f64)]) -> Vec<f64> {
349        let x = state[0];
350        let y = state[1];
351        let z = state[2];
352
353        let sigma = params.iter().find(|(n, _)| n == "sigma").map(|(_, v)| *v).unwrap_or(10.0);
354        let rho = params.iter().find(|(n, _)| n == "rho").map(|(_, v)| *v).unwrap_or(28.0);
355        let beta = params.iter().find(|(n, _)| n == "beta").map(|(_, v)| *v).unwrap_or(8.0/3.0);
356
357        vec![
358            sigma * (y - x),
359            x * (rho - z) - y,
360            x * y - beta * z,
361        ]
362    }
363
364    /// FitzHugh-Nagumo model
365    pub fn fitzhugh_nagumo(a: f64, b: f64, epsilon: f64) -> XppModel {
366        let mut model = XppModel::new("FitzHugh-Nagumo", vec!["v".into(), "w".into()]);
367        model.add_parameter("a", a);
368        model.add_parameter("b", b);
369        model.add_parameter("epsilon", epsilon);
370        model
371    }
372}
373
374#[cfg(test)]
375mod tests {
376    use super::*;
377
378    #[test]
379    fn test_lorenz_model() {
380        let model = examples::lorenz(10.0, 28.0, 8.0/3.0);
381        assert_eq!(model.variables.len(), 3);
382        assert_eq!(model.get_parameter("sigma"), Some(10.0));
383    }
384
385    #[test]
386    fn test_lorenz_rhs() {
387        let model = examples::lorenz(10.0, 28.0, 8.0/3.0);
388        let state = vec![1.0, 1.0, 1.0];
389        let rhs = examples::lorenz_rhs(&state, &model.parameters);
390
391        assert_eq!(rhs.len(), 3);
392        assert!((rhs[0] - 0.0).abs() < 1e-10); // sigma*(y-x) = 10*(1-1) = 0
393    }
394
395    #[test]
396    fn test_fixed_point_classification() {
397        // Stable node (all negative real)
398        let eig = vec![Complex64::new(-1.0, 0.0), Complex64::new(-2.0, 0.0)];
399        assert_eq!(classify_fixed_point(&eig), FixedPointType::StableNode);
400
401        // Saddle (mixed)
402        let eig = vec![Complex64::new(-1.0, 0.0), Complex64::new(1.0, 0.0)];
403        assert_eq!(classify_fixed_point(&eig), FixedPointType::Saddle);
404    }
405}
oldies_xppaut/lib.rs

oldies_xppaut/
lib.rs
