koopman_dmd/
types.rs

1use faer::Mat;
2
3use crate::lifting::{LiftingConfig, LiftingInfo};
4
5/// Error types for DMD operations.
6#[derive(Debug, thiserror::Error)]
7pub enum DmdError {
8    #[error("invalid input: {0}")]
9    InvalidInput(String),
10
11    #[error("SVD computation failed: {0}")]
12    SvdFailed(String),
13
14    #[error("eigendecomposition failed: {0}")]
15    EigenFailed(String),
16
17    #[error("linear solve failed: {0}")]
18    SolveFailed(String),
19
20    #[error("numerical error: {0}")]
21    NumericalError(String),
22}
23
24/// Configuration for DMD computation.
25#[derive(Debug, Clone)]
26pub struct DmdConfig {
27    /// Truncation rank. None for automatic selection (99% variance).
28    pub rank: Option<usize>,
29    /// Whether to center the data (subtract row means).
30    pub center: bool,
31    /// Time step between snapshots.
32    pub dt: f64,
33    /// Optional lifting transformation for Extended DMD.
34    pub lifting: Option<LiftingConfig>,
35}
36
37impl Default for DmdConfig {
38    fn default() -> Self {
39        Self {
40            rank: None,
41            center: false,
42            dt: 1.0,
43            lifting: None,
44        }
45    }
46}
47
48/// Components of the truncated SVD.
49#[derive(Debug, Clone)]
50pub struct SvdComponents {
51    /// Left singular vectors (m × r).
52    pub u: Mat<f64>,
53    /// Singular values (r), stored as a column vector.
54    pub s: Vec<f64>,
55    /// Right singular vectors (n × r), columns are right singular vectors.
56    pub v: Mat<f64>,
57}
58
59/// Complex number type (re, im).
60#[derive(Debug, Clone, Copy)]
61pub struct C64 {
62    pub re: f64,
63    pub im: f64,
64}
65
66impl C64 {
67    /// Create a new complex number.
68    pub fn new(re: f64, im: f64) -> Self {
69        Self { re, im }
70    }
71
72    /// Magnitude |z| = sqrt(re² + im²).
73    pub fn norm(&self) -> f64 {
74        (self.re * self.re + self.im * self.im).sqrt()
75    }
76
77    /// Squared magnitude re² + im².
78    pub fn norm_sqr(&self) -> f64 {
79        self.re * self.re + self.im * self.im
80    }
81
82    /// Phase angle atan2(im, re).
83    pub fn arg(&self) -> f64 {
84        self.im.atan2(self.re)
85    }
86
87    /// Complex conjugate (re, -im).
88    pub fn conj(&self) -> Self {
89        Self {
90            re: self.re,
91            im: -self.im,
92        }
93    }
94
95    /// Raise to a real power: (r e^{iθ})^p = r^p e^{ipθ}.
96    pub fn powf(&self, p: f64) -> Self {
97        let r = self.norm();
98        let theta = self.arg();
99        let rp = r.powf(p);
100        Self {
101            re: rp * (p * theta).cos(),
102            im: rp * (p * theta).sin(),
103        }
104    }
105
106    /// The zero complex number (0 + 0i).
107    pub fn zero() -> Self {
108        Self { re: 0.0, im: 0.0 }
109    }
110}
111
112impl std::ops::Add for C64 {
113    type Output = Self;
114    fn add(self, rhs: Self) -> Self {
115        Self {
116            re: self.re + rhs.re,
117            im: self.im + rhs.im,
118        }
119    }
120}
121
122impl std::ops::AddAssign for C64 {
123    fn add_assign(&mut self, rhs: Self) {
124        self.re += rhs.re;
125        self.im += rhs.im;
126    }
127}
128
129impl std::ops::Sub for C64 {
130    type Output = Self;
131    fn sub(self, rhs: Self) -> Self {
132        Self {
133            re: self.re - rhs.re,
134            im: self.im - rhs.im,
135        }
136    }
137}
138
139impl std::ops::Mul for C64 {
140    type Output = Self;
141    fn mul(self, rhs: Self) -> Self {
142        Self {
143            re: self.re * rhs.re - self.im * rhs.im,
144            im: self.re * rhs.im + self.im * rhs.re,
145        }
146    }
147}
148
149impl std::ops::Mul<f64> for C64 {
150    type Output = Self;
151    fn mul(self, rhs: f64) -> Self {
152        Self {
153            re: self.re * rhs,
154            im: self.im * rhs,
155        }
156    }
157}
158
159impl std::ops::Div for C64 {
160    type Output = Self;
161    fn div(self, rhs: Self) -> Self {
162        let denom = rhs.norm_sqr();
163        Self {
164            re: (self.re * rhs.re + self.im * rhs.im) / denom,
165            im: (self.im * rhs.re - self.re * rhs.im) / denom,
166        }
167    }
168}
169
170impl std::ops::Div<f64> for C64 {
171    type Output = Self;
172    fn div(self, rhs: f64) -> Self {
173        Self {
174            re: self.re / rhs,
175            im: self.im / rhs,
176        }
177    }
178}
179
180/// Result of a DMD computation.
181#[derive(Debug, Clone)]
182pub struct DmdResult {
183    /// Full Koopman operator approximation (m × m), complex.
184    pub a_matrix: Vec<Vec<C64>>,
185    /// DMD modes Φ (m × r), columns are modes.
186    pub modes: Vec<Vec<C64>>,
187    /// Eigenvalues λ (r).
188    pub eigenvalues: Vec<C64>,
189    /// Initial amplitudes b (r).
190    pub amplitudes: Vec<C64>,
191    /// Truncation rank used.
192    pub rank: usize,
193    /// Truncated SVD components.
194    pub svd: SvdComponents,
195    /// Reduced DMD matrix Ã (r × r), complex.
196    pub a_tilde: Vec<Vec<C64>>,
197    /// First snapshot (m).
198    pub x_first: Vec<f64>,
199    /// Last snapshot (m).
200    pub x_last: Vec<f64>,
201    /// Data dimensions (n_vars, n_time).
202    pub data_dim: (usize, usize),
203    /// Whether data was centered.
204    pub center: bool,
205    /// Row means (if centered).
206    pub x_mean: Option<Vec<f64>>,
207    /// Time step.
208    pub dt: f64,
209    /// Lifting metadata (if lifting was applied).
210    pub lifting_info: Option<LiftingInfo>,
211}
212
213impl DmdResult {
214    /// Number of original (pre-lifting) state variables.
215    pub fn n_vars_original(&self) -> usize {
216        match &self.lifting_info {
217            Some(info) => info.n_vars_original,
218            None => self.data_dim.0,
219        }
220    }
221
222    /// Whether lifting was applied.
223    pub fn is_lifted(&self) -> bool {
224        self.lifting_info.is_some()
225    }
226    /// Get mode column j as a slice of C64.
227    pub fn mode(&self, j: usize) -> Vec<C64> {
228        let n_vars = self.data_dim.0;
229        (0..n_vars).map(|i| self.modes[i][j]).collect()
230    }
231
232    /// Number of state variables.
233    pub fn n_vars(&self) -> usize {
234        self.data_dim.0
235    }
236}
237
238/// Information about a single DMD mode.
239#[derive(Debug, Clone)]
240pub struct ModeInfo {
241    /// Mode index.
242    pub index: usize,
243    /// Complex eigenvalue.
244    pub eigenvalue: C64,
245    /// Eigenvalue magnitude |λ|.
246    pub magnitude: f64,
247    /// Eigenvalue phase angle (radians).
248    pub phase: f64,
249    /// Oscillation frequency (cycles per dt).
250    pub frequency: f64,
251    /// Oscillation period (in dt units).
252    pub period: f64,
253    /// Growth rate (log|λ|/dt).
254    pub growth_rate: f64,
255    /// Half-life for decaying modes (positive), doubling time for growing (negative).
256    pub half_life: Option<f64>,
257    /// Stability classification.
258    pub stability: Stability,
259    /// Mode amplitude |b|.
260    pub amplitude: f64,
261}
262
263/// Stability classification of a mode or system.
264#[derive(Debug, Clone, Copy, PartialEq, Eq)]
265pub enum Stability {
266    Decaying,
267    Neutral,
268    Growing,
269}
270
271impl std::fmt::Display for Stability {
272    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
273        match self {
274            Stability::Decaying => write!(f, "decaying"),
275            Stability::Neutral => write!(f, "neutral"),
276            Stability::Growing => write!(f, "growing"),
277        }
278    }
279}
280
281/// Result of stability analysis.
282#[derive(Debug, Clone)]
283pub struct StabilityResult {
284    /// Whether all modes are decaying or neutral.
285    pub is_stable: bool,
286    /// Whether any mode is growing.
287    pub is_unstable: bool,
288    /// Whether any mode is exactly neutral (within tolerance).
289    pub is_marginal: bool,
290    /// Maximum eigenvalue magnitude.
291    pub spectral_radius: f64,
292    /// Per-mode stability classification.
293    pub mode_stability: Vec<Stability>,
294}
295
296/// Error metrics for reconstruction quality.
297#[derive(Debug, Clone)]
298pub struct ErrorMetrics {
299    /// Root mean square error.
300    pub rmse: f64,
301    /// Mean absolute error.
302    pub mae: f64,
303    /// Mean absolute percentage error.
304    pub mape: f64,
305    /// Relative error (Frobenius norm ratio).
306    pub relative_error: f64,
307    /// Per-variable RMSE.
308    pub per_variable_rmse: Vec<f64>,
309}
310
311/// Criterion for selecting dominant modes.
312#[derive(Debug, Clone, Copy)]
313pub enum DominantCriterion {
314    /// Sort by amplitude |b|.
315    Amplitude,
316    /// Sort by energy |b| × |λ|.
317    Energy,
318    /// Sort by stability (closest to unit circle first).
319    Stability,
320}
321
322/// Result of residual analysis.
323#[derive(Debug, Clone)]
324pub struct ResidualResult {
325    /// Overall residual Frobenius norm.
326    pub residual_norm: f64,
327    /// Relative residual (residual_norm / data_norm).
328    pub residual_relative: f64,
329    /// Per-step residual norms.
330    pub per_step_residual: Vec<f64>,
331    /// Per-mode residual contributions.
332    pub per_mode_residual: Vec<f64>,
333}
334
335/// Result of pseudospectrum computation.
336#[derive(Debug, Clone)]
337pub struct PseudospectrumResult {
338    /// Real axis grid points.
339    pub x: Vec<f64>,
340    /// Imaginary axis grid points.
341    pub y: Vec<f64>,
342    /// Minimum singular values at each grid point (grid_n × grid_n, row-major).
343    pub sigma_min: Vec<Vec<f64>>,
344    /// DMD eigenvalues for reference.
345    pub eigenvalues: Vec<C64>,
346    /// Epsilon contour levels.
347    pub epsilon: Vec<f64>,
348}
349
350/// Result of convergence analysis.
351#[derive(Debug, Clone)]
352pub struct ConvergenceResult {
353    /// Sample sizes used.
354    pub sample_sizes: Vec<usize>,
355    /// Eigenvalues at each sample size.
356    pub eigenvalues: Vec<Vec<C64>>,
357    /// Max eigenvalue magnitude changes between successive fits.
358    pub eigenvalue_changes: Vec<f64>,
359    /// Estimated convergence rate (O(1/m^alpha)), None if insufficient data.
360    pub convergence_estimate: Option<f64>,
361}
koopman_dmd/types.rs

koopman_dmd/
types.rs
