gam_geometry/
integrator.rs

1use ndarray::{Array1, Array2, ArrayView1};
2
3use crate::manifold::{GeometryResult, RiemannianManifold, check_len};
4
5#[derive(Debug, Clone, Copy, PartialEq)]
6pub struct GeodesicIntegrator {
7    pub steps: usize,
8    pub step_size: f64,
9}
10
11impl Default for GeodesicIntegrator {
12    fn default() -> Self {
13        Self {
14            steps: 32,
15            step_size: 1.0 / 32.0,
16        }
17    }
18}
19
20impl GeodesicIntegrator {
21    /// Integrate the geodesic ODE on the manifold by composing the closed-form
22    /// geodesic flow `exp_p(v · h)` with parallel transport of the velocity
23    /// along that segment. This is the mathematically correct flow on every
24    /// Riemannian manifold whose `RiemannianManifold` impl supplies an
25    /// `exp_map` and `parallel_transport` consistent with its metric — it
26    /// reduces to a Christoffel-driven leap-frog integrator on coordinates
27    /// where Γ ≡ 0 (e.g. Euclidean / flat tori) but stays on the manifold
28    /// and conserves the Riemannian kinetic energy on positively curved
29    /// spaces such as the sphere, which previously drifted off because the
30    /// extrinsic Christoffel symbols are not zero.
31    pub fn integrate(
32        &self,
33        manifold: &dyn RiemannianManifold,
34        point: ArrayView1<'_, f64>,
35        tangent: ArrayView1<'_, f64>,
36    ) -> GeometryResult<Array1<f64>> {
37        let d = manifold.ambient_dim();
38        check_len("geodesic integrator point", point.len(), d)?;
39        check_len("geodesic integrator tangent", tangent.len(), d)?;
40        // Integrate the geodesic from parameter t=0 to t=1 so the result is
41        // the canonical unit-time endpoint exp_p(v). The IVP energy ½|v(t)|²
42        // is invariant along any geodesic; the BVP quantity ½|log_p(γ(1))|²
43        // matches ½|v|² only when total integration time equals 1. We honour
44        // both `steps` (substep count) and `step_size` (substep duration) by
45        // taking the finer of the two — `n = max(steps, ⌈1/step_size⌉)` —
46        // each of duration `h = 1/n`. This gives strictly more accuracy when
47        // either bound is tightened and never overshoots the unit-time
48        // endpoint, so closed-form exp_map/parallel_transport composition is
49        // numerically exact for any caller-provided refinement.
50        let mut x = point.to_owned();
51        let mut v = manifold.project_tangent(point, tangent)?;
52        let steps_from_size = if self.step_size.is_finite() && self.step_size > 0.0 {
53            (1.0 / self.step_size).ceil() as usize
54        } else {
55            1
56        };
57        let n_substeps = self.steps.max(1).max(steps_from_size);
58        let h = 1.0 / n_substeps as f64;
59        for _ in 0..n_substeps {
60            let step_tangent = &v * h;
61            let x_next = manifold.exp_map(x.view(), step_tangent.view())?;
62            let mut path = Array2::<f64>::zeros((2, d));
63            path.row_mut(0).assign(&x);
64            path.row_mut(1).assign(&x_next);
65            let v_next = manifold.parallel_transport(path.view(), v.view())?;
66            x = x_next;
67            v = manifold.project_tangent(x.view(), v_next.view())?;
68        }
69        Ok(x)
70    }
71}
72
73#[cfg(test)]
74mod tests {
75    use super::*;
76    use crate::manifolds::euclidean::EuclideanManifold;
77    use ndarray::array;
78
79    fn integrator_default() -> GeodesicIntegrator {
80        GeodesicIntegrator::default()
81    }
82
83    #[test]
84    fn default_has_expected_steps_and_step_size() {
85        let g = integrator_default();
86        assert_eq!(g.steps, 32);
87        assert!((g.step_size - 1.0 / 32.0).abs() < 1e-15);
88    }
89
90    #[test]
91    fn euclidean_geodesic_is_straight_line() {
92        // On R^n the geodesic is a straight line: exp_p(v) = p + v.
93        // The integrator should recover p + v exactly regardless of step count.
94        let m = EuclideanManifold::new(2);
95        let p = array![1.0_f64, 2.0];
96        let v = array![3.0_f64, 4.0];
97        let result = integrator_default()
98            .integrate(&m, p.view(), v.view())
99            .unwrap();
100        assert!((result[0] - 4.0).abs() < 1e-12, "x: {}", result[0]);
101        assert!((result[1] - 6.0).abs() < 1e-12, "y: {}", result[1]);
102    }
103
104    #[test]
105    fn euclidean_zero_tangent_returns_same_point() {
106        let m = EuclideanManifold::new(3);
107        let p = array![5.0_f64, -1.0, 2.0];
108        let v = array![0.0_f64, 0.0, 0.0];
109        let result = integrator_default()
110            .integrate(&m, p.view(), v.view())
111            .unwrap();
112        for i in 0..3 {
113            assert!((result[i] - p[i]).abs() < 1e-12);
114        }
115    }
116
117    #[test]
118    fn dimension_mismatch_point_returns_error() {
119        let m = EuclideanManifold::new(3);
120        let p = array![1.0_f64, 2.0]; // wrong length
121        let v = array![0.0_f64, 0.0, 0.0];
122        assert!(integrator_default().integrate(&m, p.view(), v.view()).is_err());
123    }
124
125    #[test]
126    fn single_step_integrator_still_recovers_euclidean_geodesic() {
127        let m = EuclideanManifold::new(2);
128        let g = GeodesicIntegrator { steps: 1, step_size: 1.0 };
129        let p = array![0.0_f64, 0.0];
130        let v = array![2.0_f64, -3.0];
131        let result = g.integrate(&m, p.view(), v.view()).unwrap();
132        assert!((result[0] - 2.0).abs() < 1e-12);
133        assert!((result[1] - (-3.0)).abs() < 1e-12);
134    }
135}
gam_geometry/integrator.rs

gam_geometry/
integrator.rs
