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}
gam_geometry/integrator.rs

gam_geometry/
integrator.rs
