Skip to main content

cartan_stochastic/
development.rs

1//! Stochastic development (Eells-Elworthy-Malliavin).
2//!
3//! The rolling-without-slipping construction: drive `O(M)` along the
4//! horizontal lift of Euclidean Brownian motion `W_t ∈ R^n`. The base-point
5//! projection is a diffusion on `M` with generator `½ Δ_M` (Laplace-Beltrami).
6//!
7//! Anti-development is the inverse: given a curve on `M`, read off its
8//! representation in any chosen initial frame via repeated parallel-transport
9//! and inner products. That is not implemented here because it is an
10//! application rather than a primitive — call `Manifold::log` to get
11//! consecutive tangent vectors, then `Manifold::inner` against each frame
12//! basis vector.
13
14use cartan_core::Real;
15use rand::Rng;
16use rand_distr::{Distribution, StandardNormal};
17
18use crate::error::StochasticError;
19use crate::frame::OrthonormalFrame;
20use crate::sde::{stratonovich_step, StratonovichDevelopment};
21
22/// A stored stochastic-development trajectory on `M`.
23///
24/// Contains the base-point path and the final frame. Intermediate frames
25/// are not stored by default because most downstream consumers (BEL weight
26/// computation, Bismut integration-by-parts) only need the base-point path
27/// and a terminal frame. Callers needing the full frame history can call
28/// `stratonovich_step` directly in a loop.
29#[derive(Debug, Clone)]
30pub struct DevelopmentPath<M: cartan_core::Manifold> {
31    /// The base-point path, length `n_steps + 1`, starting at the caller's
32    /// initial point.
33    pub path: Vec<M::Point>,
34    /// The orthonormal frame at the terminal point, parallel-transported and
35    /// re-orthonormalised at each step.
36    pub final_frame: OrthonormalFrame<M>,
37}
38
39/// Integrate a Brownian motion on `M` via stochastic development.
40///
41/// Draws `n_steps × dim(M)` i.i.d. standard normals and steps the Stratonovich
42/// SDE on the orthonormal frame bundle. Returns the full base-point path.
43///
44/// `dt` is the per-step time increment; standard practice is `T / n_steps`
45/// for total horizon `T`. `tol` is the Gram-Schmidt rank-deficiency threshold
46/// used to re-orthonormalise the frame after each step.
47pub fn stochastic_development<M: StratonovichDevelopment, R: Rng + ?Sized>(
48    manifold: &M,
49    p0: &M::Point,
50    frame0: OrthonormalFrame<M>,
51    n_steps: usize,
52    dt: Real,
53    rng: &mut R,
54    tol: Real,
55) -> Result<DevelopmentPath<M>, StochasticError>
56where
57    StandardNormal: Distribution<Real>,
58{
59    let n = manifold.dim();
60    let mut path: Vec<M::Point> = Vec::with_capacity(n_steps + 1);
61    path.push(p0.clone());
62    let mut p = p0.clone();
63    let mut frame = frame0;
64    let mut dw = vec![0.0_f64; n];
65    for _ in 0..n_steps {
66        for w in dw.iter_mut() {
67            *w = StandardNormal.sample(rng);
68        }
69        let (p_next, frame_next) = stratonovich_step(manifold, &p, &frame, &dw, dt, tol)?;
70        path.push(p_next.clone());
71        p = p_next;
72        frame = frame_next;
73    }
74    Ok(DevelopmentPath { path, final_frame: frame })
75}