modewise_mvmd/
params.rs

1use modewise_method::MethodError;
2
3/// Starting centre frequencies, in normalised frequency (cycles per sample).
4#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
5pub enum FrequencyInit {
6    #[default]
7    Zero,
8    /// Evenly spaced over [0, 0.5].
9    Linear,
10    /// Log-spaced over [0.0005, 0.5].
11    Exponential,
12    /// (0.5 / K)·k, as the original MATLAB VMD and MVMD codes and vmdpy (`init = 1`).
13    Uniform,
14}
15
16impl FrequencyInit {
17    pub(crate) fn omegas(self, num_modes: usize) -> Vec<f64> {
18        let step = |i: usize| i as f64 / (num_modes - 1).max(1) as f64;
19        (0..num_modes)
20            .map(|i| match self {
21                Self::Zero => 0.0,
22                Self::Linear => 0.5 * step(i),
23                Self::Exponential => 0.5 * 10f64.powf(-3.0 + 3.0 * step(i)),
24                Self::Uniform => 0.5 * i as f64 / num_modes as f64,
25            })
26            .collect()
27    }
28}
29
30/// Defaults follow the reference implementation.
31#[derive(Debug, Clone, PartialEq)]
32pub struct MvmdParams {
33    pub num_modes: usize,
34    /// Bandwidth penalty. Larger values give narrower modes.
35    pub alpha: f64,
36    pub init: FrequencyInit,
37    /// Stop when the mean squared mode update falls below this.
38    pub tolerance: f64,
39    /// Dual ascent step. 0 allows noise to stay out of the modes.
40    pub tau: f64,
41    pub max_iterations: usize,
42}
43
44impl MvmdParams {
45    pub fn new(num_modes: usize, alpha: f64) -> Self {
46        Self {
47            num_modes,
48            alpha,
49            init: FrequencyInit::Zero,
50            tolerance: 1e-3,
51            tau: 1e-2,
52            max_iterations: 500,
53        }
54    }
55
56    pub(crate) fn validate(&self) -> Result<(), MethodError> {
57        if self.num_modes == 0 {
58            return Err(MethodError::invalid("num_modes", "must be at least 1"));
59        }
60        if !(self.alpha.is_finite() && self.alpha > 0.0) {
61            return Err(MethodError::invalid("alpha", "must be finite and positive"));
62        }
63        if !(self.tolerance.is_finite() && self.tolerance > 0.0) {
64            return Err(MethodError::invalid(
65                "tolerance",
66                "must be finite and positive",
67            ));
68        }
69        if !(self.tau.is_finite() && self.tau >= 0.0) {
70            return Err(MethodError::invalid(
71                "tau",
72                "must be finite and non-negative",
73            ));
74        }
75        if self.max_iterations == 0 {
76            return Err(MethodError::invalid("max_iterations", "must be at least 1"));
77        }
78        Ok(())
79    }
80}
81
82#[cfg(test)]
83mod tests {
84    use super::*;
85
86    #[test]
87    fn init_matches_numpy_spacing() {
88        assert_eq!(FrequencyInit::Linear.omegas(3), [0.0, 0.25, 0.5]);
89        let exp = FrequencyInit::Exponential.omegas(4);
90        let expected = [0.0005, 0.005, 0.05, 0.5];
91        for (a, b) in exp.iter().zip(expected) {
92            assert!((a - b).abs() < 1e-15, "{a} != {b}");
93        }
94        assert_eq!(FrequencyInit::Exponential.omegas(1), [0.0005]);
95    }
96
97    #[test]
98    fn rejects_invalid_values() {
99        let base = MvmdParams::new(3, 1000.0);
100        let cases = [
101            MvmdParams {
102                num_modes: 0,
103                ..base.clone()
104            },
105            MvmdParams {
106                alpha: 0.0,
107                ..base.clone()
108            },
109            MvmdParams {
110                alpha: f64::NAN,
111                ..base.clone()
112            },
113            MvmdParams {
114                tolerance: 0.0,
115                ..base.clone()
116            },
117            MvmdParams {
118                tau: -1.0,
119                ..base.clone()
120            },
121            MvmdParams {
122                max_iterations: 0,
123                ..base.clone()
124            },
125        ];
126        for params in cases {
127            assert!(params.validate().is_err(), "accepted {params:?}");
128        }
129        assert!(base.validate().is_ok());
130    }
131}
modewise_mvmd/params.rs

modewise_mvmd/
params.rs
