atelier_quant 0.0.12

//! Numerical quadrature for conditional-mean forecasting.
//!
//! The expected gap to the next event given current kernel state is:
//!
//! $$E[\Delta] = \int_0^\infty S(s)\,ds$$
//!
//! where the survival function is:
//!
//! $$S(s) = \exp\!\bigl(-\mu\,s - K(s)\bigr)$$
//!
//! and $K(s)$ is the integrated intensity contribution from the kernel
//! ([`ExcitationKernel::integrated_intensity_contribution`](super::kernel::ExcitationKernel::integrated_intensity_contribution)).
//!
//! We apply the substitution $x = \mu \cdot s$ so the integration
//! domain becomes $[0, x_{\max}]$ with $x_{\max} \approx 40$ (since
//! $e^{-40} \approx 4 \times 10^{-18}$), regardless of the scale of $\mu$.
//!
//! The integrand under the substitution is:
//!
//! $$\frac{1}{\mu}\,\exp\!\bigl(-x - K(x/\mu)\bigr)$$
//!
//! We use adaptive Simpson's rule for robust accuracy.

/// Adaptive Simpson's quadrature of `f` on `[a, b]`.
///
/// Subdivides until the estimated error is below `tol` or the recursion
/// depth exceeds `max_depth`.
pub fn adaptive_simpson<F: Fn(f64) -> f64>(
    f: &F,
    a: f64,
    b: f64,
    tol: f64,
    max_depth: usize,
) -> f64 {
    let mid = 0.5 * (a + b);
    let whole = simpson_rule(f, a, b);
    let left = simpson_rule(f, a, mid);
    let right = simpson_rule(f, mid, b);
    adaptive_step(f, a, b, tol, whole, left, right, max_depth)
}

fn simpson_rule<F: Fn(f64) -> f64>(f: &F, a: f64, b: f64) -> f64 {
    let h = b - a;
    (h / 6.0) * (f(a) + 4.0 * f(0.5 * (a + b)) + f(b))
}

fn adaptive_step<F: Fn(f64) -> f64>(
    f: &F,
    a: f64,
    b: f64,
    tol: f64,
    whole: f64,
    left: f64,
    right: f64,
    depth: usize,
) -> f64 {
    let sum = left + right;
    let err = (sum - whole) / 15.0;

    if depth == 0 || err.abs() < tol {
        return sum + err; // Richardson extrapolation correction
    }

    let mid = 0.5 * (a + b);
    let half_tol = 0.5 * tol;

    let ml = 0.5 * (a + mid);
    let mr = 0.5 * (mid + b);

    let left_left = simpson_rule(f, a, ml);
    let left_right = simpson_rule(f, ml, mid);
    let right_left = simpson_rule(f, mid, mr);
    let right_right = simpson_rule(f, mr, b);

    adaptive_step(f, a, mid, half_tol, left, left_left, left_right, depth - 1)
        + adaptive_step(
            f,
            mid,
            b,
            half_tol,
            right,
            right_left,
            right_right,
            depth - 1,
        )
}

/// Compute the conditional expected gap $E[\Delta \mid \text{state}]$.
///
/// Uses the substitution $x = \mu \cdot s$ so the integral becomes:
///
/// $$E[\Delta] = \frac{1}{\mu} \int_0^{x_{\max}} \exp\!\bigl(-x - K(x/\mu)\bigr)\,dx$$
///
/// where $K(s)$ is provided by `kernel_integrated_intensity`, a closure
/// that computes `kernel.integrated_intensity_contribution(state, s)`.
///
/// # Arguments
///
/// * `mu` — baseline intensity.
/// * `kernel_integrated_intensity` — closure `|s| -> f64` returning the
///   integrated kernel contribution from the current state over interval $s$.
/// * `tol` — absolute tolerance for adaptive quadrature (default: 1e-10).
///
/// # Returns
///
/// The conditional expected gap in the same time units as $\mu$.
pub fn conditional_mean_gap<F: Fn(f64) -> f64>(
    mu: f64,
    kernel_integrated_intensity: F,
    tol: f64,
) -> f64 {
    if mu <= 0.0 || !mu.is_finite() {
        return f64::INFINITY;
    }

    // Integration in the substituted variable x = μ·s, domain [0, x_max].
    // S(x/μ) = exp(-x - K(x/μ)),  integrand = (1/μ) · S(x/μ)
    //
    // x_max = 40 gives exp(-40) ≈ 4e-18, well below machine epsilon
    // for the survival function even before the kernel contribution.
    let x_max = 40.0_f64;

    let inv_mu = 1.0 / mu;

    let integrand = |x: f64| -> f64 {
        let s = x * inv_mu;
        let k = kernel_integrated_intensity(s);
        (-x - k).exp()
    };

    // The (1/μ) factor comes from the ds = dx/μ substitution.
    inv_mu * adaptive_simpson(&integrand, 0.0, x_max, tol, 20)
}