gam-terms 0.3.152

use super::*;

/// Relative native-penalty ridge kept on Duchon's affine trend directions.
///
/// The affine block is also exposed as a `DoublePenaltyNullspace` component so
/// REML can shrink unused slopes.  Letting that component be driven all the way
/// to zero, however, can make the realized unpenalized design rank-deficient
/// after `OrthogonalToParametric` has projected constants/trends against the
/// constrained kernel block: the slopes then count as unpenalized columns even
/// when they are aliased by the realized design.  A machine-scale ridge in the
/// always-present native penalty keeps those slope directions structurally
/// penalized without changing the Duchon Hilbert scale at statistical scale.
///
/// This is a RELATIVE coefficient (`√ε ≈ 2⁻²⁶`): the ridge placed on each affine
/// slope column is this fraction of the curvature block's mean diagonal, NOT an
/// absolute constant. An absolute floor over-penalizes a low-magnitude curvature
/// Gram (few centers / small support): once the whole penalty is Frobenius-
/// normalized the absolute ridge no longer sits below the curvature scale, so the
/// affine slopes leave the penalty's null space and the smooth loses the affine
/// trend it is supposed to leave free at statistical scale (gam#880). Scaling by
/// the curvature magnitude keeps the ridge machine-scale relative to the penalty
/// in every configuration, so the affine trend stays in the effective null space
/// while the slopes remain structurally (non-zero) penalized.
pub(crate) const DUCHON_AFFINE_NATIVE_RIDGE_REL: f64 = 1.490_116_119_384_765_6e-8;

/// N-D periodic-cyclic-B-spline first-derivative jet `∂Φ̃/∂t` per row.
///
/// One-dimensional periodic B-spline basis (one latent axis). `t` is the
/// `(n_rows, 1)` latent matrix; each row evaluates a length-`num_basis`
/// derivative stencil w.r.t. the scalar latent coordinate. The result is
/// `(n_rows, num_basis, 1)`. This is the derivative of the row-normalized
/// design returned by [`build_periodic_bspline_basis_1d`]. The raw
/// derivative formula `B'_i(x) = (B_{i,k−1}(x) − B_{i+1,k−1}(x)) / h` is
/// evaluated alongside the unnormalized basis row `Φ`; the returned row uses
/// the quotient rule for `Φ̃ = Φ / S`, where `S = Σ_j Φ_j`.
pub fn periodic_bspline_first_derivative_nd(
    t: ArrayView2<'_, f64>,
    data_range: (f64, f64),
    degree: usize,
    num_basis: usize,
) -> Result<Array3<f64>, BasisError> {
    if t.ncols() != 1 {
        crate::bail_invalid_basis!(
            "periodic_bspline_first_derivative_nd: t must have exactly 1 column; got {}",
            t.ncols()
        );
    }
    if degree == 0 {
        crate::bail_invalid_basis!("periodic_bspline_first_derivative_nd requires degree >= 1");
    }
    if num_basis < degree + 1 {
        crate::bail_invalid_basis!(
            "periodic_bspline_first_derivative_nd requires num_basis >= degree + 1 (got num_basis={num_basis}, degree={degree})"
        );
    }
    let (start, end) = data_range;
    if !(start.is_finite() && end.is_finite()) || end <= start {
        crate::bail_invalid_basis!(
            "periodic_bspline_first_derivative_nd: data_range must be finite and ordered, got {data_range:?}"
        );
    }
    let period = end - start;
    let n_rows = t.nrows();
    let t_col = t.column(0);

    let mut phi = vec![0.0_f64; num_basis];
    let mut dphi = vec![0.0_f64; num_basis];
    let mut out = Array3::<f64>::zeros((n_rows, num_basis, 1));
    for row in 0..n_rows {
        let xi = t_col[row];
        if !xi.is_finite() {
            crate::bail_invalid_basis!(
                "periodic_bspline_first_derivative_nd: non-finite latent at row {row}"
            );
        }
        let rowsum =
            fill_periodic_bspline_unnormalized_value_row(xi, start, period, degree, &mut phi);
        if !rowsum.is_finite() || rowsum <= 0.0 {
            crate::bail_invalid_basis!(
                "periodic_bspline_first_derivative_nd: non-positive rowsum at row {row}: {rowsum}"
            );
        }
        let rowsum_derivative =
            fill_periodic_bspline_unnormalized_derivative_row(xi, start, period, degree, &mut dphi);
        if !rowsum_derivative.is_finite() {
            crate::bail_invalid_basis!(
                "periodic_bspline_first_derivative_nd: non-finite rowsum derivative at row {row}: {rowsum_derivative}"
            );
        }
        let rowsum_squared = rowsum * rowsum;
        for i in 0..num_basis {
            out[[row, i, 0]] = dphi[i] / rowsum - phi[i] * rowsum_derivative / rowsum_squared;
        }
    }
    Ok(out)
}

/// Tensor-product 1-D-B-spline first-derivative jet `∂Φ/∂t` per row.
///
/// `t` is the `(n_rows, n_axes)` latent matrix and each axis carries its
/// own `(knots, degree)` univariate B-spline. The tensor-product basis is
///
/// ```text
///     Φ_{n, k}(t_n) = ∏_a B^{(a)}_{j_a(k)}(t_{n,a}),
/// ```
///
/// where `k` enumerates the row-major tensor product
/// `j_0 ∈ [0, K_0) × … × j_{n_axes−1} ∈ [0, K_{n_axes−1})`. The product
/// rule then gives, for the partial w.r.t. axis `axis`:
///
/// ```text
///     ∂Φ_{n,k} / ∂t_{n, axis}
///         = (B^{(axis)}_{j_axis})'(t_{n, axis})
///           · ∏_{a ≠ axis} B^{(a)}_{j_a}(t_{n,a}).
/// ```
///
/// Returned tensor shape: `(n_rows, K_total, n_axes)` where
/// `K_total = ∏_a K_a` and `K_a = knots[a].len() − degree[a] − 1`.
pub fn bspline_tensor_first_derivative(
    t: ArrayView2<'_, f64>,
    knots_per_axis: &[ArrayView1<'_, f64>],
    degrees: &[usize],
) -> Result<Array3<f64>, BasisError> {
    let n_axes = t.ncols();
    if knots_per_axis.len() != n_axes || degrees.len() != n_axes {
        crate::bail_invalid_basis!(
            "bspline_tensor_first_derivative: t has {n_axes} axes but received \
             {} knot vectors and {} degrees",
            knots_per_axis.len(),
            degrees.len(),
        );
    }
    if n_axes == 0 {
        crate::bail_invalid_basis!(
            "bspline_tensor_first_derivative: t must have at least one axis".into(),
        );
    }
    let n_rows = t.nrows();
    // Per-axis basis sizes and total tensor size.
    let mut k_per_axis = Vec::<usize>::with_capacity(n_axes);
    let mut total = 1usize;
    for a in 0..n_axes {
        let k = knots_per_axis[a]
            .len()
            .checked_sub(degrees[a] + 1)
            .ok_or_else(|| {
                BasisError::InvalidInput(format!(
                    "bspline_tensor_first_derivative: axis {a} knot vector too short \
                     for degree {}",
                    degrees[a]
                ))
            })?;
        k_per_axis.push(k);
        total = total.checked_mul(k).ok_or_else(|| {
            BasisError::InvalidInput(
                "bspline_tensor_first_derivative: tensor-product basis size overflow".into(),
            )
        })?;
    }
    let mut out = Array3::<f64>::zeros((n_rows, total, n_axes));
    // Scratch per row: per-axis value vector and derivative vector.
    let mut values_per_axis: Vec<Vec<f64>> = k_per_axis.iter().map(|&k| vec![0.0; k]).collect();
    let mut derivs_per_axis: Vec<Vec<f64>> = k_per_axis.iter().map(|&k| vec![0.0; k]).collect();
    // Hoist per-axis scratch allocations outside the row loop. Previously each
    // row reallocated a fresh `BsplineScratch` for the value path and (via
    // `evaluate_bspline_derivative_scalar`) a fresh lower-basis `Vec<f64>` and
    // lower-degree `BsplineScratch` for the derivative path on every axis,
    // turning the tensor evaluator into O(n_rows · n_axes) heap traffic.
    let mut value_scratch_per_axis: Vec<internal::BsplineScratch> = degrees
        .iter()
        .map(|&d| internal::BsplineScratch::new(d))
        .collect();
    let mut lower_basis_per_axis: Vec<Vec<f64>> = knots_per_axis
        .iter()
        .zip(degrees.iter())
        .map(|(knots, &d)| vec![0.0; knots.len().saturating_sub(d)])
        .collect();
    let mut lower_scratch_per_axis: Vec<internal::BsplineScratch> = degrees
        .iter()
        .map(|&d| internal::BsplineScratch::new(d.saturating_sub(1)))
        .collect();
    let mut idx = vec![0usize; n_axes];
    let mut prefix = vec![1.0; n_axes + 1];
    let mut suffix = vec![1.0; n_axes + 1];
    for n in 0..n_rows {
        // Evaluate B^{(a)} and (B^{(a)})' at t_{n, a} for each axis.
        for a in 0..n_axes {
            internal::evaluate_splines_at_point_into(
                t[[n, a]],
                degrees[a],
                knots_per_axis[a],
                &mut values_per_axis[a],
                &mut value_scratch_per_axis[a],
            );
            evaluate_bspline_derivative_scalar_into(
                t[[n, a]],
                knots_per_axis[a],
                degrees[a],
                &mut derivs_per_axis[a],
                &mut lower_basis_per_axis[a],
                &mut lower_scratch_per_axis[a],
            )?;
        }
        // Enumerate tensor product in row-major order matching
        // `j = j_0 * (K_1 K_2 … K_{n_axes-1}) + j_1 * (K_2 … K_{n_axes-1}) + … + j_{n_axes-1}`.
        for k in 0..total {
            // Reconstruct multi-index `idx` from flat `k`.
            let mut rem = k;
            for a in (0..n_axes).rev() {
                idx[a] = rem % k_per_axis[a];
                rem /= k_per_axis[a];
            }

            prefix[0] = 1.0;
            for a in 0..n_axes {
                prefix[a + 1] = prefix[a] * values_per_axis[a][idx[a]];
            }
            suffix[n_axes] = 1.0;
            for a in (0..n_axes).rev() {
                suffix[a] = suffix[a + 1] * values_per_axis[a][idx[a]];
            }

            // For each output axis, derivative of axis-`axis` factor times
            // values of the others.
            for axis in 0..n_axes {
                let leave_one_out = prefix[axis] * suffix[axis + 1];
                out[[n, k, axis]] = derivs_per_axis[axis][idx[axis]] * leave_one_out;
            }
        }
    }
    Ok(out)
}

#[inline]
pub(crate) fn periodic_distance_1d(x: f64, c: f64, period: f64) -> f64 {
    let dx = (x - c).rem_euclid(period).abs();
    dx.min(period - dx).abs()
}

/// 2m-th Bernoulli polynomial ``B_{2m}(t)``, evaluated on ``t ∈ [0, 1]``.
///
/// Closed forms for the orders the Duchon stack actually uses:
///   * ``B₂(t)  = t² − t + 1/6``
///   * ``B₄(t)  = t⁴ − 2t³ + t² − 1/30``
///   * ``B₆(t)  = t⁶ − 3t⁵ + (5/2)t⁴ − (1/2)t² + 1/42``
///   * ``B₈(t)  = t⁸ − 4t⁷ + (14/3)t⁶ − (7/3)t⁴ + (2/3)t² − 1/30``
///
/// Defined for ``t ∈ [0, 1]`` then extended periodically (the caller has
/// already reduced ``r/period`` modulo 1).
pub(crate) fn even_bernoulli_polynomial(degree: usize, t: f64) -> Result<f64, BasisError> {
    let t2 = t * t;
    match degree {
        2 => Ok(t2 - t + 1.0 / 6.0),
        4 => Ok(t2 * t2 - 2.0 * t2 * t + t2 - 1.0 / 30.0),
        6 => {
            let t4 = t2 * t2;
            let t6 = t4 * t2;
            Ok(t6 - 3.0 * t4 * t + 2.5 * t4 - 0.5 * t2 + 1.0 / 42.0)
        }
        8 => {
            let t4 = t2 * t2;
            let t6 = t4 * t2;
            let t8 = t4 * t4;
            Ok(
                t8 - 4.0 * t6 * t + (14.0 / 3.0) * t6 - (7.0 / 3.0) * t4 + (2.0 / 3.0) * t2
                    - 1.0 / 30.0,
            )
        }
        other => Err(BasisError::InvalidInput(format!(
            "periodic Duchon Bernoulli kernel only implemented for B_{{2m}} with m ∈ {{1, 2, 3, 4}}; got degree {other}"
        ))),
    }
}

/// Periodic Green's function of the iterated 1D Laplacian ``(d²/dx²)^m`` on
/// the circle of circumference ``period``, modulo the constant nullspace.
///
/// Returns ``(-1)^(m+1) · B_{2m}(r / period)`` where ``B_{2m}`` is the
/// ``2m``-th Bernoulli polynomial extended periodically. The Fourier series
/// is
///
/// ```text
///     2 · (-1)^(m+1) · (2π)^{2m} / (2m)! · Σ_{n≥1} cos(2π n t) / n^{2m}
/// ```
///
/// so every nonzero harmonic carries weight ``∝ 1/n^{2m}`` and the kernel
/// matrix is full rank (modulo the constant direction) on **any** lattice of
/// ``K`` distinct circle points — uniform or not, even or odd ``K``. The
/// sign ``(-1)^(m+1)`` makes every Fourier coefficient positive, so the
/// kernel matrix is positive semidefinite with rank ``K − 1`` (a single
/// zero eigenvalue along the constants).
///
/// **Contrast with the polyharmonic kernel evaluated at wrapped distance**:
/// for ``m = 2`` the polyharmonic path computes ``φ(r) = c · r``, which is
/// the triangle wave on the circle. The triangle wave's Fourier series
/// carries only **odd** harmonics; sampled on a uniform K-lattice with even
/// K, the discrete DFT lands exactly on the zero (even-harmonic) modes and
/// the kernel matrix loses ``K/2 − 1`` singular values. The Bernoulli
/// kernel is the actual Green's function the operator demands and does not
/// suffer that lattice-parity degeneracy.
pub(crate) fn periodic_duchon_kernel_bernoulli(
    r: f64,
    m: usize,
    period: f64,
) -> Result<f64, BasisError> {
    if !period.is_finite() || period <= 0.0 {
        crate::bail_invalid_basis!(
            "periodic Duchon kernel requires positive finite period; got {period}"
        );
    }
    if m == 0 {
        crate::bail_invalid_basis!("periodic Duchon order m must be at least 1");
    }
    let t = (r / period).rem_euclid(1.0);
    let sign = if m % 2 == 1 { 1.0 } else { -1.0 };
    Ok(sign * even_bernoulli_polynomial(2 * m, t)?)
}

/// First and second derivatives ``(B'_{2m}(s), B''_{2m}(s))`` of the even
/// Bernoulli polynomial w.r.t. its argument ``s``, for the orders the Duchon
/// stack uses (``m ∈ {1, 2, 3, 4}``).
///
/// Obtained by differentiating the closed forms in [`even_bernoulli_polynomial`]
/// (each is a plain polynomial in ``s``), so they are the EXACT derivatives of
/// the forward kernel value — the analytic backward of the periodic Bernoulli
/// Green's-function design (gam#580).
pub(crate) fn even_bernoulli_polynomial_derivatives(
    degree: usize,
    s: f64,
) -> Result<(f64, f64), BasisError> {
    let s2 = s * s;
    match degree {
        2 => Ok((2.0 * s - 1.0, 2.0)),
        4 => {
            let d1 = 4.0 * s2 * s - 6.0 * s2 + 2.0 * s;
            let d2 = 12.0 * s2 - 12.0 * s + 2.0;
            Ok((d1, d2))
        }
        6 => {
            let s3 = s2 * s;
            let s4 = s2 * s2;
            let s5 = s4 * s;
            let d1 = 6.0 * s5 - 15.0 * s4 + 10.0 * s3 - s;
            let d2 = 30.0 * s4 - 60.0 * s3 + 30.0 * s2 - 1.0;
            Ok((d1, d2))
        }
        8 => {
            let s3 = s2 * s;
            let s4 = s2 * s2;
            let s5 = s4 * s;
            let s6 = s4 * s2;
            let s7 = s6 * s;
            let d1 = 8.0 * s7 - 28.0 * s6 + 28.0 * s5 - (28.0 / 3.0) * s3 + (4.0 / 3.0) * s;
            let d2 = 56.0 * s6 - 168.0 * s5 + 140.0 * s4 - 28.0 * s2 + 4.0 / 3.0;
            Ok((d1, d2))
        }
        other => Err(BasisError::InvalidInput(format!(
            "periodic Duchon Bernoulli kernel derivative only implemented for B_{{2m}} with m ∈ {{1, 2, 3, 4}}; got degree {other}"
        ))),
    }
}

/// Radial jet ``(φ, dφ/dr, d²φ/dr²)`` of the periodic Bernoulli Green's-function
/// kernel ``φ(r) = (−1)^{m+1} · B_{2m}(r / period)``.
///
/// The forward design uses ``periodic_duchon_kernel_bernoulli``; this is its
/// EXACT radial derivative so the analytic backward (the position-API VJP) is
/// consistent with the Bernoulli forward, mirroring how the polyharmonic
/// triplet feeds the non-periodic derivative path. The caller already reduces
/// the signed offset to ``[−period/2, period/2]`` and passes ``r = |offset|``
/// with the sign applied separately, so ``s = r / period ∈ [0, 1/2]`` needs no
/// further modular reduction. Each ``d/dr`` brings a ``1/period`` factor by the
/// chain rule.
pub(crate) fn periodic_duchon_kernel_bernoulli_triplet(
    r: f64,
    m: usize,
    period: f64,
) -> Result<(f64, f64, f64), BasisError> {
    if !period.is_finite() || period <= 0.0 {
        crate::bail_invalid_basis!(
            "periodic Duchon kernel requires positive finite period; got {period}"
        );
    }
    if m == 0 {
        crate::bail_invalid_basis!("periodic Duchon order m must be at least 1");
    }
    let s = (r / period).rem_euclid(1.0);
    let sign = if m % 2 == 1 { 1.0 } else { -1.0 };
    let phi = sign * even_bernoulli_polynomial(2 * m, s)?;
    let (b1, b2) = even_bernoulli_polynomial_derivatives(2 * m, s)?;
    let dphi_dr = sign * b1 / period;
    let d2phi_dr2 = sign * b2 / (period * period);
    Ok((phi, dphi_dr, d2phi_dr2))
}

// ── Exact circular periodization of the hybrid Duchon–Matérn kernel ──────────
//
// The line-space hybrid Duchon kernel has spectral density
// ``f(ρ) = ρ^{−2p}(κ²+ρ²)^{−s}`` (κ = 1/length_scale). Evaluating that kernel at
// wrapped ("cut-and-wrap") distance and building a Gram from it is NOT a
// positive-definite construction on the circle: the circular Fourier
// coefficients of ``K(d_wrap)`` are not samples of the spectral density and go
// negative once the length scale is a material fraction of the period, tripping
// the penalty PSD guard with a genuinely indefinite Gram (gam#2372,
// `min_eigenvalue ≈ −1.6e-4`).
//
// The exact periodic kernel is instead the Fourier series whose coefficients ARE
// the spectral density sampled on the circle lattice ``ρ_k = 2πk/P``:
//
// ```text
//     K_per(r) = (2/P) Σ_{k≥1} f(ρ_k) cos(2πk r/P)     (mod the constant mode).
// ```
//
// Because every coefficient ``f(ρ_k) > 0`` this is positive semidefinite on ANY
// set of circle points (its center Gram is a nonnegative combination of
// rank-one ``cos``/``sin`` outer products), full rank modulo constants. Summing
// the series in closed form uses the SAME partial-fraction split the real-space
// kernel already uses (`duchon_partial_fraction_coeffs`):
//
// ```text
//     f(ρ) = Σ_m a_m ρ^{−2m} + Σ_n b_n (κ²+ρ²)^{−n}.
// ```
//
// * each ``ρ^{−2m}`` block periodizes to the Bernoulli polynomial
//   ``S_{2m}(r) = P^{2m−1}/(2m)! · (−1)^{m+1} B_{2m}(r/P)``
//   (`periodic_duchon_kernel_bernoulli`), which is κ-independent;
// * each ``(κ²+ρ²)^{−n}`` block periodizes to the "periodic Sobolev block"
//   ``T_n(r) = (2/P) Σ_{k≥1} (κ²+ρ_k²)^{−n} cos(2πk r/P)``.
//
// The whole ``T_n`` tower comes from ONE scalar
// ``ψ₁(u) = Σ_{k≥1} cos(2πk r/P)/(u+ρ_k²)`` (``u = κ²``): since
// ``(u+ρ²)^{−n} = (−1)^{n−1}/(n−1)! ∂^{n−1}/∂u^{n−1}(u+ρ²)^{−1}`` we have
// ``ψ_n = (−1)^{n−1} c_{n−1}`` where ``c_j`` are the Taylor coefficients of ψ₁ in
// ``(u−u₀)`` and ``T_n = (2/P) ψ_n``. The naive closed form
// ``ψ₁ = (P/4κ)·cosh(κ(P/2−r))/sinh(κP/2) − 1/(2κ²)`` cancels catastrophically as
// κ→0, so ψ₁ is instead evaluated from the branch-free entire-function series
//
// ```text
//     ψ₁(u) = P̂(u) / (2 Ŝ(u)),
//     Ŝ(u) = Σ_{i≥0} h^{2i} u^i/(2i+1)!            [= sinh(√u h)/(√u h)],
//     P̂(u) = Σ_{j≥0} [ b^{2(j+1)}/(2(j+1))! − h^{2(j+1)}/(2(j+1)+1)! ] u^j,
// ```
//
// with ``h = P/2`` and ``b = P/2 − r``; the leading ``u⁰`` term of the two
// P̂ contributions cancels analytically, so no float cancellation remains. Every
// series term is generated by a bounded ratio recurrence (no factorial/`√u`
// overflow), giving machine-precision kernel values across the reachable
// length-scale range (verified against the direct spectral sum for
// κ ∈ [0.5, 10³]).

/// Order-`order` truncated Taylor series in one variable, storing the
/// coefficients `c[j]` of `(x − x₀)^j`. The ops used to periodize the hybrid
/// Duchon kernel — Cauchy product and series reciprocal — are exact term
/// recurrences, so the extracted derivatives are analytic (no finite
/// differencing).
#[derive(Clone)]
struct PeriodicKernelJet {
    c: Vec<f64>,
}

impl PeriodicKernelJet {
    fn zeros(order: usize) -> Self {
        Self {
            c: vec![0.0; order + 1],
        }
    }

    /// Jet of the identity `x` seeded at `x₀` (`[x₀, 1, 0, …]`).
    fn seed(order: usize, x0: f64) -> Self {
        let mut c = vec![0.0; order + 1];
        c[0] = x0;
        if order >= 1 {
            c[1] = 1.0;
        }
        Self { c }
    }

    fn order(&self) -> usize {
        self.c.len() - 1
    }

    /// Truncated Cauchy product `self · other`.
    fn mul(&self, other: &Self) -> Self {
        let n = self.order();
        let mut out = vec![0.0; n + 1];
        for i in 0..=n {
            let a = self.c[i];
            if a == 0.0 {
                continue;
            }
            for j in 0..=(n - i) {
                out[i + j] += a * other.c[j];
            }
        }
        Self { c: out }
    }

    /// Multiply every coefficient by the scalar `k` (`self · k` as a series).
    fn scaled(&self, k: f64) -> Self {
        Self {
            c: self.c.iter().map(|&v| v * k).collect(),
        }
    }

    /// Series reciprocal `1 / self`; requires a non-zero constant term.
    fn recip(&self) -> Self {
        let n = self.order();
        let mut h = vec![0.0; n + 1];
        h[0] = 1.0 / self.c[0];
        for k in 1..=n {
            let mut s = 0.0;
            for j in 1..=k {
                s += self.c[j] * h[k - j];
            }
            h[k] = -s * h[0];
        }
        Self { c: h }
    }
}

/// Taylor coefficients (in `u − κ²`) of the circular Sobolev scalar
/// ``ψ₁(u) = Σ_{k≥1} cos(2πk r/P)/(u + ρ_k²)`` up to `order`, from which the
/// periodized Matérn block ``T_n(r) = (2/P)·(−1)^{n−1}·c[n−1]`` is read off for
/// ``n = 1 ..= order+1`` (see the module note above). `r` is the circular
/// distance (already reduced to `[0, P/2]`).
fn periodic_sobolev_block_psi_jet(r: f64, kappa: f64, period: f64, order: usize) -> Vec<f64> {
    let h = 0.5 * period;
    let b = h - r;
    let u0 = kappa * kappa;
    let uj = PeriodicKernelJet::seed(order, u0);

    // Ŝ(u) = Σ_{i≥0} h^{2i}/(2i+1)! u^i ; running term st_i, ratio h²/((2i)(2i+1)).
    // P̂(u) = Σ_{j≥0} w_{j+1} u^j with w_m = b^{2m}/(2m)! − h^{2m}/(2m+1)!,
    //   split into PB (the b term, j=0 seed b²/2!) and PH (the h term, seed h²/3!)
    //   so each is a single positive series advanced by a bounded ratio.
    let mut sh = PeriodicKernelJet::zeros(order);
    let mut pb = PeriodicKernelJet::zeros(order);
    let mut ph = PeriodicKernelJet::zeros(order);
    let mut st = PeriodicKernelJet::zeros(order);
    let mut pbt = PeriodicKernelJet::zeros(order);
    let mut pht = PeriodicKernelJet::zeros(order);
    st.c[0] = 1.0; // i = 0: h⁰/1! = 1
    pbt.c[0] = 0.5 * b * b; // j = 0: b²/2!
    pht.c[0] = h * h / 6.0; // j = 0: h²/3!
    for t in 0..=order {
        sh.c[t] += st.c[t];
        pb.c[t] += pbt.c[t];
        ph.c[t] += pht.c[t];
    }
    let h2 = h * h;
    let b2 = b * b;
    // Enough terms to converge sinh(√u₀ h)/(√u₀ h); the guard below exits early.
    const MAX_TERMS: usize = 8192;
    for i in 1..=MAX_TERMS {
        let fi = i as f64;
        let ratio_s = h2 / ((2.0 * fi) * (2.0 * fi + 1.0));
        let ratio_pb = b2 / ((2.0 * fi + 1.0) * (2.0 * fi + 2.0));
        let ratio_ph = h2 / ((2.0 * fi + 2.0) * (2.0 * fi + 3.0));
        st = st.mul(&uj).scaled(ratio_s);
        pbt = pbt.mul(&uj).scaled(ratio_pb);
        pht = pht.mul(&uj).scaled(ratio_ph);
        for t in 0..=order {
            sh.c[t] += st.c[t];
            pb.c[t] += pbt.c[t];
            ph.c[t] += pht.c[t];
        }
        if i > 4
            && st.c[0].abs() < 1e-18 * sh.c[0].abs().max(f64::MIN_POSITIVE)
            && pbt.c[0].abs() < 1e-18 * pb.c[0].abs().max(f64::MIN_POSITIVE)
            && pht.c[0].abs() < 1e-18 * ph.c[0].abs().max(f64::MIN_POSITIVE)
        {
            break;
        }
    }
    let p_hat = PeriodicKernelJet {
        c: (0..=order).map(|t| pb.c[t] - ph.c[t]).collect(),
    };
    p_hat.mul(&sh.recip()).scaled(0.5).c
}

/// Value of the exact circular (periodized) hybrid Duchon–Matérn kernel of
/// spectral orders `(p_order, s_order)` at circular distance `r`, with
/// ``κ = 1/length_scale``. This is
/// ``(2/P) Σ_{k≥1} ρ_k^{−2p}(κ²+ρ_k²)^{−s} cos(2πk r/P)`` in closed form — PSD by
/// construction (nonnegative spectral samples) and exact (no image-sum
/// truncation). `p_order` covers the Bernoulli blocks (`m ∈ {1..=4}`).
pub(crate) fn periodic_hybrid_duchon_kernel_value(
    r: f64,
    kappa: f64,
    p_order: usize,
    s_order: usize,
    period: f64,
) -> Result<f64, BasisError> {
    if !period.is_finite() || period <= 0.0 {
        crate::bail_invalid_basis!(
            "periodic hybrid Duchon kernel requires positive finite period; got {period}"
        );
    }
    if !kappa.is_finite() || kappa <= 0.0 {
        crate::bail_invalid_basis!(
            "periodic hybrid Duchon kernel requires positive finite κ; got {kappa}"
        );
    }
    let coeffs = duchon_partial_fraction_coeffs(p_order, s_order, kappa);
    let mut val = 0.0_f64;
    // κ-independent polyharmonic (Bernoulli) blocks: a_m · S_{2m}.
    for (m, &a_m) in coeffs.a.iter().enumerate().skip(1) {
        if a_m == 0.0 {
            continue;
        }
        let s_2m = periodic_bernoulli_block(r, m, period)?;
        val += a_m * s_2m;
    }
    // Periodic Matérn blocks: b_n · T_n, T_n = (2/P)(−1)^{n−1} ψ_n.
    if s_order >= 1 {
        let psi = periodic_sobolev_block_psi_jet(r, kappa, period, s_order - 1);
        for (n, &b_n) in coeffs.b.iter().enumerate().skip(1) {
            if b_n == 0.0 {
                continue;
            }
            let sign = if (n - 1) % 2 == 0 { 1.0 } else { -1.0 };
            let t_n = (2.0 / period) * sign * psi[n - 1];
            val += b_n * t_n;
        }
    }
    if !val.is_finite() {
        crate::bail_invalid_basis!("periodic hybrid Duchon kernel produced a non-finite value");
    }
    Ok(val)
}

/// The κ-independent periodized ``ρ^{−2m}`` block
/// ``S_{2m}(r) = P^{2m−1}/(2m)! · (−1)^{m+1} B_{2m}(r/P)``, expressed via
/// [`periodic_duchon_kernel_bernoulli`] (which already carries the
/// ``(−1)^{m+1} B_{2m}`` sign/value).
fn periodic_bernoulli_block(r: f64, m: usize, period: f64) -> Result<f64, BasisError> {
    let bern = periodic_duchon_kernel_bernoulli(r, m, period)?;
    let scale = period.powi((2 * m - 1) as i32) / factorial_f64(2 * m);
    Ok(scale * bern)
}

/// Value and log-κ derivatives ``(K_per, ∂K_per/∂ψ, ∂²K_per/∂ψ²)`` of the
/// periodic hybrid Duchon kernel, with ``ψ = ln κ``. The ``ρ^{−2p}`` blocks are
/// κ-independent and only the ``(κ²+ρ²)^{−s}`` factor carries ψ; differentiating
/// the spectrum gives the exact tower relations
///
/// ```text
///     ∂K^{(p,s)}/∂ψ   = −2s κ² · K^{(p,s+1)},
///     ∂²K^{(p,s)}/∂ψ² = −4s κ² · K^{(p,s+1)} + 4s(s+1) κ⁴ · K^{(p,s+2)},
/// ```
///
/// so the derivative path reuses the SAME closed-form periodization at spectral
/// orders `s`, `s+1`, `s+2` — no separate differentiation of the cosh/sinh
/// chains.
pub(crate) fn periodic_hybrid_duchon_kernel_psi_triplet(
    r: f64,
    kappa: f64,
    p_order: usize,
    s_order: usize,
    period: f64,
) -> Result<(f64, f64, f64), BasisError> {
    let k2 = kappa * kappa;
    let value = periodic_hybrid_duchon_kernel_value(r, kappa, p_order, s_order, period)?;
    if s_order == 0 {
        // Pure polyharmonic spectrum: κ-independent, so both ψ-derivatives vanish.
        return Ok((value, 0.0, 0.0));
    }
    let s = s_order as f64;
    let k_s1 = periodic_hybrid_duchon_kernel_value(r, kappa, p_order, s_order + 1, period)?;
    let k_s2 = periodic_hybrid_duchon_kernel_value(r, kappa, p_order, s_order + 2, period)?;
    let d_psi = -2.0 * s * k2 * k_s1;
    let d_psi_psi = -4.0 * s * k2 * k_s1 + 4.0 * s * (s + 1.0) * k2 * k2 * k_s2;
    Ok((value, d_psi, d_psi_psi))
}

/// Scaled Bernoulli function ``kᵥ(t) = Bᵥ(t) / ν!`` and its first derivative
/// ``k'ᵥ(t) = B'ᵥ(t)/ν! = Bᵥ₋₁(t)/(ν−1)! = kᵥ₋₁(t)`` for the degrees the
/// mixed-periodicity Sobolev kernel needs (``ν ∈ {0,1,2,3,4}``).
///
/// These are the standard Sobolev-spline reproducing-kernel building blocks
/// (Wahba 1990; Gu, *Smoothing Spline ANOVA*). ``Bᵥ`` is the (ordinary, not
/// periodised) Bernoulli polynomial:
///   ``B₀=1``, ``B₁=t−½``, ``B₂=t²−t+1/6``, ``B₃=t³−(3/2)t²+(1/2)t``,
///   ``B₄=t⁴−2t³+t²−1/30``.
fn scaled_bernoulli_value_and_derivative(nu: usize, t: f64) -> Result<(f64, f64), BasisError> {
    // (value of Bᵥ(t), value of B'ᵥ(t)=ν·Bᵥ₋₁(t)); we then divide by ν!.
    let (bv, dbv) = match nu {
        0 => (1.0, 0.0),
        1 => (t - 0.5, 1.0),
        2 => (t * t - t + 1.0 / 6.0, 2.0 * t - 1.0),
        3 => {
            let t2 = t * t;
            (t2 * t - 1.5 * t2 + 0.5 * t, 3.0 * t2 - 3.0 * t + 0.5)
        }
        4 => {
            let t2 = t * t;
            (
                t2 * t2 - 2.0 * t2 * t + t2 - 1.0 / 30.0,
                4.0 * t2 * t - 6.0 * t2 + 2.0 * t,
            )
        }
        other => {
            crate::bail_invalid_basis!(
                "mixed-periodicity Sobolev kernel needs Bernoulli degree ν ≤ 4; got {other}"
            )
        }
    };
    let factorial = (1..=nu).map(|v| v as f64).product::<f64>();
    Ok((bv / factorial, dbv / factorial))
}

/// Penalised part of the 1-D Sobolev (smoothing-spline) reproducing kernel of
/// order ``m`` on ``[0, 1]`` for a NON-periodic axis:
///
/// ```text
///     R(x, y) = kₘ(x) kₘ(y) + (−1)^{m−1} k_{2m}(|x − y|)
/// ```
///
/// with ``kᵥ(t) = Bᵥ(t)/ν!`` (Wahba 1990; Gu 2013, the ``R₁`` reproducing
/// kernel of the seminorm ``∫ (f^{(m)})²``). This kernel is **positive
/// semidefinite** and its null space is exactly the polynomials of degree
/// ``< m`` — precisely the unpenalised directions the cylinder/torus Duchon
/// nullspace must contain on a non-periodic axis (gam#1423). Using it as the
/// per-axis factor (instead of the conditionally-PD chord-polyharmonic kernel)
/// is what restores positive semidefiniteness to the mixed-periodicity penalty
/// (gam#1422).
///
/// The caller passes the RAW axis coordinates ``x, y`` together with the
/// per-axis ``(lo, hi)`` from the centers; both are affine-mapped to ``[0, 1]``
/// so ``Bᵥ`` is evaluated on its canonical domain and the same map is replayed
/// identically at prediction time.
fn nonperiodic_sobolev_kernel_1d(
    x: f64,
    y: f64,
    m: usize,
    lo: f64,
    hi: f64,
) -> Result<f64, BasisError> {
    if m == 0 {
        crate::bail_invalid_basis!("non-periodic Sobolev kernel order m must be at least 1");
    }
    let span = (hi - lo).max(1e-300);
    let xs = ((x - lo) / span).clamp(0.0, 1.0);
    let ys = ((y - lo) / span).clamp(0.0, 1.0);
    let (kx, _) = scaled_bernoulli_value_and_derivative(m, xs)?;
    let (ky, _) = scaled_bernoulli_value_and_derivative(m, ys)?;
    let diff = (xs - ys).abs();
    let sign = if m % 2 == 1 { 1.0 } else { -1.0 };
    let k2m = even_bernoulli_polynomial(2 * m, diff)? / factorial_f64(2 * m);
    Ok(kx * ky + sign * k2m)
}

/// Radial-style jet ``(R, ∂R/∂x, ∂²R/∂x²)`` of
/// [`nonperiodic_sobolev_kernel_1d`] w.r.t. the first (data) coordinate ``x``,
/// for the prediction/position-API path. Derivatives carry the chain-rule
/// ``1/span`` factor from the affine map to ``[0, 1]``; the ``|x − y|`` term's
/// first derivative picks up ``sign(x − y)`` (its second derivative is the even
/// Bernoulli second derivative, continuous across ``x = y``).
pub(crate) fn nonperiodic_sobolev_kernel_1d_triplet(
    x: f64,
    y: f64,
    m: usize,
    lo: f64,
    hi: f64,
) -> Result<(f64, f64, f64), BasisError> {
    if m == 0 {
        crate::bail_invalid_basis!("non-periodic Sobolev kernel order m must be at least 1");
    }
    let span = (hi - lo).max(1e-300);
    let xs = ((x - lo) / span).clamp(0.0, 1.0);
    let ys = ((y - lo) / span).clamp(0.0, 1.0);
    let (kx, dkx) = scaled_bernoulli_value_and_derivative(m, xs)?;
    let (ky, _) = scaled_bernoulli_value_and_derivative(m, ys)?;
    // d/dx [kₘ(xs)] = k'ₘ(xs)/span; d²/dx² = k''ₘ(xs)/span². k''ₘ = kₘ₋₂ etc.;
    // reuse the even-Bernoulli second-derivative only via the |x−y| term and
    // build kₘ's second derivative from the (m−1) scaled-Bernoulli derivative.
    let (_, d2kx_inner) = if m >= 1 {
        scaled_bernoulli_value_and_derivative(m.saturating_sub(1), xs)?
    } else {
        (0.0, 0.0)
    };
    let sign = if m % 2 == 1 { 1.0 } else { -1.0 };
    let diff = xs - ys;
    let adiff = diff.abs();
    let (b1, b2) = even_bernoulli_polynomial_derivatives(2 * m, adiff)?;
    let fac2m = factorial_f64(2 * m);
    let k2m = even_bernoulli_polynomial(2 * m, adiff)? / fac2m;
    let dsign = if diff >= 0.0 { 1.0 } else { -1.0 };

    let value = kx * ky + sign * k2m;
    // ∂/∂x: kₘ'(xs)·ky/span + sign·(B'_{2m}(|d|)/((2m)!))·sgn(d)/span
    let d1 = (dkx * ky / span) + sign * (b1 / fac2m) * dsign / span;
    // ∂²/∂x²: kₘ''(xs)·ky/span² + sign·B''_{2m}(|d|)/((2m)!)/span²
    // kₘ''(xs) = (d/dxs) k'ₘ(xs) = (d/dxs) kₘ₋₁(xs) = k'ₘ₋₁(xs) = d2kx_inner.
    let d2 = (d2kx_inner * ky / (span * span)) + sign * (b2 / fac2m) / (span * span);
    Ok((value, d1, d2))
}

#[inline]
fn factorial_f64(n: usize) -> f64 {
    (1..=n).map(|v| v as f64).product::<f64>()
}

/// Per-axis ``[lo, hi]`` bounds of the centers along every NON-periodic axis,
/// used to affine-map that axis to ``[0, 1]`` for the Sobolev kernel. Periodic
/// axes carry a placeholder (their kernel uses the period, not these bounds).
/// Mirrored byte-for-byte between the forward builder and the prediction/jet
/// path so the realized design and the frozen-center replay agree.
pub(crate) fn mixed_periodicity_axis_bounds(
    centers: ArrayView2<'_, f64>,
    periodic_per_axis: &[bool],
) -> Vec<(f64, f64)> {
    let d = centers.ncols();
    (0..d)
        .map(|j| {
            if periodic_per_axis[j] {
                (0.0, 1.0)
            } else {
                let col = centers.column(j);
                let lo = col.iter().fold(f64::INFINITY, |a, &b| a.min(b));
                let hi = col.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));
                (lo, hi)
            }
        })
        .collect()
}

/// Additive mixed-periodicity Duchon kernel value
/// ``K(x, c) = Σ_j R_j(x_j, c_j)``, the sum of per-axis positive-semidefinite
/// reproducing kernels: the periodic Bernoulli Green's function on periodic
/// axes and the 1-D Sobolev kernel on non-periodic axes. As a SUM of PSD
/// kernels it is PSD, and its null space is the polynomials of degree ``< m``
/// in the non-periodic coordinates only (constants on the periodic axes),
/// which is exactly the correct cylinder/torus Duchon null space
/// (gam#1422 / gam#1423). This replaces the conditionally-PD
/// polyharmonic-of-chord-distance kernel.
pub(crate) fn mixed_periodicity_additive_kernel(
    x: ArrayView1<'_, f64>,
    c: ArrayView1<'_, f64>,
    m: usize,
    periodic_per_axis: &[bool],
    periods: &[f64],
    axis_bounds: &[(f64, f64)],
) -> Result<f64, BasisError> {
    let d = x.len();
    let mut acc = 0.0_f64;
    for j in 0..d {
        acc += if periodic_per_axis[j] {
            let r = periodic_distance_1d(x[j], c[j], periods[j]);
            periodic_duchon_kernel_bernoulli(r, m, periods[j])?
        } else {
            let (lo, hi) = axis_bounds[j];
            nonperiodic_sobolev_kernel_1d(x[j], c[j], m, lo, hi)?
        };
    }
    Ok(acc)
}

/// Per-axis value + first/second self-derivative of the additive
/// mixed-periodicity kernel for one ``(x, c)`` pair. Because the kernel is the
/// SUM of per-axis 1-D kernels, the gradient w.r.t. ``x`` is per-axis
/// (``∂K/∂x_a = R_a'(x_a, c_a)``, no cross terms) and the Hessian is DIAGONAL
/// (``∂²K/∂x_a∂x_c = δ_{ac} R_a''``). Returns ``(value, grad_a, hess_aa)`` with
/// length-``d`` per-axis vectors, so the caller can assemble the input-location
/// jet/Hessian directly. This is the exact analytic derivative of
/// [`mixed_periodicity_additive_kernel`].
pub(crate) fn mixed_periodicity_additive_kernel_jet(
    x: ArrayView1<'_, f64>,
    c: ArrayView1<'_, f64>,
    m: usize,
    periodic_per_axis: &[bool],
    periods: &[f64],
    axis_bounds: &[(f64, f64)],
) -> Result<(f64, Vec<f64>, Vec<f64>), BasisError> {
    let d = x.len();
    let mut value = 0.0_f64;
    let mut grad = vec![0.0_f64; d];
    let mut hess = vec![0.0_f64; d];
    for a in 0..d {
        let (v, d1, d2) = if periodic_per_axis[a] {
            // The Bernoulli triplet differentiates w.r.t. the unsigned radial
            // distance `r = |x_a − c_a|` reduced mod period; convert to the
            // derivative w.r.t. `x_a` via the chain rule (sign of the reduced
            // signed offset). The second derivative is sign-independent.
            let p = periods[a];
            let signed = {
                let raw = (x[a] - c[a]).rem_euclid(p);
                if raw > 0.5 * p { raw - p } else { raw }
            };
            let r = signed.abs();
            let (phi, dphi_dr, d2phi_dr2) = periodic_duchon_kernel_bernoulli_triplet(r, m, p)?;
            let dsign = if signed >= 0.0 { 1.0 } else { -1.0 };
            (phi, dphi_dr * dsign, d2phi_dr2)
        } else {
            let (lo, hi) = axis_bounds[a];
            nonperiodic_sobolev_kernel_1d_triplet(x[a], c[a], m, lo, hi)?
        };
        value += v;
        grad[a] = d1;
        hess[a] = d2;
    }
    Ok((value, grad, hess))
}

/// Polynomial side-condition block for the mixed-periodicity Duchon null space:
/// monomials of total degree ``< m`` in the NON-periodic coordinates only
/// (periodic axes contribute only the constant, which is the degree-0 monomial
/// shared by all axes). For the cylinder ``(θ periodic, y free)`` with ``m = 2``
/// this yields ``{1, y}`` — so ``f(θ, y) = a + b y`` is correctly unpenalised
/// (gam#1423).
pub(crate) fn mixed_periodicity_nullspace_poly_block(
    points: ArrayView2<'_, f64>,
    m: usize,
    periodic_per_axis: &[bool],
) -> Array2<f64> {
    let n = points.nrows();
    let nonperiodic_axes: Vec<usize> = (0..points.ncols())
        .filter(|&j| !periodic_per_axis[j])
        .collect();
    let max_degree = m.saturating_sub(1);
    // Monomial exponents of total degree ≤ (m−1) over the non-periodic axes.
    let exps = monomial_exponents(nonperiodic_axes.len(), max_degree);
    let mut block = Array2::<f64>::zeros((n, exps.len()));
    for (col, exp) in exps.iter().enumerate() {
        for row in 0..n {
            let mut value = 1.0_f64;
            for (local_axis, &power) in exp.iter().enumerate() {
                let axis = nonperiodic_axes[local_axis];
                value *= points[[row, axis]].powi(power as i32);
            }
            block[[row, col]] = value;
        }
    }
    block
}

/// Drop centers that periodically identify with the leftmost anchor.
///
/// When the user describes a closed periodic lattice by including BOTH
/// endpoints of ``[left, left+period]``, the right endpoint is the same
/// circle point as ``left`` and produces an identical kernel column. We
/// remove every such duplicate (tested under the periodic metric with a
/// tolerance scaled to ``period``); the remaining centers correspond to
/// geometrically distinct points on the circle.
pub(crate) fn collapse_periodic_endpoint(
    centers: Array2<f64>,
    left: f64,
    period: f64,
) -> Array2<f64> {
    if period <= 0.0 || !period.is_finite() {
        return centers;
    }
    // Tolerance: relative to ``period``, well below any reasonable lattice
    // spacing (mgcv's smallest practical periodic ``k`` is around 3, giving a
    // spacing of ``period/3``).
    let tol = period.max(1.0) * 1.0e-10;
    let col = centers.column(0);
    let n_rows = col.len();
    // Keep the first center that maps to the circle point of ``left`` and
    // drop every subsequent center in the same equivalence class. A
    // naive "always keep index 0, drop other left-equivalents" rule loses
    // the geometric point entirely when the user passes centers in
    // unsorted order — e.g. ``[5, 0, period]`` would collapse to ``[5]``
    // because both ``0`` and ``period`` are left-equivalents at indices
    // ``> 0``.
    let mut seen_left = false;
    let keep: Vec<usize> = (0..n_rows)
        .filter(|&i| {
            if periodic_distance_1d(col[i], left, period) <= tol {
                if seen_left {
                    return false;
                }
                seen_left = true;
            }
            true
        })
        .collect();
    if keep.len() == n_rows {
        return centers;
    }
    let mut trimmed = Array2::<f64>::zeros((keep.len(), centers.ncols()));
    for (out_row, &src_row) in keep.iter().enumerate() {
        for c in 0..centers.ncols() {
            trimmed[[out_row, c]] = centers[[src_row, c]];
        }
    }
    trimmed
}

pub(crate) fn build_periodic_duchon_basis_1d(
    data: ArrayView2<'_, f64>,
    spec: &DuchonBasisSpec,
    centers: Array2<f64>,
    workspace: &mut BasisWorkspace,
) -> Result<BasisBuildResult, BasisError> {
    if data.ncols() != 1 {
        crate::bail_invalid_basis!(
            "periodic Duchon smooths currently require exactly one covariate"
        );
    }
    // ``left + period`` is the same circle point as ``left``. If the user
    // supplied centers spanning ``[left, left+period]`` (the natural way to
    // describe a closed periodic lattice and what the position-API validator
    // requires) the rightmost point is a duplicate of the leftmost under
    // periodic identification. Two identical kernel columns make the design
    // ``rank K−1`` instead of ``K``; ``X'X`` becomes singular (cond ~10¹⁷)
    // and the REML whitening transform amplifies machine noise into a ~10⁻⁶
    // negative eigenvalue, tripping the solver's PSD check.
    //
    // ``prepare_periodic_duchon_centers_1d_with_period`` validates the center
    // matrix, computes ``(left, period)`` and drops the periodically duplicate
    // center, in one place that every periodic Duchon code path shares. When
    // ``spec.periodic`` carries an explicit per-axis period (the position-API
    // half-open lattice path — gam#580), honor it as the domain wrap; otherwise
    // derive it from the center span (the closed lattice the formula DSL emits).
    let explicit_period = spec
        .periodic
        .as_ref()
        .and_then(|axes| axes.first().copied().flatten());
    let (centers, left, period) =
        prepare_periodic_duchon_centers_1d_with_period(centers, explicit_period)?;
    // The user encodes the Duchon order ``m`` in ``spec.nullspace_order``
    // (``Zero → m=1``, ``Linear → m=2``, ``Degree(d) → m=d+1``). Periodicity
    // forces the *constraint* nullspace to ``{constants}`` (the only
    // polynomial that is itself periodic), but the *kernel* must still
    // encode full ``m``-th-order smoothness. The right kernel for that is
    // the periodic Green's function of ``(d²/dx²)^m`` — the Bernoulli
    // polynomial ``B_{2m}(r/period)`` — not the polyharmonic kernel
    // ``r^{2p+2s-d}`` evaluated at wrapped distance (which collapses to the
    // triangle wave ``r^1`` after the periodic constraint forces ``p=1`` and
    // produces zero singular values on even-K uniform lattices).
    let user_m = duchon_p_from_nullspace_order(spec.nullspace_order);
    let effective_nullspace_order = DuchonNullspaceOrder::Zero;
    let p_order = duchon_p_from_nullspace_order(effective_nullspace_order);
    let s_order = spec.power_as_usize();
    // Validate against the INTEGER `s` the hybrid kernel actually evaluates
    // (`power_as_usize` truncates a fractional `spec.power`), so the
    // well-posedness gate matches the realized kernel rather than the raw power.
    validate_duchon_kernel_orders(spec.length_scale, p_order, s_order as f64, 1)?;
    let z = kernel_constraint_nullspace(
        centers.view(),
        effective_nullspace_order,
        &mut workspace.cache,
    )?;
    let kernel_cols = z.ncols();
    let mut basis = Array2::<f64>::zeros((data.nrows(), kernel_cols + 1));
    let coeffs = spec
        .length_scale
        .map(|ls| duchon_partial_fraction_coeffs(p_order, s_order, 1.0 / ls.max(1e-300)));
    let pure_poly_coeff = if spec.length_scale.is_none() {
        Some(PolyharmonicBlockCoeff::new(
            (pure_duchon_block_order(p_order, s_order as f64)) as f64,
            1,
        ))
    } else {
        None
    };
    let kernel_amp = duchon_kernel_amplification(
        centers.view(),
        spec.length_scale,
        p_order,
        s_order,
        1,
        None,
        coeffs.as_ref(),
        pure_poly_coeff.as_ref(),
    );
    // Step 1: build the N×K raw kernel matrix in parallel (each row is
    // independent; no shared writes). Step 2: design[:, :kernel_cols] =
    // K @ z via fast_ab (BLAS), which beats a hand-rolled per-row matvec
    // loop both at small K (compiler vectorizes the inner loop) and at
    // large K (one big matmul vs. many small ones).
    let centers_col0: Vec<f64> = centers.column(0).to_vec();
    let n_data = data.nrows();
    let k_centers = centers_col0.len();
    let len_scale = spec.length_scale;
    let mut raw_kernel = Array2::<f64>::zeros((n_data, k_centers));
    let err_flag = std::sync::atomic::AtomicBool::new(false);
    // Hoist the kernel-form choice out of the inner row × center loop. The
    // pure-Duchon vs. hybrid-Matern branch is the same for every row, so a
    // single-time dispatch saves N·K conditional branches at large scale.
    let amp = kernel_amp;
    if pure_poly_coeff.is_some() {
        // Pure polyharmonic case (no Matern length-scale). Use the periodic
        // Green's function — Bernoulli ``B_{2m}(r/period)`` — directly. This
        // is the actual Green's function of ``(d²/dx²)^m`` on the circle
        // modulo constants. Every Fourier mode contributes with weight
        // ``∝ 1/n^{2m}``, so the kernel matrix is full rank (modulo the
        // constant direction) on any K-point lattice — uniform or not, even
        // or odd K. The triangle-wave kernel ``r`` that the polyharmonic
        // dispatch would emit here only has odd Fourier modes and collapses
        // on even-K uniform lattices.
        raw_kernel
            .axis_chunks_iter_mut(ndarray::Axis(0), 1024)
            .into_par_iter()
            .enumerate()
            .for_each(|(chunk_idx, mut block)| {
                let row_offset = chunk_idx * 1024;
                for (local_i, mut out_row) in block.outer_iter_mut().enumerate() {
                    let i = row_offset + local_i;
                    let x = wrap_to_period(data[[i, 0]], left, period);
                    for j in 0..k_centers {
                        let r = periodic_distance_1d(x, centers_col0[j], period);
                        match periodic_duchon_kernel_bernoulli(r, user_m, period) {
                            Ok(v) => out_row[j] = v,
                            Err(_) => {
                                err_flag.store(true, std::sync::atomic::Ordering::Relaxed);
                                return;
                            }
                        }
                    }
                }
            });
    } else {
        // Hybrid Matérn-blended case. The line-space kernel evaluated at wrapped
        // distance is NOT positive-definite on the circle (gam#2372); use the
        // exact Fourier-series periodization of the spectral density, which is
        // PSD by construction. `kappa` is well defined here because this branch
        // runs iff `length_scale` is `Some`.
        let kappa = 1.0
            / len_scale
                .expect("hybrid branch requires length_scale")
                .max(1e-300);
        raw_kernel
            .axis_chunks_iter_mut(ndarray::Axis(0), 1024)
            .into_par_iter()
            .enumerate()
            .for_each(|(chunk_idx, mut block)| {
                let row_offset = chunk_idx * 1024;
                for (local_i, mut out_row) in block.outer_iter_mut().enumerate() {
                    let i = row_offset + local_i;
                    let x = wrap_to_period(data[[i, 0]], left, period);
                    for j in 0..k_centers {
                        let r = periodic_distance_1d(x, centers_col0[j], period);
                        match periodic_hybrid_duchon_kernel_value(
                            r, kappa, p_order, s_order, period,
                        ) {
                            Ok(v) => out_row[j] = v * amp,
                            Err(_) => {
                                err_flag.store(true, std::sync::atomic::Ordering::Relaxed);
                                return;
                            }
                        }
                    }
                }
            });
    }
    if err_flag.load(std::sync::atomic::Ordering::Relaxed) {
        crate::bail_invalid_basis!("periodic Duchon kernel evaluation produced a non-finite value");
    }
    // design[:, :kernel_cols] = raw_kernel @ z; design[:, kernel_cols] = 1
    let design_kernel = fast_ab(&raw_kernel, &z);
    basis
        .slice_mut(s![.., 0..kernel_cols])
        .assign(&design_kernel);
    basis.column_mut(kernel_cols).fill(1.0);
    let mut center_kernel = Array2::<f64>::zeros((centers.nrows(), centers.nrows()));
    fill_symmetric_from_row_kernel(&mut center_kernel, |i, j| {
        let r = periodic_distance_1d(centers[[i, 0]], centers[[j, 0]], period);
        if pure_poly_coeff.is_some() {
            // Same Bernoulli Green's function the design uses — keeps the
            // penalty ``ω = z' K_centers z`` exactly the Gram matrix of the
            // smoother in its native basis, with no scale mismatch.
            periodic_duchon_kernel_bernoulli(r, user_m, period)
        } else {
            // Same exact circular periodization the design uses, so
            // ``ω = z' K_centers z`` is the PSD Gram of the periodic smoother.
            let kappa = 1.0
                / spec
                    .length_scale
                    .expect("hybrid branch requires length_scale")
                    .max(1e-300);
            Ok(
                periodic_hybrid_duchon_kernel_value(r, kappa, p_order, s_order, period)?
                    * kernel_amp,
            )
        }
    })?;
    let omega = fast_ab(&fast_atb(&z, &center_kernel), &z);
    let mut penalty = Array2::<f64>::zeros((basis.ncols(), basis.ncols()));
    penalty
        .slice_mut(s![0..kernel_cols, 0..kernel_cols])
        .assign(&omega);
    let raw_primary =
        ConstructiveQuadratic::try_from_dense_psd(penalty, "periodic Duchon raw primary penalty")?;
    let base_design = DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(basis));
    let identifiability_transform = spatial_identifiability_transform_from_design_matrix(
        data,
        &base_design,
        &spec.identifiability,
        "periodic Duchon",
    )?;
    let (design, primary) = if let Some(transform) = identifiability_transform.as_ref() {
        let design = wrap_dense_design_with_transform(base_design, transform, "periodic Duchon")?;
        let gauge = gam_problem::Gauge::from_block_transforms(&[transform.clone()]);
        let transformed =
            raw_primary.restricted(&gauge, "periodic Duchon identified primary penalty")?;
        (design, transformed)
    } else {
        (base_design, raw_primary)
    };
    let candidates = vec![normalize_constructive_penalty_candidate(
        primary,
        PenaltySource::Primary,
    )?];
    let filtered = filter_penalty_candidates(candidates)?;
    Ok(BasisBuildResult {
        design,
        affine_offset: None,
        active_penalties: filtered.active,
        dropped_penalties: filtered.dropped,
        joint_null_rotation: None,
        metadata: BasisMetadata::Duchon {
            centers,
            // `input_scale: ONE` below; see the Duchon builder in
            // `duchon_thinplate.rs` for why that makes this original units.
            length_scale: spec.length_scale.map(crate::OriginalUnits::new),
            periodic: Some(vec![Some(period)]),
            power: spec.power,
            nullspace_order: effective_nullspace_order,
            identifiability_transform,
            input_scale: crate::IsotropicScale::ONE,
            aniso_log_scales: None,
            operator_collocation_points: None,
            radial_reparam: None,
            spectral_basis: None,
        },
        kronecker_factored: None,
    })
}

/// Build a multi-dimensional Duchon basis with per-axis periodicity
/// (cylinder ``(True, False)``, torus ``(True, True)``, …).
///
/// ## Construction (gam#1422 / gam#1423)
///
/// The penalty is built from an **additive tensor (ANOVA) reproducing
/// kernel** — the sum of per-axis positive-semidefinite 1-D reproducing
/// kernels — NOT the polyharmonic kernel evaluated at the cylinder/torus
/// chord distance. The chord-polyharmonic kernel is only *conditionally*
/// positive-definite (PD orthogonal to the chord-embedding linear span), so
/// projecting out only the constants leaves indefinite linear modes and the
/// center Gram ``Ω = Zᵀ K Z`` carries large negative eigenvalues (gam#1422).
///
///   * **periodic** axis (period ``P_j``): the periodic Bernoulli Green's
///     function ``(−1)^{m+1} B_{2m}(Δ/P_j)`` ([`periodic_duchon_kernel_bernoulli`]),
///     which is PSD (every Fourier coefficient ``∝ 1/n^{2m} > 0``) with null
///     space ``{constants}``.
///   * **non-periodic** axis: the 1-D Sobolev smoothing-spline reproducing
///     kernel ([`nonperiodic_sobolev_kernel_1d`]), which is PSD with null
///     space the polynomials of degree ``< m``.
///
/// The total center kernel ``K_CC = Σ_j R_j`` is PSD (sum of PSD), and its
/// null space is the polynomials of total degree ``< m`` in the **non-periodic
/// coordinates only** (periodic coordinates contribute only constants) — the
/// correct cylinder/torus Duchon null space (gam#1423). We build
/// ``Z = null(Pᵀ)`` from that polynomial block ([`mixed_periodicity_nullspace_poly_block`]),
/// append the matching unpenalised polynomial columns to the design, and form
/// the single Primary penalty ``Ω = Zᵀ K_CC Z``, which is PSD by congruence.
/// The per-axis kernels are evaluated on each axis's own coordinate, so the
/// design wraps cleanly at every periodic seam.
pub(crate) fn build_duchon_basis_mixed_periodicity(
    data: ArrayView2<'_, f64>,
    spec: &DuchonBasisSpec,
    centers: Array2<f64>,
    periodic_per_axis: &[bool],
    periods: &[f64],
) -> Result<BasisBuildResult, BasisError> {
    let d = data.ncols();
    if d == 0 {
        crate::bail_invalid_basis!("Duchon basis requires at least one covariate dimension");
    }
    if periodic_per_axis.len() != d {
        crate::bail_invalid_basis!(
            "periodic_per_axis must have length d={d}, got {}",
            periodic_per_axis.len()
        );
    }
    if periods.len() != d {
        crate::bail_invalid_basis!("periods must have length d={d}, got {}", periods.len());
    }
    for (j, (&per, &period)) in periodic_per_axis.iter().zip(periods.iter()).enumerate() {
        if per && !(period.is_finite() && period > 0.0) {
            crate::bail_invalid_basis!(
                "axis {j} is periodic but period={period} is not finite & positive"
            );
        }
    }
    if centers.ncols() != d {
        crate::bail_invalid_basis!(
            "centers ncols={} does not match data ncols={d}",
            centers.ncols()
        );
    }

    // Hybrid Matérn (length_scale = Some) is not supported on the
    // cylinder/torus path yet; the generalized chord distance plus the
    // partial-fraction Matérn chain has not been validated for periodic
    // axes. Surface a clear error instead of silently producing nonsense.
    if spec.length_scale.is_some() {
        crate::bail_invalid_basis!(
            "mixed-periodicity Duchon basis currently only supports the pure polyharmonic spectrum (length_scale=None)"
        );
    }
    // s_order > 0 (the Sobolev tail) is similarly unvalidated for periodic
    // axes — gate to s = 0 (pure polyharmonic).
    if spec.power != 0.0 {
        crate::bail_invalid_basis!(
            "mixed-periodicity Duchon basis currently requires power = 0 (pure polyharmonic); got power={}",
            spec.power
        );
    }

    let user_m = duchon_p_from_nullspace_order(spec.nullspace_order);
    let s_order_int = 0usize;
    validate_duchon_kernel_orders(None, user_m, s_order_int as f64, d)?;

    // gam#1422 / gam#1423 — PSD mixed-periodicity Duchon via an ADDITIVE
    // tensor (ANOVA) reproducing kernel, NOT the conditionally-PD
    // polyharmonic-of-chord-distance kernel. The center Gram is the sum of
    // per-axis positive-semidefinite reproducing kernels: the periodic
    // Bernoulli Green's function on periodic axes (PSD, null = constants) and
    // the 1-D Sobolev kernel on non-periodic axes (PSD, null = polynomials of
    // degree < m). The sum is PSD (sum of PSD), and its null space is the
    // polynomials of degree < m in the NON-periodic coordinates only — exactly
    // the cylinder/torus Duchon null space. We build `Z` from that null space
    // (`{1, y, …}` on the cylinder), append the matching polynomial columns to
    // the design (so `a + b·y` is representable AND unpenalised), and form
    // `Ω = Zᵀ K_CC Z`, which is PSD by congruence.
    let axis_bounds = mixed_periodicity_axis_bounds(centers.view(), periodic_per_axis);

    let centers_owned = centers.clone();
    let k_centers = centers_owned.nrows();
    let n_data = data.nrows();

    // Non-periodic-only polynomial side condition → translation-aware null
    // space `Z = null(Pᵀ)`. For the cylinder (m=2) `P = [1, y]`.
    let poly_block_centers =
        mixed_periodicity_nullspace_poly_block(centers_owned.view(), user_m, periodic_per_axis);
    let z = kernel_constraint_nullspace_from_matrix(poly_block_centers.view())?;
    let kernel_cols = z.ncols();
    let n_poly = poly_block_centers.ncols();

    // Row-parallel additive kernel: K[i, j] = Σ_a R_a(x_i[a], c_j[a]).
    let mut raw_kernel = Array2::<f64>::zeros((n_data, k_centers));
    let kernel_err: std::sync::Mutex<Option<BasisError>> = std::sync::Mutex::new(None);
    raw_kernel
        .axis_chunks_iter_mut(ndarray::Axis(0), 1024)
        .into_par_iter()
        .enumerate()
        .for_each(|(chunk_idx, mut block)| {
            let row_offset = chunk_idx * 1024;
            for (local_i, mut out_row) in block.outer_iter_mut().enumerate() {
                let i = row_offset + local_i;
                let x_row = data.row(i);
                for j in 0..k_centers {
                    match mixed_periodicity_additive_kernel(
                        x_row,
                        centers_owned.row(j),
                        user_m,
                        periodic_per_axis,
                        periods,
                        &axis_bounds,
                    ) {
                        Ok(v) => out_row[j] = v,
                        Err(e) => {
                            *kernel_err.lock().expect(
                                "kernel-error slot is poisoned: a worker panicked \
                                         while recording a kernel failure",
                            ) = Some(e);
                            return;
                        }
                    }
                }
            }
        });
    if let Some(e) = kernel_err.into_inner().expect(
        "kernel-error slot is poisoned: a worker panicked while recording a \
                 kernel failure",
    ) {
        return Err(e);
    }

    // Design = [K @ Z, P(data)] — the kernel columns plus the explicit
    // unpenalised polynomial columns (in the non-periodic coordinates only).
    let design_kernel = fast_ab(&raw_kernel, &z);
    let poly_block_data = mixed_periodicity_nullspace_poly_block(data, user_m, periodic_per_axis);
    let mut basis = Array2::<f64>::zeros((n_data, kernel_cols + n_poly));
    basis
        .slice_mut(s![.., 0..kernel_cols])
        .assign(&design_kernel);
    basis
        .slice_mut(s![.., kernel_cols..kernel_cols + n_poly])
        .assign(&poly_block_data);

    // Penalty: Ω = Zᵀ K_CC Z (kernel-Gram identity in the projected basis),
    // padded with zero rows/cols for the unpenalised polynomial columns. PSD
    // because K_CC is PSD and Z is real (congruence preserves PSD).
    let mut center_kernel = Array2::<f64>::zeros((k_centers, k_centers));
    fill_symmetric_from_row_kernel(&mut center_kernel, |i, j| {
        mixed_periodicity_additive_kernel(
            centers_owned.row(i),
            centers_owned.row(j),
            user_m,
            periodic_per_axis,
            periods,
            &axis_bounds,
        )
    })?;
    let omega = fast_ab(&fast_atb(&z, &center_kernel), &z);
    let mut penalty = Array2::<f64>::zeros((basis.ncols(), basis.ncols()));
    penalty
        .slice_mut(s![0..kernel_cols, 0..kernel_cols])
        .assign(&omega);
    let raw_primary = ConstructiveQuadratic::try_from_dense_psd(
        penalty,
        "mixed-periodicity Duchon raw primary penalty",
    )?;

    let base_design = DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(basis));
    let identifiability_transform = spatial_identifiability_transform_from_design_matrix(
        data,
        &base_design,
        &spec.identifiability,
        "mixed-periodicity Duchon",
    )?;
    let (design, primary) = if let Some(transform) = identifiability_transform.as_ref() {
        let design =
            wrap_dense_design_with_transform(base_design, transform, "mixed-periodicity Duchon")?;
        let gauge = gam_problem::Gauge::from_block_transforms(&[transform.clone()]);
        let transformed = raw_primary.restricted(
            &gauge,
            "mixed-periodicity Duchon identified primary penalty",
        )?;
        (design, transformed)
    } else {
        (base_design, raw_primary)
    };
    let candidates = vec![normalize_constructive_penalty_candidate(
        primary,
        PenaltySource::Primary,
    )?];
    let filtered = filter_penalty_candidates(candidates)?;
    Ok(BasisBuildResult {
        design,
        affine_offset: None,
        active_penalties: filtered.active,
        dropped_penalties: filtered.dropped,
        joint_null_rotation: None,
        metadata: BasisMetadata::Duchon {
            centers: centers_owned,
            length_scale: None,
            // `periods[j]` is always present; the metadata convention is
            // `Some(period)` only for axes the caller marked periodic.
            periodic: Some(
                periodic_per_axis
                    .iter()
                    .zip(periods.iter())
                    .map(|(&is_periodic, &period)| if is_periodic { Some(period) } else { None })
                    .collect(),
            ),
            power: spec.power,
            // Record the user's requested order so the prediction/jet replay
            // rebuilds the SAME non-periodic-only polynomial null space and the
            // SAME order-`m` additive kernel (gam#1423).
            nullspace_order: spec.nullspace_order,
            identifiability_transform,
            input_scale: crate::IsotropicScale::ONE,
            aniso_log_scales: None,
            operator_collocation_points: None,
            radial_reparam: None,
            spectral_basis: None,
        },
        kronecker_factored: None,
    })
}

/// Public driver for the mixed-periodicity Duchon basis: derives per-axis
/// ``(left_j, period_j)`` from the supplied centers (mirroring how the 1D
/// periodic path infers the period from min/max), then dispatches into
/// `build_duchon_basis_mixed_periodicity`.
///
/// `periods` may be `None` (auto-derive from centers along every periodic
/// axis) or `Some(vec![...])` (length == data.ncols(); entries for
/// non-periodic axes are ignored).
pub fn build_duchon_basis_mixed_periodicity_auto(
    data: ArrayView2<'_, f64>,
    spec: &DuchonBasisSpec,
    periodic_per_axis: &[bool],
    periods: Option<&[f64]>,
) -> Result<BasisBuildResult, BasisError> {
    let mut workspace = BasisWorkspace::default();
    let centers = select_centers_by_strategy(data, &spec.center_strategy)?;
    assert_spatial_centers_below_large_scale_cap(data.ncols(), centers.view())?;
    let d = data.ncols();
    if periodic_per_axis.len() != d {
        crate::bail_invalid_basis!(
            "periodic_per_axis must have length d={d}, got {}",
            periodic_per_axis.len()
        );
    }
    let resolved_periods: Vec<f64> = match periods {
        Some(p) => {
            if p.len() != d {
                crate::bail_invalid_basis!("periods must have length d={d}, got {}", p.len());
            }
            p.to_vec()
        }
        None => {
            // Auto-derive: along each periodic axis use (max - min) over centers.
            // Non-periodic axes get a placeholder 1.0 (unused).
            let mut out = vec![1.0_f64; d];
            for j in 0..d {
                if periodic_per_axis[j] {
                    let col = centers.column(j);
                    let left = col.iter().fold(f64::INFINITY, |a, &b| a.min(b));
                    let right = col.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));
                    if !left.is_finite() || !right.is_finite() || left >= right {
                        return Err(BasisError::InvalidRange(left, right));
                    }
                    out[j] = right - left;
                }
            }
            out
        }
    };
    // The 1D periodic circle is NOT a mixed-periodicity cylinder/torus: the
    // chord-embedding polyharmonic kernel ``φ(r) = c·r^{2m−d}`` is only
    // CONDITIONALLY positive-definite on ℝ and is genuinely indefinite under
    // the chord metric on the circle (its periodised Gram carries large
    // negative eigenvalues), so it cannot serve as a PSD penalty (gam#580).
    // The actual Green's function of ``(d²/dx²)^m`` on the circle is the
    // Bernoulli kernel built by ``build_periodic_duchon_basis_1d`` — full rank
    // modulo constants, PSD by construction. Route the 1D periodic case there
    // for EVERY caller (basis design and function-norm penalty alike) so the
    // two stay consistent; reserve the chord builder for true ``d ≥ 2``
    // cylinder/torus products where it is the right object.
    if d == 1 && periodic_per_axis[0] {
        let mut periodic_spec = spec.clone();
        periodic_spec.periodic = Some(vec![Some(resolved_periods[0])]);
        return build_periodic_duchon_basis_1d(data, &periodic_spec, centers, &mut workspace);
    }
    build_duchon_basis_mixed_periodicity(data, spec, centers, periodic_per_axis, &resolved_periods)
}

/// The magic *request-layer* default `(nullspace_order, power)` for a
/// non-periodic Euclidean Duchon basis of dimension `d`: the cubic polyharmonic
/// kernel in every dimension.
///
/// Returns an affine (`Linear`, `d+1` polynomial columns) null space and the
/// fractional spectral power `s = (d − 1)/2`. With `m = p + s = 2 + (d−1)/2` the
/// pure kernel exponent `2m − d = 3`, i.e. `φ(r) = r³` for every `d` — no order
/// escalation, no even/odd-`d` log special case. The smoothing structure is the
/// analytic native reproducing-norm Gram (`PenaltySource::Primary`) plus a
/// null-space ridge; only the global mean is left free.
///
/// This is applied by the FRONT-ENDS (formula / CLI / pyffi) when the user gives
/// no explicit `power`. The basis builder itself treats `spec.power` literally,
/// so an explicit `power = 0` is honored as `s = 0` — the integer-order Duchon
/// kernel `r²·log r` (≡ the thin-plate kernel) in even `d` — rather than being
/// upgraded to the cubic default.
pub fn duchon_cubic_default(dim: usize) -> (DuchonNullspaceOrder, f64) {
    (DuchonNullspaceOrder::Linear, (dim as f64 - 1.0) / 2.0)
}

/// Build the **analytic** Duchon penalty for a non-periodic Euclidean Duchon
/// basis: the native reproducing-norm Gram `ω = α²·Zᵀ K_CC Z` (the kernel
/// evaluated at center pairs, projected through the polynomial-constraint null
/// space `Z`) plus an analytic null-space shrinkage ridge. This is the exact
/// `(m+s)`-order Duchon seminorm — pure closed form, no quadrature — the same
/// object mgcv `bs="ds"` uses, mirroring the Matérn `double_penalty` path. The
/// design scales its kernel columns by the underflow amplification `α`, so the
/// coefficient-space penalty scales by `α²`. The null-space ridge penalizes the
/// affine trend's slope (mean-free: the constant is absorbed by the model
/// intercept) so the trend is not left fully unpenalized.
/// The constrained native bending penalty `Ω_c = α² · Zᵀ K_CC Z` (m×m, in the
/// kernel-coefficient frame, pre-identifiability). This is exactly the `omega`
/// block that `duchon_native_penalty_candidates` builds; it is exposed so the
/// data-metric radial reparameterization (#1355) can solve the generalized
/// eigenproblem `Ω_c v = μ G_c v` against the realized design Gram `G_c`.
pub(crate) fn duchon_constrained_bending_penalty(
    centers: ArrayView2<'_, f64>,
    length_scale: Option<f64>,
    power: f64,
    nullspace_order: DuchonNullspaceOrder,
    aniso_log_scales: Option<&[f64]>,
    kernel_transform: &Array2<f64>,
) -> Result<Array2<f64>, BasisError> {
    let (center_kernel, kernel_amp) = duchon_center_kernel_value_matrix(
        centers,
        length_scale,
        power,
        nullspace_order,
        aniso_log_scales,
    )?;
    duchon_constrained_bending_penalty_from_kernel(&center_kernel, kernel_amp, kernel_transform)
}

/// Exact center-pair kernel values and the chart amplification applied by the
/// Duchon design. Keeping this assembly in one place guarantees that native
/// roughness and function-metric penalties see precisely the same center chart.
pub(crate) fn duchon_center_kernel_value_matrix(
    centers: ArrayView2<'_, f64>,
    length_scale: Option<f64>,
    power: f64,
    nullspace_order: DuchonNullspaceOrder,
    aniso_log_scales: Option<&[f64]>,
) -> Result<(Array2<f64>, f64), BasisError> {
    let dim = centers.ncols();
    if dim == 0 {
        crate::bail_invalid_basis!(
            "Duchon center kernel requires centers with at least one column"
        );
    }
    let k = centers.nrows();
    let p_order = duchon_p_from_nullspace_order(nullspace_order);
    let s_int = duchon_power_to_usize(power);
    let pure = length_scale.is_none();
    let pure_poly_coeff = if pure {
        Some(PolyharmonicBlockCoeff::new(
            pure_duchon_block_order(p_order, power),
            dim,
        ))
    } else {
        None
    };
    let coeffs =
        length_scale.map(|ls| duchon_partial_fraction_coeffs(p_order, s_int, 1.0 / ls.max(1e-300)));
    let kernel_amp = duchon_kernel_amplification(
        centers,
        length_scale,
        p_order,
        s_int,
        dim,
        aniso_log_scales,
        coeffs.as_ref(),
        pure_poly_coeff.as_ref(),
    );
    let axis_scales = aniso_log_scales.map(aniso_axis_scales);

    // K_CC: kernel value at every center pair (anisotropic distance when set).
    let mut center_kernel = Array2::<f64>::zeros((k, k));
    fill_symmetric_from_row_kernel(&mut center_kernel, |i, j| {
        let r = if let Some(scales) = axis_scales.as_deref() {
            aniso_distance_rows_with_scales(centers, i, centers, j, scales)
        } else {
            euclidean_distance_rows(centers, i, centers, j)
        };
        if let Some(ppc) = pure_poly_coeff.as_ref() {
            Ok(ppc.eval(r))
        } else {
            duchon_matern_kernel_general_from_distance(
                r,
                length_scale,
                p_order,
                s_int,
                dim,
                coeffs.as_ref(),
            )
        }
    })?;

    Ok((center_kernel, kernel_amp))
}

fn duchon_constrained_bending_penalty_from_kernel(
    center_kernel: &Array2<f64>,
    kernel_amp: f64,
    kernel_transform: &Array2<f64>,
) -> Result<Array2<f64>, BasisError> {
    let amp2 = kernel_amp * kernel_amp;
    let zt_k = fast_atb(kernel_transform, center_kernel);
    let omega = fast_ab(&zt_k, kernel_transform).mapv(|value| value * amp2);

    // gam#1424 — the hybrid (Duchon–Matérn) kernel's exact spectral density
    // `ρ^{-2p}(κ²+ρ²)^{-s}` is nonnegative, so the constrained bending Gram
    // `Ω_c = α²·Zᵀ K_CC Z` is positive semidefinite in exact arithmetic. The
    // historical failure was numerical: `duchon_matern_kernel_general_from_distance`
    // assembled the kernel from an ALTERNATING partial-fraction expansion
    // (polyharmonic `r^{2m−d}` minus Matérn `r^ν K_ν(κr)` blocks) whose
    // individually-enormous terms cancel to the float noise floor at high
    // dimension / spectral power (d=16, s=7: largest block ~1e3, true value
    // ~1e-13), pushing λ_min to ≈ −0.26 after normalization. That kernel
    // evaluation now routes through the cancellation-free single-integral form
    // (`duchon_hybrid_kernel_stable_integral`), so the constrained spectrum is
    // genuinely nonnegative to machine precision. Rather than silently
    // projecting (which would mask a true loss of positive-definiteness), we
    // REJECT a materially-negative spectrum and only clamp float-noise
    // negatives — per gam#1424's required PSD check before normalization.
    reject_nonpsd_then_clamp_noise(&symmetrize_penalty(&omega))
}

/// Enforce the PSD contract on a constrained Duchon bending penalty before
/// normalization (gam#1424).
///
/// The kernel's spectral density is nonnegative, so the constrained Gram must
/// be PSD in exact arithmetic. A genuinely PSD matrix only ever carries
/// negative eigenvalues at the float noise floor; those are clamped to zero. A
/// *materially* negative eigenvalue means the numerical kernel has stopped
/// representing the stated kernel — that is rejected with a clear, actionable
/// error rather than masked, because clamping a −0.26 mode silently fabricates
/// a different penalty.
fn reject_nonpsd_then_clamp_noise(matrix: &Array2<f64>) -> Result<Array2<f64>, BasisError> {
    use faer::Side;
    use gam_linalg::faer_ndarray::FaerEigh;
    let sym = symmetrize_penalty(matrix);
    let n = sym.nrows();
    if n == 0 || n != sym.ncols() {
        return Ok(sym);
    }
    let (evals, _) = FaerEigh::eigh(&sym, Side::Lower)
        .map_err(|e| BasisError::InvalidInput(format!("Duchon penalty PSD check failed: {e}")))?;
    if evals.is_empty() {
        return Ok(sym);
    }
    let max_abs_ev = evals
        .iter()
        .copied()
        .fold(0.0_f64, |acc, v| acc.max(v.abs()));
    let min_ev = evals.iter().copied().fold(f64::INFINITY, f64::min);
    // Noise-floor tolerance in eigenvalue units, so uniform scaling of the
    // penalty does not change the PSD decision. Read from the canonical
    // penalty-spectrum cutoff itself: this block is scored for rank downstream
    // against that same cutoff, so a PSD verdict taken against a private copy
    // of its formula could disagree with the rank it is later assigned.
    let tol = spectral_tolerance(&evals);
    if min_ev < -tol {
        crate::bail_invalid_basis!(
            "Duchon constrained penalty is not positive semidefinite: λ_min={min_ev:.6e} \
             (tol=−{tol:.6e}, λ_max={max_abs_ev:.6e}). The hybrid kernel's spectral density is \
             nonnegative, so a materially-negative mode indicates the kernel evaluation lost \
             positive-definiteness numerically (see gam#1424)."
        );
    }
    // λ_min is at the noise floor: clamp the harmless negative residue to zero.
    Ok(project_penalty_to_psd_cone(&sym))
}

/// The STRUCTURAL null frame of the Duchon curvature seminorm, in the chart
/// the emitted penalties live in (#2445).
///
/// The seminorm annihilates every polynomial-block direction — a theorem of
/// the RKHS construction, not a property of the shipped matrix, which
/// deliberately carries a `√ε`-relative conditioning ridge on the affine
/// slope columns (gam#880/#1816). In the raw `(kernel | poly)` frame the null
/// space is therefore exactly the `poly_cols` trailing coordinate axes. Under
/// an outer identifiability transform `T` it becomes
/// `{γ : Tγ ∈ span(poly axes)} = null(T[..kernel_cols, :])` — the null space
/// of the KERNEL-block rows of the chart, a rank decision on an orthonormal
/// matrix whose gaps are principal-angle-sized, never the Gram's conditioning.
///
/// The full polynomial block (constant AND slopes) is what must be carried:
/// intersecting the chart with the trend (slope-only) subspace is generically
/// `{0}`, while the chart itself is what removes the constant when an
/// intercept constraint is present.
pub(crate) fn duchon_structural_trend_null_frame(
    kernel_cols: usize,
    total_cols: usize,
    outer_identifiability: Option<&Array2<f64>>,
) -> Result<Array2<f64>, BasisError> {
    let poly_cols = total_cols.saturating_sub(kernel_cols);
    match outer_identifiability {
        None => {
            let mut frame = Array2::<f64>::zeros((total_cols, poly_cols));
            for column in 0..poly_cols {
                frame[[kernel_cols + column, column]] = 1.0;
            }
            Ok(frame)
        }
        Some(transform) => {
            if transform.nrows() != total_cols {
                crate::bail_dim_basis!(
                    "Duchon structural null frame: identifiability transform has {} rows \
                     but the pre-identifiability frame has {total_cols} columns",
                    transform.nrows()
                );
            }
            let kernel_rows_t = transform.slice(s![..kernel_cols, ..]).t().to_owned();
            let (frame, _rank) =
                gam_linalg::faer_ndarray::rrqr_nullspace_basis(&kernel_rows_t, 1.0)
                    .map_err(BasisError::LinalgError)?;
            Ok(frame)
        }
    }
}

pub(crate) fn duchon_native_penalty_candidates(
    centers: ArrayView2<'_, f64>,
    length_scale: Option<f64>,
    power: f64,
    nullspace_order: DuchonNullspaceOrder,
    aniso_log_scales: Option<&[f64]>,
    kernel_transform: &Array2<f64>,
    outer_identifiability: Option<&Array2<f64>>,
) -> Result<Vec<PenaltyCandidate>, BasisError> {
    duchon_native_penalty_candidates_with_curvature(
        centers,
        length_scale,
        power,
        nullspace_order,
        aniso_log_scales,
        kernel_transform,
        outer_identifiability,
        None,
    )
}

pub(crate) fn duchon_native_penalty_candidates_with_curvature(
    centers: ArrayView2<'_, f64>,
    length_scale: Option<f64>,
    power: f64,
    nullspace_order: DuchonNullspaceOrder,
    aniso_log_scales: Option<&[f64]>,
    kernel_transform: &Array2<f64>,
    outer_identifiability: Option<&Array2<f64>>,
    reduced_curvature: Option<&Array2<f64>>,
) -> Result<Vec<PenaltyCandidate>, BasisError> {
    let dim = centers.ncols();
    if dim == 0 {
        crate::bail_invalid_basis!(
            "duchon_native_penalty_candidates: centers must have at least one column"
        );
    }
    let z = kernel_transform;
    let n_kernel = z.ncols();

    // ω = α² · Zᵀ K_CC Z, embedded in the kernel block of the
    // (n_kernel + poly) pre-identifiability frame (polynomial columns carry no
    // native roughness), then mapped through the outer identifiability `T`.
    let (center_kernel, kernel_amp) = duchon_center_kernel_value_matrix(
        centers,
        length_scale,
        power,
        nullspace_order,
        aniso_log_scales,
    )?;
    let omega = match reduced_curvature {
        Some(curvature) => {
            if curvature.dim() != (n_kernel, n_kernel) {
                crate::bail_dim_basis!(
                    "Duchon reduced spectral curvature has shape {:?}; expected ({n_kernel}, {n_kernel})",
                    curvature.dim()
                );
            }
            reject_nonpsd_then_clamp_noise(&symmetrize_penalty(curvature))?
        }
        None => duchon_constrained_bending_penalty_from_kernel(&center_kernel, kernel_amp, z)?,
    };
    let center_mean: Vec<f64> = (0..dim)
        .map(|axis| centers.column(axis).sum() / centers.nrows().max(1) as f64)
        .collect();
    let mut centered = centers.to_owned();
    for axis in 0..dim {
        let mean = center_mean[axis];
        centered.column_mut(axis).mapv_inplace(|value| value - mean);
    }
    let center_poly = polynomial_block_from_order(centered.view(), nullspace_order);
    let poly_cols = center_poly.ncols();
    let n_pre = n_kernel + poly_cols;
    // Range-floor the ill-conditioned curvature spectrum so its numerical null
    // space is exactly the polynomial null space (#1815): without this, the
    // Duchon Gram's low-curvature tail sits below `analyze_penalty_block`'s
    // `nrows·1e-10·λmax` cutoff and those retained kernel modes are classed
    // UNPENALIZED, so no `λ` collapses them and the smooth cannot reach the null
    // on an irrelevant covariate. `n_pre` is the embedded penalty dimension the
    // assembled block is later scored against.
    let omega = duchon_range_floor_curvature(&omega, n_pre)?;
    let mut primary_pre = Array2::<f64>::zeros((n_pre, n_pre));
    primary_pre
        .slice_mut(s![..n_kernel, ..n_kernel])
        .assign(&omega);
    if poly_cols > 1 {
        // Machine-scale ridge on the affine SLOPE columns (the constant column
        // `n_kernel` stays free — it is the model intercept). Scale by the
        // curvature block's mean diagonal so the ridge is `√ε`-relative to the
        // penalty, not an absolute floor that would survive Frobenius
        // normalization and push the affine slopes out of the null space on a
        // low-curvature Gram (gam#880). `omega` is PSD, so its diagonal is
        // non-negative; the mean diagonal is a scale-faithful proxy for the
        // curvature magnitude and is bounded above by `‖omega‖_F`, so the
        // normalized ridge stays ≤ `√ε/√n_kernel < 1e-8` — below the statistical
        // scale while remaining strictly positive (structurally penalized).
        let curvature_scale = if n_kernel > 0 {
            let trace: f64 = (0..n_kernel).map(|i| omega[[i, i]].abs()).sum();
            trace / n_kernel as f64
        } else {
            0.0
        };
        // Fall back to the bare relative constant only for a degenerate all-zero
        // curvature block (no kernel columns / no curvature), so the slope
        // columns still carry a strictly positive ridge.
        let affine_ridge = if curvature_scale > 0.0 {
            DUCHON_AFFINE_NATIVE_RIDGE_REL * curvature_scale
        } else {
            DUCHON_AFFINE_NATIVE_RIDGE_REL
        };
        for col in (n_kernel + 1)..n_pre {
            primary_pre[[col, col]] = affine_ridge;
        }
    }
    let primary = symmetrize(&project_penalty_matrix(&primary_pre, outer_identifiability));

    let shrink = if poly_cols > 1 {
        // Evaluate the active coefficient chart on its frozen center support.
        // This compact Gram is the function metric for the represented Duchon
        // space; it is independent of training-row multiplicities and remains
        // available on every n-free κ re-key.
        let center_kernel_design = fast_ab(&center_kernel, z).mapv(|value| value * kernel_amp);
        let mut center_design = Array2::<f64>::zeros((centers.nrows(), n_pre));
        center_design
            .slice_mut(s![.., 0..n_kernel])
            .assign(&center_kernel_design);
        center_design
            .slice_mut(s![.., n_kernel..])
            .assign(&center_poly);

        // Construct the physical trend functional in the RAW coefficient
        // chart first, then restrict it through the same outer gauge as the
        // Primary. This order is essential. The collection builder receives
        // the raw Primary/ridge pair, restricts both by its global coefficient
        // gauge, and only then rebuilds the complementary ridge from the
        // constrained Primary. Inventing a fresh trend frame after the gauge
        // is a different functional and can preserve a fifth penalty that the
        // authoritative collection correctly eliminated (#2433).
        let mut trend_frame = Array2::<f64>::zeros((n_pre, poly_cols - 1));
        for column in 1..poly_cols {
            trend_frame[[n_kernel + column, column - 1]] = 1.0;
        }
        let function_gram = symmetrize_penalty(&fast_ata(&center_design));
        // Complementary metric ridge `N(NᵀGN)Nᵀ` (range = span(trend frame)),
        // NOT the leaky metric projector `GN(NᵀGN)⁻¹NᵀG` (range = span(GN)),
        // so the constant stays in the joint null space (gam#2372).
        let raw = function_space_subspace_trend_ridge(&trend_frame, &function_gram)?;
        Some(project_penalty_matrix(&raw, outer_identifiability))
    } else {
        None
    };
    let mut out = Vec::new();
    let mut primary_candidate = normalize_penalty_candidate(primary, PenaltySource::Primary)?;
    // Declare the seminorm's structural null frame on the shipped Primary
    // (#2445): the polynomial block is null by theorem, and the `√ε` affine
    // conditioning ridge deliberately present in the matrix must not be able
    // to move that decision. The frame is what the metric-consistent
    // double-penalty rebuild consumes instead of a rank test, both in the
    // frozen-chart replay below and at the term-collection chokepoint.
    let structural_frame =
        duchon_structural_trend_null_frame(n_kernel, n_pre, outer_identifiability)?;
    primary_candidate.matrix = primary_candidate.matrix.with_structural_null_frame(
        structural_frame,
        "Duchon primary structural null declaration",
    )?;
    out.push(primary_candidate);
    if let Some(shrink) = shrink {
        out.push(normalize_penalty_candidate(
            shrink,
            PenaltySource::DoublePenaltyNullspace,
        )?);
    }
    // A frozen outer-identifiability chart must reproduce the collection
    // builder's FINAL penalty topology, not merely congruence-transform the
    // raw trend ridge. The chart can remove the last structural null direction
    // of the primary penalty; in that case a carried raw ridge no longer
    // represents `null(S_primary)` and must disappear. The collection path
    // performs this same metric-consistent rebuild after applying its global
    // gauge. Repeating it here is what makes frozen single-term κ rebuilds and
    // n-free penalty re-keys share that authoritative topology (#2433).
    if outer_identifiability.is_some() {
        let primary_candidate = out
            .iter()
            .find(|candidate| matches!(candidate.source, PenaltySource::Primary))
            .ok_or_else(|| {
                BasisError::InvalidInput(
                    "Duchon trend ridge has no primary penalty in the final coefficient chart"
                        .to_string(),
                )
            })?;
        let primary_physical = primary_candidate.matrix.scaled(
            primary_candidate.normalization_scale,
            "physical frozen-chart Duchon primary",
        )?;
        let width = primary_physical.nrows();
        for candidate in &mut out {
            if !matches!(candidate.source, PenaltySource::DoublePenaltyNullspace) {
                continue;
            }
            let ridge_physical = candidate.matrix.scaled(
                candidate.normalization_scale,
                "physical frozen-chart Duchon trend ridge",
            )?;
            match rebuild_metric_consistent_ridge(&primary_physical, &ridge_physical)? {
                Some(rebuilt) => {
                    let normalized = normalize_constructive_penalty_candidate(
                        rebuilt,
                        PenaltySource::DoublePenaltyNullspace,
                    )?;
                    candidate.matrix = normalized.matrix;
                    candidate.normalization_scale = normalized.normalization_scale;
                }
                None => {
                    candidate.matrix = ConstructiveQuadratic::zero(width);
                    candidate.normalization_scale = 1.0;
                }
            }
            candidate.kronecker_factors = None;
            candidate.op = None;
        }
    }
    Ok(out)
}

/// Farthest-point collocation points per basis center for the lower-order
/// (mass / tension) operator penalties. The sample is space-filling over the
/// data SUPPORT (density-blind — sparse and dense regions weighted alike, which
/// is the regularization you want), `m = OVERSAMPLE·k` capped at `n`: dense
/// enough to resolve the `k`-bump basis, independent of `n`.
pub(crate) const DUCHON_COLLOCATION_OVERSAMPLE: usize = 3;

/// The lower two rungs of the Hilbert scale for a Duchon smooth, as FUNCTION
/// penalties collocated on a density-blind `O(k)` farthest-point sample of the
/// data support:
///   * `mass    = Σ(f−f̄)²` — centered value-design Gram (amplitude / distance
///     from the mean; kernel block only — the affine trend's slope is governed
///     by the null-space ridge, so only the global mean stays free).
///   * `tension = Σ‖∇f‖²`  — gradient-design Gram (first-order roughness).
///
/// Curvature is intentionally NOT here: it is the EXACT RKHS reproducing-norm
/// `Primary` Gram (`duchon_native_penalty_candidates`). These two orders have no
/// convergent continuous integral for the growing polyharmonic kernel, so the
/// data-support quadrature *is* their definition — and it is `O(k)`-in-`n` (the
/// sample size does not grow with the data). Each is a plain penalty (`op = None`)
/// with its own REML λ; REML drives an unhelpful one to zero. Stiffness (`D2`) is
/// absent on purpose — `Primary` is the exact, superior curvature.
/// Emit the lower-order Hilbert-scale penalties — mass `Σ(f−f̄)²` (q=0),
/// tension `Σ‖∇f‖²` (q=1), stiffness `Σ‖∇²f‖²` (q=2) — for a Duchon smooth.
///
/// Each active order routes through the shared closed-form factory, which uses
/// the EXACT continuous reproducing-norm Gram wherever the polyharmonic
/// integral converges (UV/IR + CPD adequacy — `n`-free, the high-`d` accuracy
/// and scale win) and falls back to the `D_qᵀ D_q` quadrature otherwise. That
/// quadrature is collocated on a density-blind, space-filling `O(k)`
/// farthest-point sample of the DATA SUPPORT (`select_thin_plate_knots(data,
/// 3k)`) — never the `k` sparse centers (which under-resolve a `k`-bump basis
/// and made these penalties explode), and never all `n` (which would scale with
/// the data). The collocation `D_q` is built with `max_op = max active order`,
/// so a disabled higher order never allocates its `O(d²)`-row Hessian.
///
/// The operators use the ISOTROPIC metric (`aniso = None`): the anisotropy
/// lives entirely in the curvature (`Primary`) RKHS Gram, which carries its own
/// exact `η`-derivative. Keeping these low-order stabilizers isotropic makes
/// their `η`-gradient identically zero, so the REML anisotropy optimization
/// stays consistent without per-axis operator derivatives.
pub(crate) fn duchon_operator_penalty_candidates(
    collocation_points: ArrayView2<'_, f64>,
    centers: ArrayView2<'_, f64>,
    operator_penalties: &DuchonOperatorPenaltySpec,
    length_scale: Option<f64>,
    power: f64,
    nullspace_order: DuchonNullspaceOrder,
    per_axis_relevance: bool,
    identifiability_transform: Option<&Array2<f64>>,
    radial_reparam: Option<&Array2<f64>>,
    workspace: &mut BasisWorkspace,
) -> Result<Vec<PenaltyCandidate>, BasisError> {
    let want_mass = matches!(operator_penalties.mass, OperatorPenaltySpec::Active { .. });
    let mut want_tension = matches!(
        operator_penalties.tension,
        OperatorPenaltySpec::Active { .. }
    );
    let mut want_stiffness = matches!(
        operator_penalties.stiffness,
        OperatorPenaltySpec::Active { .. }
    );
    // Collocation validity: the gradient (D1) and Hessian (D2) operator
    // quadratures are defined only when `2(p+s) > d+1` / `> d+2` respectively
    // (mass/D0 needs only kernel existence, `2(p+s) > d`, guaranteed upstream).
    // Outside that regime the operator's radial limit is undefined, so the
    // order is SKIPPED — the higher Hilbert rungs (Primary curvature, mass,
    // trend) still regularize — rather than failing the whole basis build. E.g.
    // order=0, d=3, s=1 gives `2(p+s)=4`, so tension and stiffness drop out
    // cleanly and the smooth is curvature + mass + trend.
    let effective_order = duchon_effective_nullspace_order(centers, nullspace_order);
    let p_order = duchon_p_from_nullspace_order(effective_order);
    let dim = centers.ncols();
    let two_pps = 2.0 * (p_order as f64 + power);
    want_tension = want_tension && two_pps > dim as f64 + 1.0;
    want_stiffness = want_stiffness && two_pps > dim as f64 + 2.0;
    if !want_mass && !want_tension && !want_stiffness {
        return Ok(Vec::new());
    }
    // Effective spec carrying only the collocation-valid active orders.
    let mut effective_spec = operator_penalties.clone();
    if !want_tension {
        effective_spec.tension = OperatorPenaltySpec::Disabled;
    }
    if !want_stiffness {
        effective_spec.stiffness = OperatorPenaltySpec::Disabled;
    }
    let max_op = duchon_max_active_operator_derivative_order(&effective_spec);
    let ops = build_duchon_collocation_operator_matriceswithworkspace(
        centers,
        collocation_points,
        None,
        length_scale,
        power,
        nullspace_order,
        None,
        identifiability_transform.map(|t| t.view()),
        max_op,
        radial_reparam.map(|v| v.view()),
        workspace,
    )?;
    let kernel_nullspace = ops.kernel_nullspace_transform.as_ref();
    let poly_cols = ops.polynomial_block_cols;
    // When per-axis relevance is requested (`scale_dims`) and tension is a
    // collocation-valid active order, the single isotropic gradient penalty
    // `Σ‖∇f‖²` is REPLACED by `dim` per-axis penalties `Σ(∂f/∂x_a)²`, each its
    // own REML λ_a (ARD: REML shrinks an axis's contribution toward flat only
    // when it does not earn its keep). The isotropic-order penalties
    // (mass, stiffness) still route through the shared factory; tension is
    // removed from its spec here and re-emitted per-axis below. The D1 block
    // acts on the whole function basis, including polynomial slopes.
    let split_tension = per_axis_relevance && want_tension;
    let factory_spec = if split_tension {
        let mut spec = effective_spec.clone();
        spec.tension = OperatorPenaltySpec::Disabled;
        spec
    } else {
        effective_spec
    };
    // The collocation `D_q` already carry the kernel CPD nullspace `Z`, the
    // polynomial padding, and the identifiability transform (final β-basis), so
    // the factory's quadrature fallback `fast_ata(d_q)` is β-basis. Its
    // closed-form branch rebuilds the same β-basis from `centers` via the SAME
    // `kernel_nullspace` + `poly_cols` + `outer_identifiability`, so both
    // branches agree. q=0 mass is always the centered quadrature Gram.
    let mut candidates = if let Some(length_scale) = length_scale {
        operator_penalty_candidates_closed_form(
            centers,
            &ops.d0,
            &ops.d1,
            &ops.d2,
            &factory_spec,
            p_order,
            duchon_power_to_usize(power),
            length_scale,
            None,
            kernel_nullspace,
            poly_cols,
            identifiability_transform,
        )?
    } else {
        operator_penalty_candidates_closed_form_pure(
            centers,
            &ops.d0,
            &ops.d1,
            &ops.d2,
            &factory_spec,
            p_order,
            power,
            None,
            kernel_nullspace,
            poly_cols,
            identifiability_transform,
        )?
    };
    if split_tension {
        // `D1` rows are indexed `collocation_i · dim + axis`, so axis `a` owns
        // the strided row set `a, a+dim, a+2·dim, …`. `fast_ata` of that slice
        // is the density-blind support quadrature of `∫(∂f/∂x_a)²` in the final
        // β-basis.
        for axis in 0..dim {
            let d1_axis = ops.d1.slice(s![axis..; dim, ..]).to_owned();
            candidates.push(normalize_penalty_candidate(
                symmetrize(&fast_ata(&d1_axis)),
                PenaltySource::OperatorRelevance { axis },
            )?);
        }
    }
    Ok(candidates)
}

#[cfg(test)]
mod mixed_periodicity_psd_tests {
    //! Regression tests for gam#1422 (mixed-periodicity Duchon penalty must be
    //! PSD — the additive ANOVA reproducing kernel replaces the conditionally-PD
    //! polyharmonic-of-chord-distance kernel) and gam#1423 (the cylinder null
    //! space must contain polynomials of total degree `< m` in the NON-periodic
    //! coordinates, not just constants).
    use super::*;
    use faer::Side;
    use gam_linalg::faer_ndarray::FaerEigh;
    use ndarray::{Array2, array};

    fn cylinder_spec() -> DuchonBasisSpec {
        // m = 2 ⇒ Linear null-space order; pure polyharmonic (no length scale,
        // power = 0) — the only spectrum the mixed-periodicity path supports.
        DuchonBasisSpec {
            center_strategy: CenterStrategy::UserProvided(Array2::<f64>::zeros((0, 0))),
            periodic: None,
            length_scale: None,
            power: 0.0,
            nullspace_order: DuchonNullspaceOrder::Linear,
            identifiability: SpatialIdentifiability::None,
            aniso_log_scales: None,
            operator_penalties: DuchonOperatorPenaltySpec::default(),
            boundary: OneDimensionalBoundary::Open,
            radial_reparam: None,
        }
    }

    /// Build the cylinder (`θ` periodic on `[0, 2π]`, `y` free on `[0, 1]`)
    /// Primary penalty via the public mixed-periodicity driver and return it.
    fn cylinder_primary_penalty() -> (Array2<f64>, usize) {
        // Anchor the periodic span to exactly [0, 2π] so the auto-derived period
        // is the geometric period; this mirrors the Python cylinder fixture.
        let two_pi = std::f64::consts::TAU;
        let theta = [0.0, 0.6, 1.2, 1.9, 2.5, 3.1, 3.8, 4.4, 5.0, 5.6, two_pi];
        let y = [0.5, 0.1, 0.9, 0.3, 0.7, 0.2, 0.8, 0.4, 0.6, 0.15, 0.5];
        let mut centers = Array2::<f64>::zeros((theta.len(), 2));
        for i in 0..theta.len() {
            centers[[i, 0]] = theta[i];
            centers[[i, 1]] = y[i];
        }
        let mut spec = cylinder_spec();
        spec.center_strategy = CenterStrategy::UserProvided(centers.clone());
        let periodic_per_axis = [true, false];
        let built = build_duchon_basis_mixed_periodicity_auto(
            centers.view(),
            &spec,
            &periodic_per_axis,
            None,
        )
        .expect("cylinder mixed-periodicity basis must build");
        let penalty = built
            .active_penalties
            .iter()
            .find(|penalty| matches!(penalty.info.source, PenaltySource::Primary))
            .expect("cylinder build must emit a Primary penalty");
        (penalty.matrix.clone(), centers.nrows())
    }

    #[test]
    fn native_primary_keeps_affine_trend_structurally_penalized_gam1816() {
        // gam#1816: REML may deselect the explicit DoublePenaltyNullspace slope
        // ridge.  The Primary Duchon component must therefore retain a tiny
        // native ridge on affine trend columns so those columns never become
        // realized-unpenalized aliases of the constrained kernel block.
        let centers = array![
            [0.0, 0.0],
            [1.0, 0.0],
            [0.0, 1.0],
            [1.0, 1.0],
            [0.5, 0.2],
            [0.2, 0.7],
        ];
        let order = DuchonNullspaceOrder::Linear;
        let poly_cols = polynomial_block_from_order(centers.view(), order).ncols();
        let mut workspace = BasisWorkspace::default();
        let z = kernel_constraint_nullspace(centers.view(), order, &mut workspace.cache)
            .expect("kernel constraint nullspace must build");
        let n_kernel = z.ncols();
        let candidates =
            duchon_native_penalty_candidates(centers.view(), None, 0.0, order, None, &z, None)
                .expect("native Duchon penalties must build");
        let primary = candidates
            .iter()
            .find(|candidate| matches!(candidate.source, PenaltySource::Primary))
            .expect("Primary candidate must be present");
        for col in (n_kernel + 1)..(n_kernel + poly_cols) {
            assert!(
                primary.matrix[[col, col]] > 0.0,
                "affine trend column {col} must carry the native ridge floor"
            );
        }
    }

    /// gam#2372: the Duchon trend block is the COMPLEMENTARY metric ridge
    /// `R = N(NᵀGN)Nᵀ`, not the metric projector `GN(NᵀGN)⁻¹NᵀG`. The two agree
    /// on how they weight the trend directions (both carry `NᵀGN`, the center
    /// FUNCTION metric restricted to the structural trends — this is what makes
    /// the block covariant under a center-chart reparameterization rather than a
    /// Euclidean coefficient shrinkage `NNᵀ`), but only the ridge has range
    /// exactly `span(N)`, so it annihilates the constant and the kernel block
    /// and keeps `null(Σ λ_k S_k) = span{1}`. The projector's range `span(GN)`
    /// leaks onto the constant. This test pins the ridge identity and the
    /// constant-annihilation property, and still discriminates against a
    /// Euclidean selector (for which `NᵀGN` would be the identity).
    #[test]
    fn native_trend_ridge_acts_as_center_function_metric_on_structural_trends() {
        let centers = array![
            [0.0, 0.0],
            [1.0, 0.0],
            [0.0, 1.0],
            [1.0, 1.0],
            [0.5, 0.2],
            [0.2, 0.7],
        ];
        let order = DuchonNullspaceOrder::Linear;
        let mut workspace = BasisWorkspace::default();
        let z = kernel_constraint_nullspace(centers.view(), order, &mut workspace.cache)
            .expect("kernel constraint nullspace");
        let n_kernel = z.ncols();
        let candidates =
            duchon_native_penalty_candidates(centers.view(), None, 0.0, order, None, &z, None)
                .expect("native Duchon penalties");
        let ridge = candidates
            .iter()
            .find(|candidate| matches!(candidate.source, PenaltySource::DoublePenaltyNullspace))
            .expect("affine Duchon basis must emit a trend ridge");
        let ridge = ridge.matrix.dense() * ridge.normalization_scale;

        let (center_kernel, amplification) =
            duchon_center_kernel_value_matrix(centers.view(), None, 0.0, order, None)
                .expect("center kernel");
        let center_kernel_design = fast_ab(&center_kernel, &z).mapv(|value| value * amplification);
        let center_mean: Vec<f64> = (0..centers.ncols())
            .map(|axis| centers.column(axis).sum() / centers.nrows() as f64)
            .collect();
        let mut centered = centers.clone();
        for axis in 0..centers.ncols() {
            let mean = center_mean[axis];
            centered.column_mut(axis).mapv_inplace(|value| value - mean);
        }
        let poly = polynomial_block_from_order(centered.view(), order);
        let mut center_design = Array2::<f64>::zeros((centers.nrows(), n_kernel + poly.ncols()));
        center_design
            .slice_mut(s![.., 0..n_kernel])
            .assign(&center_kernel_design);
        center_design.slice_mut(s![.., n_kernel..]).assign(&poly);
        let gram = symmetrize_penalty(&fast_ata(&center_design));
        let mut trend_frame = Array2::<f64>::zeros((center_design.ncols(), poly.ncols() - 1));
        for column in 1..poly.ncols() {
            trend_frame[[n_kernel + column, column - 1]] = 1.0;
        }

        // (1) The shipped block is exactly the complementary metric ridge
        // `R = N(NᵀGN)Nᵀ`.
        let trend_metric = trend_frame.t().dot(&gram).dot(&trend_frame);
        let reference = symmetrize_penalty(&fast_abt(
            &fast_ab(&trend_frame, &trend_metric),
            &trend_frame,
        ));
        let ridge_scale = reference
            .iter()
            .map(|value| value.abs())
            .fold(1.0_f64, f64::max);
        let ridge_err = (&ridge - &reference)
            .iter()
            .map(|value| value.abs())
            .fold(0.0_f64, f64::max);
        assert!(
            ridge_err <= 2.0e-11 * ridge_scale,
            "Duchon trend ridge must equal N(NᵀGN)Nᵀ; error={ridge_err:.3e}, scale={ridge_scale:.3e}"
        );

        // (2) The ridge annihilates the constant (n_kernel is the constant poly
        // column) and the entire kernel block — range = span(trend frame), so
        // `null(Σ λ_k S_k) = span{1}` holds (gam#2372).
        for probe_col in 0..=n_kernel {
            let mut v = Array1::<f64>::zeros(center_design.ncols());
            v[probe_col] = 1.0;
            let rv_norm = ridge.dot(&v).iter().map(|x| x * x).sum::<f64>().sqrt();
            assert!(
                rv_norm <= 2.0e-11 * ridge_scale,
                "trend ridge must annihilate column {probe_col} (constant/kernel); ||Rv||={rv_norm:.3e}"
            );
        }

        // (3) The metric is the genuine center FUNCTION metric on trends, not a
        // Euclidean coefficient selector: `NᵀGN` is materially off-identity.
        let identity = Array2::<f64>::eye(trend_metric.nrows());
        let metric_gap = (&trend_metric - &identity)
            .iter()
            .map(|value| value.abs())
            .fold(0.0_f64, f64::max);
        assert!(
            metric_gap > 1e-3,
            "trend metric NᵀGN must be non-Euclidean (guards against coefficient shrinkage); gap={metric_gap:.3e}"
        );
    }

    #[test]
    fn cylinder_penalty_is_symmetric_psd_gam1422() {
        // gam#1422: with the conditionally-PD chord-polyharmonic kernel this
        // fixture produced λ_min ≈ −0.426 (3 materially negative eigenvalues).
        // The additive ANOVA kernel (sum of PSD per-axis reproducing kernels)
        // is PSD by construction. Tolerance matches the Python
        // `_assert_symmetric_psd` slack (1e-8); do NOT weaken it.
        let (penalty, k) = cylinder_primary_penalty();
        assert_eq!(penalty.nrows(), k);
        assert_eq!(penalty.ncols(), k);
        // Symmetry.
        for i in 0..k {
            for j in 0..k {
                assert!(
                    (penalty[[i, j]] - penalty[[j, i]]).abs() <= 1e-9,
                    "cylinder penalty must be symmetric at ({i}, {j})"
                );
            }
        }
        let sym = symmetrize(&penalty);
        let (evals, _) = FaerEigh::eigh(&sym, Side::Lower).expect("eigh");
        let lambda_min = evals.iter().copied().fold(f64::INFINITY, f64::min);
        assert!(
            lambda_min > -1e-8,
            "cylinder Duchon penalty not PSD; λ_min = {lambda_min:.3e} (gam#1422)"
        );
    }

    #[test]
    fn torus_penalty_is_psd_gam1422() {
        // gam#1422: the torus fixture previously gave λ_min ≈ −0.885 (4 negative
        // eigenvalues). With both axes periodic the additive kernel is the sum of
        // two PSD Bernoulli-Green kernels, hence PSD.
        let two_pi = std::f64::consts::TAU;
        let theta = [0.0, 0.7, 1.5, 2.3, 3.0, 3.9, 4.6, 5.4, two_pi];
        let phi = [0.0, 1.1, 2.0, 0.4, 3.3, 4.8, 5.9, 2.7, two_pi];
        let mut centers = Array2::<f64>::zeros((theta.len(), 2));
        for i in 0..theta.len() {
            centers[[i, 0]] = theta[i];
            centers[[i, 1]] = phi[i];
        }
        let mut spec = cylinder_spec();
        spec.center_strategy = CenterStrategy::UserProvided(centers.clone());
        let periodic_per_axis = [true, true];
        let built = build_duchon_basis_mixed_periodicity_auto(
            centers.view(),
            &spec,
            &periodic_per_axis,
            None,
        )
        .expect("torus mixed-periodicity basis must build");
        let penalty = built
            .active_penalties
            .iter()
            .find(|penalty| matches!(penalty.info.source, PenaltySource::Primary))
            .expect("torus build must emit a Primary penalty");
        let sym = symmetrize(&penalty.matrix);
        let (evals, _) = FaerEigh::eigh(&sym, Side::Lower).expect("eigh");
        let lambda_min = evals.iter().copied().fold(f64::INFINITY, f64::min);
        assert!(
            lambda_min > -1e-8,
            "torus Duchon penalty not PSD; λ_min = {lambda_min:.3e} (gam#1422)"
        );
    }

    #[test]
    fn cylinder_nullspace_contains_one_and_y_gam1423() {
        // gam#1423: on S¹×ℝ with m = 2 the Duchon null space is {1, y} — both the
        // constant AND the linear-in-the-nonperiodic-coordinate term have zero
        // seminorm, so both must be unpenalised. The polynomial block driving the
        // null space must therefore have exactly 2 columns: a constant column and
        // a `y` column (NOT a θ column — periodic axes contribute only the
        // constant).
        let points = array![[0.0_f64, 0.10], [1.0, 0.40], [2.0, 0.70], [3.0, 0.95],];
        let periodic_per_axis = [true, false];
        let block = mixed_periodicity_nullspace_poly_block(points.view(), 2, &periodic_per_axis);
        assert_eq!(
            block.ncols(),
            2,
            "cylinder (m=2) null space must be {{1, y}} — 2 columns (gam#1423)"
        );
        // One column must be the all-ones constant; another must equal the
        // non-periodic coordinate `y` (column 1 of `points`). The block does not
        // depend on θ (column 0).
        let n = points.nrows();
        let mut has_const = false;
        let mut has_y = false;
        let mut depends_on_theta = false;
        for col in 0..block.ncols() {
            let is_const = (0..n).all(|r| (block[[r, col]] - 1.0).abs() < 1e-12);
            let is_y = (0..n).all(|r| (block[[r, col]] - points[[r, 1]]).abs() < 1e-12);
            let is_theta = (0..n).all(|r| (block[[r, col]] - points[[r, 0]]).abs() < 1e-12);
            has_const |= is_const;
            has_y |= is_y;
            depends_on_theta |= is_theta && !is_const;
        }
        assert!(has_const, "null space must include the constant 1");
        assert!(
            has_y,
            "null space must include the nonperiodic coordinate y (gam#1423)"
        );
        assert!(
            !depends_on_theta,
            "periodic axis θ must contribute only the constant, never a linear θ column"
        );
    }

    #[test]
    fn cylinder_linear_in_y_is_unpenalised_gam1423() {
        // gam#1423: build the center Gram K_CC and the non-periodic null-space
        // projector Z, then confirm the linear-in-y direction lies in the kernel
        // null space — i.e. f(θ, y) = a + b·y has EXACTLY zero penalty energy.
        let two_pi = std::f64::consts::TAU;
        let theta = [0.0, 0.7, 1.5, 2.3, 3.0, 3.9, 4.6, 5.4, two_pi];
        let y = [0.0, 0.2, 0.45, 0.6, 0.3, 0.8, 0.95, 0.1, 0.5];
        let mut centers = Array2::<f64>::zeros((theta.len(), 2));
        for i in 0..theta.len() {
            centers[[i, 0]] = theta[i];
            centers[[i, 1]] = y[i];
        }
        let periodic_per_axis = [true, false];
        let periods = [two_pi, 1.0];
        let axis_bounds = mixed_periodicity_axis_bounds(centers.view(), &periodic_per_axis);
        let k = centers.nrows();
        let m = 2usize;

        // Z = null(Pᵀ) for P = [1, y].
        let poly_block =
            mixed_periodicity_nullspace_poly_block(centers.view(), m, &periodic_per_axis);
        let z = kernel_constraint_nullspace_from_matrix(poly_block.view())
            .expect("null-space basis must build");

        // K_CC = additive kernel at center pairs.
        let mut k_cc = Array2::<f64>::zeros((k, k));
        fill_symmetric_from_row_kernel(&mut k_cc, |i, j| {
            mixed_periodicity_additive_kernel(
                centers.row(i),
                centers.row(j),
                m,
                &periodic_per_axis,
                &periods,
                &axis_bounds,
            )
        })
        .expect("center kernel must build");

        // Ω = Zᵀ K_CC Z is the realized penalty in the kernel-coefficient frame.
        let omega = fast_ab(&fast_atb(&z, &k_cc), &z);
        let sym = symmetrize(&omega);
        let (evals, _) = FaerEigh::eigh(&sym, Side::Lower).expect("eigh");
        let lambda_min = evals.iter().copied().fold(f64::INFINITY, f64::min);
        assert!(
            lambda_min > -1e-8,
            "Ω = ZᵀK_CC Z must be PSD; λ_min = {lambda_min:.3e}"
        );

        // The linear-in-y trend is carried by the explicit polynomial columns
        // (which receive a zero penalty block), so the kernel-coefficient penalty
        // Ω never sees it. Confirm directly that evaluating the penalty quadratic
        // form on the constant and linear-in-y design directions yields zero: any
        // design coefficient vector that is purely in the polynomial block has
        // zero penalty because Ω only occupies the kernel block. We model that by
        // checking that K_CC applied to the {1, y} span sits inside the column
        // space the polynomial block already spans — i.e. Pᵀ K_CC P has no energy
        // that Z can pick up. Concretely, Z is orthogonal to {1, y} by
        // construction, so Zᵀ·(linear-in-y vector) = 0.
        let ones: Vec<f64> = vec![1.0; k];
        let yvec: Vec<f64> = (0..k).map(|i| centers[[i, 1]]).collect();
        for (label, v) in [("constant", &ones), ("linear-in-y", &yvec)] {
            let mut proj = vec![0.0f64; z.ncols()];
            for c in 0..z.ncols() {
                let mut acc = 0.0;
                for r in 0..k {
                    acc += z[[r, c]] * v[r];
                }
                proj[c] = acc;
            }
            let norm: f64 = proj.iter().map(|p| p * p).sum::<f64>().sqrt();
            assert!(
                norm < 1e-9,
                "{label} direction must lie in the unpenalised null space \
                 (Zᵀv = 0); got ‖Zᵀv‖ = {norm:.3e} (gam#1423)"
            );
        }
    }
}

#[cfg(test)]
mod hybrid_high_dim_psd_tests {
    //! Regression tests for gam#1424: the high-dimensional hybrid
    //! (Duchon–Matérn) constrained bending penalty must be positive
    //! semidefinite. Before the fix, the kernel was assembled from an
    //! alternating partial-fraction sum whose individually-huge polyharmonic /
    //! Matérn blocks cancelled to the float noise floor; for d=16, p=2, s=7 the
    //! constrained spectrum was λ_min ≈ −0.264 after normalization. The
    //! cancellation-free single-integral kernel evaluation
    //! (`duchon_hybrid_kernel_stable_integral`) restores genuine PSD-ness.
    use super::*;
    use faer::Side;
    use gam_linalg::faer_ndarray::FaerEigh;
    use ndarray::Array2;

    /// Deterministic pseudo-random centers in `[-1, 1]^d` (no RNG dependency in
    /// the core crate). A 64-bit SplitMix-style generator seeded from `(d, k)`.
    fn deterministic_centers(d: usize, k: usize) -> Array2<f64> {
        let mut state: u64 = 0x9E37_79B9_7F4A_7C15u64
            .wrapping_mul(d as u64 + 1)
            .wrapping_add(k as u64);
        let mut next = || {
            state = state.wrapping_add(0x9E37_79B9_7F4A_7C15);
            let mut z = state;
            z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
            z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
            z ^= z >> 31;
            // Map to (−1, 1).
            (z as f64 / u64::MAX as f64) * 2.0 - 1.0
        };
        let mut c = Array2::<f64>::zeros((k, d));
        for i in 0..k {
            for j in 0..d {
                c[[i, j]] = next();
            }
        }
        c
    }

    /// Constrained bending-penalty spectrum λ_min for a hybrid Duchon smooth at
    /// the resolved (p, s) for the given dimension. Mirrors the construction in
    /// `duchon_native_penalty_candidates` (Z = null(Pᵀ), Ω = α²·ZᵀK_CC Z) but
    /// stops at the eigenvalues so the test can assert PSD-ness directly.
    fn hybrid_constrained_lambda_min(d: usize, s_order: usize) -> f64 {
        let centers = deterministic_centers(d, 4 * d);
        let nullspace = DuchonNullspaceOrder::Linear; // m = 2 ⇒ p = 2.
        let effective = duchon_effective_nullspace_order(centers.view(), nullspace);
        let poly_block = polynomial_block_from_order(centers.view(), effective);
        let z = kernel_constraint_nullspace_from_matrix(poly_block.view())
            .expect("kernel null-space basis must build");
        let omega = duchon_constrained_bending_penalty(
            centers.view(),
            Some(1.0),
            s_order as f64,
            effective,
            None,
            &z,
        )
        .expect("hybrid constrained bending penalty must build and pass the PSD check");
        let sym = symmetrize(&omega);
        let (evals, _) = FaerEigh::eigh(&sym, Side::Lower).expect("eigh");
        let max_abs = evals.iter().copied().fold(0.0_f64, |a, v| a.max(v.abs()));
        let lambda_min = evals.iter().copied().fold(f64::INFINITY, f64::min);
        // Report in spectral (scale-relative) units so the assertion is
        // invariant to the penalty's overall magnitude.
        lambda_min / max_abs.max(f64::MIN_POSITIVE)
    }

    #[test]
    fn hybrid_d16_m2_s7_constrained_spectrum_is_psd_gam1424() {
        // The exact spectral density ρ^{-2p}(κ²+ρ²)^{-s} is nonnegative, so the
        // constrained Gram is PSD in exact arithmetic. Before gam#1424 the
        // partial-fraction kernel evaluation lost every significant digit here
        // and the normalized λ_min was ≈ −0.264. The stable single-integral
        // kernel keeps λ_min at the float noise floor. Tolerance is the same
        // scale-relative noise floor the penalty pipeline scores this block
        // against, taken from its own constant; do NOT weaken it.
        let d = 16;
        let n = 4 * d; // penalty dimension upper bound (kernel coeff frame).
        let tol = (n as f64) * SPECTRAL_RANK_RELATIVE_TOLERANCE;
        let lambda_min_rel = hybrid_constrained_lambda_min(d, 7);
        assert!(
            lambda_min_rel >= -tol,
            "d=16, m=2, s=7 hybrid constrained spectrum not PSD: \
             λ_min/λ_max = {lambda_min_rel:.6e} (tol = −{tol:.6e}) (gam#1424)"
        );
    }

    #[test]
    fn hybrid_other_high_dims_constrained_spectrum_is_psd_gam1424() {
        // A spread of high-d hybrid orders that also lost positive-definiteness
        // through partial-fraction cancellation (d=8: ~7 digits lost; d=12: ~13;
        // d=10). Each must now be PSD to the float noise floor.
        for (d, s) in [(8usize, 3usize), (10, 4), (12, 5)] {
            let n = 4 * d;
            let tol = (n as f64) * SPECTRAL_RANK_RELATIVE_TOLERANCE;
            let lambda_min_rel = hybrid_constrained_lambda_min(d, s);
            assert!(
                lambda_min_rel >= -tol,
                "d={d}, m=2, s={s} hybrid constrained spectrum not PSD: \
                 λ_min/λ_max = {lambda_min_rel:.6e} (tol = −{tol:.6e}) (gam#1424)"
            );
        }
    }
}