phop_core/
ode.rs

1//! Governing-equation discovery: recover the right-hand side of an autonomous ODE
2//! `dx/dt = f(x)` from a sampled trajectory, using phop as the function learner.
3//!
4//! This is the "phop for dynamical systems" research direction: estimate the derivative `dx/dt`
5//! from a uniformly-sampled scalar series by central differences, then run the standard discovery
6//! pipeline on the regression `x ↦ dx/dt`. The recovered law *is* the governing equation. (It is a
7//! lightweight, self-contained cousin of SINDy — phop supplies the nonlinear function class instead
8//! of a fixed sparse library; coupling to `oxieml::symreg::{sindy, pde}` for multi-variable systems
9//! is future work.)
10
11use crate::config::Config;
12use crate::dataset::DataSet;
13use crate::discoverer::Discoverer;
14use crate::error::{PhopError, Result};
15use crate::pareto::ParetoFront;
16use scirs2_core::ndarray::{Array1, Array2};
17
18/// Discover the right-hand side `f` of an autonomous scalar ODE `dx/dt = f(x)` from a
19/// uniformly-sampled trajectory `series` with timestep `dt`.
20///
21/// Derivatives at interior samples are estimated with the central difference
22/// `(x[i+1] − x[i−1]) / (2·dt)`; the resulting `(x, dx/dt)` pairs are handed to [`Discoverer::fit`].
23/// The returned Pareto front is over candidate `f`.
24///
25/// # Errors
26/// Returns [`PhopError::ShapeMismatch`] if the series is too short (fewer than 3 samples) or `dt`
27/// is non-positive, or propagates any discovery error.
28pub fn discover_ode(series: &[f64], dt: f64, cfg: &Config) -> Result<ParetoFront> {
29    if series.len() < 3 {
30        return Err(PhopError::ShapeMismatch(
31            "need at least 3 samples to estimate a derivative".to_string(),
32        ));
33    }
34    if dt <= 0.0 || dt.is_nan() {
35        return Err(PhopError::ShapeMismatch("dt must be positive".to_string()));
36    }
37    let mut xs = Vec::with_capacity(series.len() - 2);
38    let mut dxdt = Vec::with_capacity(series.len() - 2);
39    for i in 1..series.len() - 1 {
40        xs.push(series[i]);
41        dxdt.push((series[i + 1] - series[i - 1]) / (2.0 * dt));
42    }
43    let x = Array2::from_shape_vec((xs.len(), 1), xs)
44        .map_err(|e| PhopError::ShapeMismatch(e.to_string()))?;
45    let ds = DataSet::from_arrays(x, Array1::from(dxdt))?;
46    Discoverer::new(cfg.clone()).fit(&ds)
47}
48
49#[cfg(test)]
50mod tests {
51    use super::*;
52
53    #[test]
54    fn recovers_exponential_growth_ode() {
55        // x(t) = e^t solves dx/dt = x. Sampling e^t and recovering f should give f(x) = x.
56        let dt = 0.02;
57        let series: Vec<f64> = (0..200).map(|i| (i as f64 * dt).exp()).collect();
58        let cfg = Config::default().max_depth(2).max_epochs(200).seed(0);
59        let front = discover_ode(&series, dt, &cfg).expect("ode discovery");
60        let best = front.best().expect("non-empty front");
61        // dx/dt = x is the identity in x; central-difference error is O(dt^2), so the relative fit
62        // is tight. Use a relative-MSE check against the variance of dx/dt (= variance of x).
63        let mean = series.iter().sum::<f64>() / series.len() as f64;
64        let var = series.iter().map(|v| (v - mean) * (v - mean)).sum::<f64>() / series.len() as f64;
65        assert!(
66            best.mse < var * 1e-3,
67            "ODE rhs not recovered: mse {} vs var {} ({})",
68            best.mse,
69            var,
70            best.pretty()
71        );
72    }
73}
phop_core/ode.rs

phop_core/
ode.rs
