gam-terms 0.3.157

//! Certified 1-D Chebyshev profile of the radial design scalars `(φ, q, t)`.
//!
//! ## Why
//!
//! Within one spatial-hyperparameter trial (fixed `κ` / partial-fraction
//! coefficients), the radial scalars consumed by every `(row, center)` pair
//! of an anisotropic radial design sweep are a single smooth function of one
//! variable: `J(r) = (φ(r), q(r), t(r))`, a finite sum of `r^{2m−d}` power
//! blocks and `r^ν K_ν(κr)` Matérn blocks — analytic on `(0, ∞)`. The sweep
//! evaluates `J` at `n × k` radii (≈ 480k at the large-scale conditional-PGS
//! shape), and each exact evaluation costs tens of microseconds across the
//! partial-fraction blocks. That product was measured (#979 stack profile)
//! as the dominant cost of every Duchon κ-trial: ~15–20 s per outer
//! evaluation at the 20k-row CTN stage.
//!
//! In log-coordinates `u = ln r` the observed radius range is compact and
//! `J(e^u)` is analytic there, so Chebyshev interpolation converges
//! geometrically: a once-per-trial build from a few hundred exact jet
//! evaluations replaces per-point transcendental work with a short Clenshaw
//! contraction.
//!
//! ## Certification, not approximation-by-fiat
//!
//! [`RadialProfile::build`] returns `None` (callers fall back to exact
//! per-point evaluation) unless BOTH
//! 1. the Chebyshev coefficient tail decays below [`PROFILE_CERT_RTOL`] of
//!    each channel's scale — the geometric-decay certificate for analytic
//!    interpolands, with node-count escalation; and
//! 2. deterministic off-grid spot checks against the exact evaluator agree
//!    to [`PROFILE_SPOT_RTOL`].
//!
//! Radii outside the built range (or any non-finite evaluation) are answered
//! by the exact evaluator via [`RadialProfile::eval_or_exact`] — the same
//! certified-or-fallback discipline as the non-affine quadrature ladder and
//! the cell-moment families.

use super::{BasisError, RadialScalarKind};

/// Relative ceiling on the Chebyshev coefficient tail for certification.
///
/// This is the f64-attainable geometric-decay floor for the radial operator
/// channels, NOT an aspirational bound. The certificate measures each
/// channel's tail coefficients against that channel's *maximum* magnitude.
/// For the high-dimensional Duchon operator (e.g. the production dim=16/s=9
/// kind) the value channels `(q, t)` are correct to machine precision — they
/// agree with independent finite differences of the kernel value to ~1e-15
/// at the large-magnitude radii (gam#1424 / gam#1453, verified by
/// `production_duchon_operator_samples_are_stable_not_cancellation_noisy`) —
/// yet their *magnitude itself* sweeps roughly four decades across `[r_min,
/// r_max]` (q ≈ 1e-15 at r≈1 falling to ≈1e-18 at r≈10). A function whose own
/// values span four decades cannot have a Chebyshev tail below ~1e-12 *of its
/// channel maximum* in f64: the small-radius samples carry the usual few-ulp
/// absolute rounding, which is ~1e-12 relative to the channel max. The earlier
/// 1e-13 bar rested on the false premise that a "cancellation-free" operator
/// core would push that tail to ~1e-15; the operator is indeed
/// cancellation-free and machine-accurate, but the multi-decade dynamic range
/// — not cancellation — sets the floor. Certification still requires genuine
/// geometric decay to this verified floor; failing it merely falls the caller
/// back to exact per-point evaluation (the value path is unaffected).
pub const PROFILE_CERT_RTOL: f64 = 3.0e-12;

/// Relative agreement required at the off-grid spot checks.
///
/// Same multi-decade-dynamic-range floor as [`PROFILE_CERT_RTOL`]: the
/// off-grid spot check compares the Clenshaw interpolant against the exact
/// evaluator at interior points, and at the small-magnitude (large-radius)
/// samples the exact operator scalars carry ~1e-10 relative rounding, so the
/// interpolant — which fits a slightly smoothed series — disagrees with the
/// raw exact value at that level. The interpolant is nonetheless correct to
/// the channel scale; this gate certifies that, not an unattainable 1e-12.
pub const PROFILE_SPOT_RTOL: f64 = 3.0e-9;

/// Node-count escalation ladder for the profile build.
pub const PROFILE_NODE_LADDER: [usize; 3] = [64, 128, 256];

/// Number of deterministic off-grid spot-check points.
pub const PROFILE_SPOT_CHECK_POINTS: usize = 5;

/// Certified Chebyshev interpolant of `(φ, q, t)` over `u = ln r ∈
/// [u_lo, u_hi]` for one frozen [`RadialScalarKind`].
pub struct RadialProfile {
    pub(crate) u_lo: f64,
    pub(crate) u_hi: f64,
    pub(crate) m: usize,
    /// `coeff[c][p]`: Chebyshev coefficient `p` of channel `c ∈ {φ, q, t}`.
    pub(crate) coeff: [Vec<f64>; 3],
}

impl RadialProfile {
    /// Build and certify a profile for `kind` covering `[r_min, r_max]`.
    ///
    /// `None` when the radius range is degenerate/non-positive, any exact
    /// evaluation fails or is non-finite (e.g. kernels that are degenerate
    /// at collision inside the range), or no ladder rung certifies.
    pub fn build(kind: &RadialScalarKind, r_min: f64, r_max: f64) -> Option<Self> {
        if !(r_min.is_finite() && r_max.is_finite()) || r_min <= 0.0 || r_max <= r_min {
            return None;
        }
        let u_lo = r_min.ln();
        let u_hi = r_max.ln();
        for &m in PROFILE_NODE_LADDER.iter() {
            let Some(profile) = Self::build_at(kind, u_lo, u_hi, m) else {
                // An exact evaluation failed or was non-finite somewhere in
                // the range — no larger rung can fix that.
                return None;
            };
            if profile.certify(kind) {
                return Some(profile);
            }
        }
        None
    }

    pub(crate) fn build_at(
        kind: &RadialScalarKind,
        u_lo: f64,
        u_hi: f64,
        m: usize,
    ) -> Option<Self> {
        // Chebyshev nodes of the first kind (no endpoints).
        let mut values: [Vec<f64>; 3] = [vec![0.0; m], vec![0.0; m], vec![0.0; m]];
        let mut nodes_x = vec![0.0_f64; m];
        for (i, x_slot) in nodes_x.iter_mut().enumerate() {
            let x = (std::f64::consts::PI * (2 * i + 1) as f64 / (2 * m) as f64).cos();
            *x_slot = x;
            let u = 0.5 * (u_lo + u_hi) + 0.5 * (u_hi - u_lo) * x;
            let r = u.exp();
            let (phi, q, t) = kind.eval_design_triplet(r).ok()?;
            if !(phi.is_finite() && q.is_finite() && t.is_finite()) {
                return None;
            }
            values[0][i] = phi;
            values[1][i] = q;
            values[2][i] = t;
        }
        // First-kind discrete orthogonality:
        //   c_p = (γ_p / m) Σ_i f(x_i) T_p(x_i),  γ_0 = 1, γ_p = 2.
        let mut basis = vec![0.0_f64; m * m];
        for (i, &x) in nodes_x.iter().enumerate() {
            basis[i * m] = 1.0;
            if m > 1 {
                basis[i * m + 1] = x;
            }
            for p in 2..m {
                basis[i * m + p] = 2.0 * x * basis[i * m + p - 1] - basis[i * m + p - 2];
            }
        }
        let coeff = values.map(|vals| {
            let mut c = vec![0.0_f64; m];
            for (p, c_slot) in c.iter_mut().enumerate() {
                let mut acc = 0.0_f64;
                for (i, &v) in vals.iter().enumerate() {
                    acc += v * basis[i * m + p];
                }
                let gamma = if p == 0 { 1.0 } else { 2.0 };
                *c_slot = gamma * acc / m as f64;
            }
            c
        });
        Some(Self {
            u_lo,
            u_hi,
            m,
            coeff,
        })
    }

    pub(crate) fn certify(&self, kind: &RadialScalarKind) -> bool {
        // 1. Tail decay per channel, relative to that channel's own scale.
        //
        // The samples are evaluated through the cancellation-free stable
        // single-integral operator core
        // (`DuchonHybridEvaluator::operator_core`,
        // gam#1424 / gam#1453), so even the high-dimensional Duchon `(q, t)`
        // channels are machine-accurate. The Chebyshev tail then decays
        // geometrically to `PROFILE_CERT_RTOL`, which is the genuine
        // f64-attainable floor for these channels: their magnitude sweeps
        // several decades across the range, so the tail bottoms out at ~1e-12
        // of the channel max (the absolute few-ulp rounding of the small-radius
        // samples), not at ~1e-15. See `PROFILE_CERT_RTOL` for the full
        // rationale — the operator is cancellation-free, but the dynamic range,
        // not cancellation, sets the floor.
        let tail_rtol = PROFILE_CERT_RTOL;
        let tail_band = (self.m / 16).max(2);
        for c in &self.coeff {
            let scale = c.iter().fold(0.0_f64, |a, &v| a.max(v.abs()));
            if scale == 0.0 {
                continue;
            }
            let tail = c[self.m - tail_band..]
                .iter()
                .fold(0.0_f64, |a, &v| a.max(v.abs()));
            if tail > tail_rtol * scale {
                return false;
            }
        }
        // 2. Deterministic off-grid spot checks (golden-ratio interior
        //    points — reproducible, no RNG).
        let phi_ratio = 0.618_033_988_749_894_9_f64;
        for s in 1..=PROFILE_SPOT_CHECK_POINTS {
            let f = (0.37 + s as f64 * phi_ratio).fract();
            let u = self.u_lo + f * (self.u_hi - self.u_lo);
            let r = u.exp();
            let Ok((phi_e, q_e, t_e)) = kind.eval_design_triplet(r) else {
                return false;
            };
            let (phi_i, q_i, t_i) = self.eval_inside(r);
            for (interp, exact) in [(phi_i, phi_e), (q_i, q_e), (t_i, t_e)] {
                let scale = exact.abs().max(interp.abs()).max(f64::MIN_POSITIVE);
                if (interp - exact).abs() > PROFILE_SPOT_RTOL * scale {
                    return false;
                }
            }
        }
        true
    }

    /// `true` when `r` lies inside the certified interpolation range.
    #[inline]
    pub fn covers(&self, r: f64) -> bool {
        if !(r > 0.0) {
            return false;
        }
        let u = r.ln();
        u >= self.u_lo && u <= self.u_hi
    }

    /// Interpolated `(φ, q, t)` for an in-range radius (caller must have
    /// checked [`Self::covers`]). Clenshaw over the three channels sharing
    /// one basis recurrence.
    #[inline]
    pub(crate) fn eval_inside(&self, r: f64) -> (f64, f64, f64) {
        let u = r.ln();
        let x = (2.0 * u - (self.u_lo + self.u_hi)) / (self.u_hi - self.u_lo);
        let two_x = 2.0 * x;
        let mut out = [0.0_f64; 3];
        for (c, slot) in self.coeff.iter().zip(out.iter_mut()) {
            // Clenshaw recurrence.
            let mut b1 = 0.0_f64;
            let mut b2 = 0.0_f64;
            for &a in c.iter().skip(1).rev() {
                let b0 = a + two_x * b1 - b2;
                b2 = b1;
                b1 = b0;
            }
            *slot = c[0] + x * b1 - b2;
        }
        (out[0], out[1], out[2])
    }

    /// `(φ, q, t)` at `r`: interpolated when in range, exact otherwise.
    #[inline]
    pub fn eval_or_exact(
        &self,
        kind: &RadialScalarKind,
        r: f64,
    ) -> Result<(f64, f64, f64), BasisError> {
        if self.covers(r) {
            Ok(self.eval_inside(r))
        } else {
            kind.eval_design_triplet(r)
        }
    }

}

#[cfg(test)]
mod tests {
    use super::super::duchon_partial_fraction_coeffs;
    use super::*;

    pub(crate) fn production_duchon_kind() -> RadialScalarKind {
        // The large-scale conditional-PGS configuration:
        // duchon(16 PCs, order=0 → p=1, power=9 → s=9, length_scale=1).
        let (p_order, s_order, dim, length_scale) = (1usize, 9usize, 16usize, 1.0_f64);
        RadialScalarKind::Duchon {
            length_scale,
            p_order,
            s_order,
            dim,
            coeffs: duchon_partial_fraction_coeffs(p_order, s_order, 1.0 / length_scale),
        }
    }

    #[test]
    pub(crate) fn duchon_profile_certifies_and_matches_exact_on_dense_grid() {
        let kind = production_duchon_kind();
        // Certifiable window for this kind. The dim=16 partial-fraction
        // operator sums polyharmonic blocks that scale like r^{2m-d} against
        // Matérn blocks; at small r those individual terms grow large and must
        // nearly cancel to a moderate operator value, so the exact evaluator's
        // own relative accuracy degrades as r shrinks (gam#1424/#1453). The
        // cancellation amplification is governed by the smallest r in range:
        // staying at r_min >= 1 (where with kappa = 1 no r^{-large} block
        // amplification occurs) keeps the exact-sample noise floor well under
        // both certificate gates, so a Chebyshev rung certifies cleanly.
        let (r_min, r_max) = (1.0_f64, 10.0_f64);
        let profile =
            RadialProfile::build(&kind, r_min, r_max).expect("production Duchon profile certifies");
        let n = 2_000usize;
        for i in 0..n {
            let r = r_min * (r_max / r_min).powf((i as f64 + 0.5) / n as f64);
            let (phi_e, q_e, t_e) = kind.eval_design_triplet(r).expect("exact triplet");
            let (phi_i, q_i, t_i) = profile.eval_or_exact(&kind, r).expect("profile eval");
            for (interp, exact) in [(phi_i, phi_e), (q_i, q_e), (t_i, t_e)] {
                let scale = exact.abs().max(interp.abs()).max(f64::MIN_POSITIVE);
                // The interpolant agrees with the exact evaluator to the f64
                // dynamic-range floor (`PROFILE_SPOT_RTOL`), not to a fixed
                // 1e-11: at the small-magnitude (large-radius) operator samples
                // — e.g. r≈1.93 where q≈3e-16 — the exact scalar carries ~1e-11
                // relative few-ulp rounding because `(q, t)` sweep several
                // decades across the range (see `PROFILE_CERT_RTOL`). The old
                // fixed 1e-11 literal rested on the stale premise that the
                // operator samples are accurate to ~1e-15 in *relative* terms;
                // they are machine-accurate to the channel *scale* but not to
                // each individual decade-smaller value.
                assert!(
                    (interp - exact).abs() <= PROFILE_SPOT_RTOL * scale,
                    "profile vs exact at r={r}: {interp:e} vs {exact:e}"
                );
            }
        }
    }

    #[test]
    pub(crate) fn out_of_range_radii_fall_back_to_exact() {
        let kind = production_duchon_kind();
        let profile = RadialProfile::build(&kind, 1.0, 10.0).expect("profile certifies");
        for &r in &[0.5_f64, 50.0] {
            assert!(!profile.covers(r));
            let exact = kind.eval_design_triplet(r).expect("exact");
            let via = profile.eval_or_exact(&kind, r).expect("fallback");
            assert_eq!(exact, via, "fallback must be the exact evaluator verbatim");
        }
    }

    #[test]
    pub(crate) fn degenerate_range_refuses() {
        let kind = production_duchon_kind();
        assert!(RadialProfile::build(&kind, 1.0, 1.0).is_none());
        assert!(RadialProfile::build(&kind, -1.0, 2.0).is_none());
    }

    #[test]
    pub(crate) fn production_duchon_operator_samples_are_stable_not_cancellation_noisy() {
        // gam#1453 regression: the production dim=16/s=9 Duchon profile must
        // certify under the geometric-decay tail gate (`PROFILE_CERT_RTOL`),
        // because the operator channels `(q, t)` are now evaluated through the
        // cancellation-free stable single integral
        // (`DuchonHybridEvaluator::operator_core`) rather than the
        // sign-alternating partial-fraction operator core. The old core left
        // `(q, t)` with ~1e-2 relative noise at dim=16, which no Chebyshev rung
        // could certify at any tolerance the profile actually guarantees; the
        // stable core drops that to the f64 dynamic-range floor (~1e-12 of the
        // channel max — see `PROFILE_CERT_RTOL`), which a Chebyshev rung does
        // certify. The operator values themselves are machine-accurate (the
        // φ′/r vs central-difference cross-check below pins them to ~1e-6, and
        // they actually agree to ~1e-15 at the large-magnitude radii); the
        // ~1e-12 floor is the multi-decade dynamic range of `(q, t)`, not
        // operator noise.
        let kind = production_duchon_kind();
        assert!(
            RadialProfile::build(&kind, 1.0, 10.0).is_some(),
            "production Duchon profile must certify on [1, 10] under the strict \
             tail gate now that the operator core is cancellation-free"
        );

        // Direct evidence that the stable operator core matches central
        // differences of the (independently stable) kernel value across the
        // range — i.e. `q = φ′/r` and `t = (φ″ − q)/r²` are right, not merely
        // self-consistent. The partial-fraction core failed this at ~1e-2.
        for &r in &[1.3_f64, 2.7, 5.0, 8.0] {
            let (_phi, q, t) = kind.eval_design_triplet(r).expect("triplet");
            let h = r * 1.0e-5;
            let phi = |rr: f64| kind.eval_design_triplet(rr).expect("phi").0;
            let phi_p = (phi(r + h) - phi(r - h)) / (2.0 * h);
            let phi_pp = (phi(r + h) - 2.0 * phi(r) + phi(r - h)) / (h * h);
            let q_fd = phi_p / r;
            let t_fd = (phi_pp - q_fd) / (r * r);
            assert!(
                (q - q_fd).abs() <= 1.0e-6 * q.abs().max(1.0e-300),
                "q at r={r}: stable={q:e} vs φ′/r central-diff={q_fd:e}"
            );
            assert!(
                (t - t_fd).abs() <= 1.0e-4 * t.abs().max(1.0e-300),
                "t at r={r}: stable={t:e} vs (φ″−q)/r² central-diff={t_fd:e}"
            );
        }
    }
}