fdars_core/

smooth_basis.rs

//! Basis-penalized smoothing with continuous derivative penalties.
//!
//! This module implements `smooth.basis` from R's fda package. Unlike the
//! discrete difference penalty used in P-splines (`basis.rs`), this uses
//! continuous derivative penalties: `min ||y - Φc||² + λ·∫(Lf)² dt`.
//!
//! Key capabilities:
//! - [`smooth_basis`] — Penalized least squares with continuous roughness penalty
//! - [`smooth_basis_gcv`] — GCV-optimal smoothing parameter selection
//! - [`bspline_penalty_matrix`] / [`fourier_penalty_matrix`] — Roughness penalty matrices

use crate::basis::{bspline_basis, fourier_basis_with_period};
use crate::helpers::simpsons_weights;
use crate::matrix::FdMatrix;
use nalgebra::DMatrix;
use std::f64::consts::PI;

// ─── Types ──────────────────────────────────────────────────────────────────

/// Basis type for penalized smoothing.
#[non_exhaustive]
#[derive(Debug, Clone, PartialEq)]
pub enum BasisType {
    /// B-spline basis with given order (typically 4 for cubic).
    Bspline { order: usize },
    /// Fourier basis with given period.
    Fourier { period: f64 },
}

/// Functional data parameter object (basis + penalty specification).
#[derive(Debug, Clone, PartialEq)]
pub struct FdPar {
    /// Type of basis system.
    pub basis_type: BasisType,
    /// Number of basis functions.
    pub nbasis: usize,
    /// Smoothing parameter.
    pub lambda: f64,
    /// Derivative order for the penalty (default: 2).
    pub lfd_order: usize,
    /// Precomputed K×K penalty matrix (column-major).
    pub penalty_matrix: Vec<f64>,
}

/// Result of log-domain positive smoothing.
///
/// The fitted values are guaranteed strictly positive (> 0) for any strictly-positive
/// input, because they are obtained by exponentiating the log-domain B-spline smooth.
///
/// # Bias caveat
///
/// This is a log-transform smoother: it minimises squared error in **log space**
/// (`‖log(y) − Φc‖²`), not in the original scale.  The back-transformed values
/// `exp(Φ̂c)` are therefore **not** the minimum-MSE estimate of `E[y | t]` — they
/// tend to underestimate the conditional mean by a Jensen's-inequality correction.
/// For an approximately unbiased estimate in the original scale apply
/// `exp(fitted + σ² / 2)` where `σ²` is the conditional variance of the log-transformed
/// fit (e.g. from the GCV score).  The correction is intentionally left to the caller
/// in v1.
#[derive(Debug, Clone, PartialEq)]
#[non_exhaustive]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct SmoothPositiveResult {
    /// Fitted values in the original (positive) scale (n × m, all entries > 0).
    pub fitted: FdMatrix,
    /// B-spline coefficients from the log-domain smooth (n × K).
    pub log_coefficients: FdMatrix,
    /// Effective degrees of freedom of the log-domain fit.
    pub edf: f64,
    /// GCV score of the log-domain fit.
    pub gcv: f64,
}

/// Result of basis-penalized smoothing.
#[derive(Debug, Clone, PartialEq)]
#[non_exhaustive]
pub struct SmoothBasisResult {
    /// Basis coefficients (n × K).
    pub coefficients: FdMatrix,
    /// Fitted values (n × m).
    pub fitted: FdMatrix,
    /// Effective degrees of freedom.
    pub edf: f64,
    /// Generalized cross-validation score.
    pub gcv: f64,
    /// AIC.
    pub aic: f64,
    /// BIC.
    pub bic: f64,
    /// Roughness penalty matrix (K × K, column-major).
    pub penalty_matrix: Vec<f64>,
    /// Number of basis functions used.
    pub nbasis: usize,
}

// ─── Penalty Matrices ───────────────────────────────────────────────────────

/// Compute the roughness penalty matrix for B-splines via numerical quadrature.
///
/// R\[j,k\] = ∫ D^m B_j(t) · D^m B_k(t) dt
///
/// Uses Simpson's rule on a fine sub-grid for each knot interval.
///
/// # Arguments
/// * `argvals` — Evaluation points (length m)
/// * `nbasis` — Number of basis functions
/// * `order` — B-spline order (typically 4 for cubic)
/// * `lfd_order` — Derivative order for penalty (typically 2)
///
/// # Returns
/// K × K penalty matrix in column-major layout (K = nbasis)
pub fn bspline_penalty_matrix(
    argvals: &[f64],
    nbasis: usize,
    order: usize,
    lfd_order: usize,
) -> Vec<f64> {
    if nbasis < 2 || order < 1 || lfd_order >= order || argvals.len() < 2 {
        return vec![0.0; nbasis * nbasis];
    }

    let nknots = nbasis.saturating_sub(order).max(2);

    // Create a fine quadrature grid (10 sub-points per original interval)
    let n_sub = 10;
    let t_min = argvals[0];
    let t_max = argvals[argvals.len() - 1];
    let n_quad = (argvals.len() - 1) * n_sub + 1;
    let quad_t: Vec<f64> = (0..n_quad)
        .map(|i| t_min + (t_max - t_min) * i as f64 / (n_quad - 1) as f64)
        .collect();

    // Evaluate B-spline basis on fine grid
    let basis_fine = bspline_basis(&quad_t, nknots, order);
    let actual_nbasis = basis_fine.len() / n_quad;

    // Compute derivatives of B-spline basis numerically
    let h = (t_max - t_min) / (n_quad - 1) as f64;
    let deriv_basis = differentiate_basis_columns(&basis_fine, n_quad, actual_nbasis, h, lfd_order);

    // Integration weights on fine grid
    let weights = simpsons_weights(&quad_t);

    // Compute penalty matrix: R[j,k] = ∫ D^m B_j · D^m B_k dt
    integrate_symmetric_penalty(&deriv_basis, &weights, actual_nbasis, n_quad)
}

/// Compute the roughness penalty matrix for a Fourier basis.
///
/// For Fourier basis, the penalty is diagonal with eigenvalues `(2πk/T)^(2m)`.
///
/// # Arguments
/// * `nbasis` — Number of basis functions
/// * `period` — Period of the Fourier basis
/// * `lfd_order` — Derivative order for penalty
///
/// # Returns
/// K × K penalty matrix in column-major layout
pub fn fourier_penalty_matrix(nbasis: usize, period: f64, lfd_order: usize) -> Vec<f64> {
    let k = nbasis;
    let mut penalty = vec![0.0; k * k];

    // First basis function is constant → lfd_order-th derivative is 0
    // penalty[0] = 0 (already zero)

    // For sin/cos pairs: eigenvalue is (2πk/T)^(2m)
    // Matches R's fda package convention (sqrt(2)-normalized basis)
    let mut freq = 1;
    let mut idx = 1;
    while idx < k {
        let omega = 2.0 * PI * f64::from(freq) / period;
        let eigenval = omega.powi(2 * lfd_order as i32);

        // sin component
        if idx < k {
            penalty[idx + idx * k] = eigenval;
            idx += 1;
        }
        // cos component
        if idx < k {
            penalty[idx + idx * k] = eigenval;
            idx += 1;
        }
        freq += 1;
    }

    penalty
}

// ─── Smoothing Functions ────────────────────────────────────────────────────

/// Perform basis-penalized smoothing.
///
/// Solves `(Φ'Φ + λR)c = Φ'y` per curve via Cholesky decomposition.
/// This implements `smooth.basis` from R's fda package.
///
/// # Arguments
/// * `data` — Functional data matrix (n × m)
/// * `argvals` — Evaluation points (length m)
/// * `fdpar` — Functional parameter object specifying basis and penalty
///
/// # Returns
/// [`SmoothBasisResult`] with coefficients, fitted values, and diagnostics.
pub fn smooth_basis(
    data: &FdMatrix,
    argvals: &[f64],
    fdpar: &FdPar,
) -> Result<SmoothBasisResult, crate::FdarError> {
    let (n, m) = data.shape();
    if n == 0 || m == 0 || argvals.len() != m || fdpar.nbasis < 2 {
        return Err(crate::FdarError::InvalidDimension {
            parameter: "data/argvals/fdpar",
            expected: "n > 0, m > 0, argvals.len() == m, nbasis >= 2".to_string(),
            actual: format!(
                "n={}, m={}, argvals.len()={}, nbasis={}",
                n,
                m,
                argvals.len(),
                fdpar.nbasis
            ),
        });
    }

    // Evaluate basis on argvals
    let (basis_flat, actual_nbasis) = evaluate_basis(argvals, &fdpar.basis_type, fdpar.nbasis);
    let k = actual_nbasis;

    let b_mat = DMatrix::from_column_slice(m, k, &basis_flat);
    let r_mat = DMatrix::from_column_slice(k, k, &fdpar.penalty_matrix);

    // (Φ'Φ + λR + εI) — small ridge ensures positive definiteness
    let btb = b_mat.transpose() * &b_mat;
    let ridge_eps = 1e-10;
    let system: DMatrix<f64> =
        &btb + fdpar.lambda * &r_mat + ridge_eps * DMatrix::<f64>::identity(k, k);

    // Invert the penalized system
    let system_inv =
        invert_penalized_system(&system, k).ok_or_else(|| crate::FdarError::ComputationFailed {
            operation: "matrix inversion",
            detail: "failed to invert penalized system (Φ'Φ + λR); try increasing lambda or reducing the number of basis functions".to_string(),
        })?;

    // Hat matrix: H = Φ (Φ'Φ + λR)^{-1} Φ'  →  EDF = tr(H)
    let h_mat = &b_mat * &system_inv * b_mat.transpose();
    let edf: f64 = (0..m).map(|i| h_mat[(i, i)]).sum();

    // Project all curves
    let proj = &system_inv * b_mat.transpose();
    let (all_coefs, all_fitted, total_rss) = project_all_curves(data, &b_mat, &proj, n, m, k);

    let total_points = (n * m) as f64;
    let gcv = compute_gcv(total_rss, total_points, edf, m);
    let mse = total_rss / total_points;
    // Total effective degrees of freedom = n curves * per-curve edf
    let total_edf = n as f64 * edf;
    let aic = total_points * mse.max(1e-300).ln() + 2.0 * total_edf;
    let bic = total_points * mse.max(1e-300).ln() + total_points.ln() * total_edf;

    Ok(SmoothBasisResult {
        coefficients: all_coefs,
        fitted: all_fitted,
        edf,
        gcv,
        aic,
        bic,
        penalty_matrix: fdpar.penalty_matrix.clone(),
        nbasis: k,
    })
}

/// Perform basis-penalized smoothing with GCV-optimal lambda.
///
/// Searches over a log-lambda grid and selects the lambda minimizing GCV.
///
/// # Arguments
/// * `data` — Functional data matrix (n × m)
/// * `argvals` — Evaluation points (length m)
/// * `basis_type` — Type of basis system
/// * `nbasis` — Number of basis functions
/// * `lfd_order` — Derivative order for penalty
/// * `log_lambda_range` — Range of log10(lambda) to search, e.g. (-8.0, 4.0)
/// * `n_grid` — Number of grid points for the search
pub fn smooth_basis_gcv(
    data: &FdMatrix,
    argvals: &[f64],
    basis_type: &BasisType,
    nbasis: usize,
    lfd_order: usize,
    log_lambda_range: (f64, f64),
    n_grid: usize,
) -> Option<SmoothBasisResult> {
    let m = argvals.len();
    if m == 0 || nbasis < 2 || n_grid < 2 {
        return None;
    }

    // Compute penalty matrix once
    let penalty = match basis_type {
        BasisType::Bspline { order } => bspline_penalty_matrix(argvals, nbasis, *order, lfd_order),
        BasisType::Fourier { period } => fourier_penalty_matrix(nbasis, *period, lfd_order),
    };

    let (lo, hi) = log_lambda_range;
    let mut best_gcv = f64::INFINITY;
    let mut best_result: Option<SmoothBasisResult> = None;

    for i in 0..n_grid {
        let log_lam = lo + (hi - lo) * i as f64 / (n_grid - 1) as f64;
        let lam = 10.0_f64.powf(log_lam);

        let fdpar = FdPar {
            basis_type: basis_type.clone(),
            nbasis,
            lambda: lam,
            lfd_order,
            penalty_matrix: penalty.clone(),
        };

        if let Ok(result) = smooth_basis(data, argvals, &fdpar) {
            if result.gcv < best_gcv {
                best_gcv = result.gcv;
                best_result = Some(result);
            }
        }
    }

    best_result
}

/// Perform basis-penalized smoothing with AIC-optimal lambda.
///
/// Mirrors [`smooth_basis_gcv`] but selects the lambda minimizing the `aic`
/// field of [`SmoothBasisResult`] (which `smooth_basis` computes as
/// `total_points·ln(mse) + 2·total_edf`, using the hat-matrix trace as the
/// degrees of freedom — the same df-as-trace convention as the kernel-path
/// [`aic_smoother`](crate::smoothing::aic_smoother)). It searches the identical
/// log-lambda grid, so it is a drop-in AIC alternative that leaves
/// [`smooth_basis_gcv`] untouched.
///
/// # Arguments
/// * `data` — Functional data matrix (n × m)
/// * `argvals` — Evaluation points (length m)
/// * `basis_type` — Type of basis system
/// * `nbasis` — Number of basis functions
/// * `lfd_order` — Derivative order for penalty
/// * `log_lambda_range` — Range of log10(lambda) to search, e.g. (-8.0, 4.0)
/// * `n_grid` — Number of grid points for the search
///
/// # Example
/// ```no_run
/// use fdars_core::smooth_basis::{smooth_basis_aic, BasisType};
/// use fdars_core::matrix::FdMatrix;
///
/// let argvals: Vec<f64> = (0..20).map(|i| i as f64 / 19.0).collect();
/// let data = FdMatrix::from_column_major(vec![0.0; 20], 1, 20).unwrap();
/// let result = smooth_basis_aic(
///     &data,
///     &argvals,
///     &BasisType::Bspline { order: 4 },
///     8,
///     2,
///     (-8.0, 4.0),
///     30,
/// );
/// ```
pub fn smooth_basis_aic(
    data: &FdMatrix,
    argvals: &[f64],
    basis_type: &BasisType,
    nbasis: usize,
    lfd_order: usize,
    log_lambda_range: (f64, f64),
    n_grid: usize,
) -> Option<SmoothBasisResult> {
    let m = argvals.len();
    if m == 0 || nbasis < 2 || n_grid < 2 {
        return None;
    }

    // Compute penalty matrix once (same as smooth_basis_gcv).
    let penalty = match basis_type {
        BasisType::Bspline { order } => bspline_penalty_matrix(argvals, nbasis, *order, lfd_order),
        BasisType::Fourier { period } => fourier_penalty_matrix(nbasis, *period, lfd_order),
    };

    let (lo, hi) = log_lambda_range;
    let mut best_aic = f64::INFINITY;
    let mut best_result: Option<SmoothBasisResult> = None;

    for i in 0..n_grid {
        let log_lam = lo + (hi - lo) * i as f64 / (n_grid - 1) as f64;
        let lam = 10.0_f64.powf(log_lam);

        let fdpar = FdPar {
            basis_type: basis_type.clone(),
            nbasis,
            lambda: lam,
            lfd_order,
            penalty_matrix: penalty.clone(),
        };

        if let Ok(result) = smooth_basis(data, argvals, &fdpar) {
            // Reuse the aic field smooth_basis already populates; do NOT
            // recompute AIC inline, keeping the formula consistent.
            if result.aic < best_aic {
                best_aic = result.aic;
                best_result = Some(result);
            }
        }
    }

    best_result
}

/// Positive-valued smoothing via log-domain transformation.
///
/// Smooths strictly-positive functional data by transforming to the log scale,
/// applying the existing [`smooth_basis`] B-spline smoother, and then
/// exponentiating the fitted values so the result is guaranteed strictly positive.
///
/// This is the log-domain positive-smoothing idiom; it is appropriate for responses
/// that are nonnegative by construction (densities, concentrations, intensities).
///
/// # Arguments
/// * `data` — Functional data matrix (n × m); every element must be **strictly positive**
/// * `argvals` — Evaluation points (length m)
/// * `fdpar` — Functional parameter object specifying basis and penalty (unchanged from
///   [`smooth_basis`])
///
/// # Errors
/// Returns [`crate::FdarError::InvalidDimension`] if dimensions are inconsistent or
/// `n == 0 || m == 0`.
/// Returns [`crate::FdarError::InvalidParameter`] if any element of `data` is ≤ 0 (the
/// log transform is undefined there).
/// Propagates [`crate::FdarError::ComputationFailed`] from the inner [`smooth_basis`]
/// call if matrix inversion fails.
///
/// # Bias caveat
///
/// See [`SmoothPositiveResult`] for the retransformation bias note.
///
/// # R baseline
/// Analogous to `fda::smooth.pos` (R fda package), but uses the existing fdars
/// B-spline penalty infrastructure rather than a separate estimation path.  The
/// key divergence: this implementation minimises squared error in log space rather
/// than iterating toward the conditional mean in the original scale.
///
/// # Example
/// ```no_run
/// use fdars_core::smooth_basis::{smooth_positive, BasisType, FdPar, bspline_penalty_matrix};
/// use fdars_core::matrix::FdMatrix;
///
/// let argvals: Vec<f64> = (0..41).map(|i| i as f64 / 40.0).collect();
/// let y: Vec<f64> = argvals.iter().map(|&t| 2.0 + t.sin() + 0.05).collect();
/// let data = FdMatrix::from_column_major(y, 1, 41).unwrap();
/// let penalty = bspline_penalty_matrix(&argvals, 10, 4, 2);
/// let fdpar = FdPar { basis_type: BasisType::Bspline { order: 4 }, nbasis: 10, lambda: 1e-3, lfd_order: 2, penalty_matrix: penalty };
/// let result = smooth_positive(&data, &argvals, &fdpar).unwrap();
/// assert!(result.fitted[(0, 0)] > 0.0);
/// ```
#[must_use = "expensive computation whose result should not be discarded"]
pub fn smooth_positive(
    data: &FdMatrix,
    argvals: &[f64],
    fdpar: &FdPar,
) -> Result<SmoothPositiveResult, crate::FdarError> {
    let (n, m) = data.shape();
    if n == 0 || m == 0 || argvals.len() != m {
        return Err(crate::FdarError::InvalidDimension {
            parameter: "data/argvals",
            expected: "n > 0, m > 0, argvals.len() == m".to_string(),
            actual: format!("n={}, m={}, argvals.len()={}", n, m, argvals.len()),
        });
    }

    // T-44-07: validate all data > 0 before any ln to prevent NaN/-Inf.
    for i in 0..n {
        for j in 0..m {
            if data[(i, j)] <= 0.0 {
                return Err(crate::FdarError::InvalidParameter {
                    parameter: "data",
                    message: format!(
                        "smooth_positive requires strictly positive data (log-domain smoother); \
                         found value {} <= 0 at observation {}, evaluation point {}",
                        data[(i, j)],
                        i,
                        j
                    ),
                });
            }
        }
    }

    // Build log-transformed data matrix.
    let mut log_data = FdMatrix::zeros(n, m);
    for i in 0..n {
        for j in 0..m {
            log_data[(i, j)] = data[(i, j)].ln();
        }
    }

    // Delegate to the existing smooth_basis on log(data).
    let inner = smooth_basis(&log_data, argvals, fdpar)?;

    // Exp-reconstruct: exp of any finite real is strictly positive.
    let mut fitted = FdMatrix::zeros(n, m);
    for i in 0..n {
        for j in 0..m {
            fitted[(i, j)] = inner.fitted[(i, j)].exp();
        }
    }

    Ok(SmoothPositiveResult {
        fitted,
        log_coefficients: inner.coefficients,
        edf: inner.edf,
        gcv: inner.gcv,
    })
}

/// Result of Ramsay monotone smoothing via integral-of-exponential representation.
///
/// The fitted values are structurally monotone: since `f'(t) = β₁·exp(w(t))` and
/// `exp(w(t)) > 0` always, the sign of `β₁` determines direction.  This guarantee
/// holds even if the Gauss-Newton iteration did not fully converge.
///
/// # R baseline
/// Analogous to `fda::smooth.monotone` (R fda package).
#[derive(Debug, Clone, PartialEq)]
#[non_exhaustive]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct SmoothMonotoneResult {
    /// Fitted values at the input `argvals` (length m). Structurally monotone.
    pub fitted: Vec<f64>,
    /// Intercept β₀ (value of the fitted function at `argvals[0]`).
    pub beta0: f64,
    /// Scale β₁. Positive → nondecreasing fit; negative → nonincreasing fit.
    pub beta1: f64,
    /// B-spline coefficients for the log-rate function `w(u)` (length = actual nbasis).
    pub w_coefficients: Vec<f64>,
    /// Number of Gauss-Newton iterations executed.
    pub iterations: usize,
    /// Whether the Gauss-Newton update norm fell below `1e-8` before `max_iter`.
    pub converged: bool,
}

/// Monotone smoothing via Ramsay's integral-of-exponential representation.
///
/// Fits `f(t) = β₀ + β₁ ∫₀ᵗ exp(w(u)) du` where `w(u) = Σⱼ αⱼ Ψⱼ(u)` is a
/// B-spline expansion, and `(β₀, β₁, α)` are estimated by penalised Gauss-Newton
/// nonlinear least squares.
///
/// **Monotonicity is structural**: `f'(t) = β₁·exp(w(t)) > 0` (or `< 0` if
/// `β₁ < 0`) for any parameter values, including underconverged iterates.
/// Direction is auto-detected from the data: increasing data → `β₁ > 0`,
/// decreasing data → `β₁ < 0`.
///
/// # Arguments
/// * `data` — Observed response values (length m, need not be monotone)
/// * `argvals` — Evaluation / domain points (length m, strictly increasing recommended)
/// * `nbasis` — Target number of B-spline basis functions for `w(u)` (≥ 2)
/// * `order` — B-spline order (e.g. 4 for cubic; ≥ 1)
/// * `lambda` — Roughness penalty weight on `w(u)` (≥ 0; 0 = no penalty)
/// * `max_iter` — Maximum Gauss-Newton iterations (≥ 1; typically 50 suffices)
///
/// # Errors
/// * [`crate::FdarError::InvalidDimension`] — `data.len() < 3`, `argvals.len() != data.len()`
/// * [`crate::FdarError::InvalidParameter`] — `nbasis < 2`, `order < 1`, `lambda < 0`,
///   `max_iter < 1`
/// * [`crate::FdarError::ComputationFailed`] — Cholesky of the GN normal equations failed
///   after Levenberg ridge augmentation (rank-deficient system)
///
/// # Divergence from R
/// Dense Gauss-Newton; no sparse solver.  Practical for m ≤ 500.  The
/// `λ_mono` parameter penalises roughness of `w(u)` via `bspline_penalty_matrix`,
/// reusing the existing fdars infrastructure.  Direction flip (decreasing data)
/// is handled by allowing `β₁ < 0` rather than by negating the data.
///
/// # Example
/// ```no_run
/// use fdars_core::smooth_basis::smooth_monotone;
///
/// let t: Vec<f64> = (0..41).map(|i| i as f64 / 40.0).collect();
/// let y: Vec<f64> = t.iter().map(|&ti| ti * ti).collect(); // t² is increasing
/// let result = smooth_monotone(&y, &t, 8, 4, 1e-3, 50).unwrap();
/// assert!(result.beta1 > 0.0);
/// for i in 1..t.len() {
///     assert!(result.fitted[i] >= result.fitted[i - 1] - 1e-9);
/// }
/// ```
#[must_use = "expensive nonlinear least-squares computation whose result should not be discarded"]
pub fn smooth_monotone(
    data: &[f64],
    argvals: &[f64],
    nbasis: usize,
    order: usize,
    lambda: f64,
    max_iter: usize,
) -> Result<SmoothMonotoneResult, crate::FdarError> {
    // ── Validation ────────────────────────────────────────────────────────────
    let m = data.len();
    if m < 3 || argvals.len() != m {
        return Err(crate::FdarError::InvalidDimension {
            parameter: "data/argvals",
            expected: "data.len() >= 3, argvals.len() == data.len()".to_string(),
            actual: format!("data.len()={}, argvals.len()={}", m, argvals.len()),
        });
    }
    if nbasis < 2 {
        return Err(crate::FdarError::InvalidParameter {
            parameter: "nbasis",
            message: format!("smooth_monotone requires nbasis >= 2, got {}", nbasis),
        });
    }
    if order < 1 {
        return Err(crate::FdarError::InvalidParameter {
            parameter: "order",
            message: format!("smooth_monotone requires order >= 1, got {}", order),
        });
    }
    if lambda < 0.0 {
        return Err(crate::FdarError::InvalidParameter {
            parameter: "lambda",
            message: format!("smooth_monotone requires lambda >= 0, got {}", lambda),
        });
    }
    if max_iter < 1 {
        return Err(crate::FdarError::InvalidParameter {
            parameter: "max_iter",
            message: format!("smooth_monotone requires max_iter >= 1, got {}", max_iter),
        });
    }

    // ── B-spline basis setup ──────────────────────────────────────────────────
    // Use the same nknots convention as bspline_penalty_matrix so that
    // actual_k == penalty matrix dimension.
    let nknots = nbasis.saturating_sub(order).max(2);
    // basis is n*actual_k column-major: basis[ti + j*m] = Ψ_j(t_i)
    let basis_flat = crate::basis::bspline_basis(argvals, nknots, order);
    let actual_k = basis_flat.len() / m;

    // Precompute roughness penalty R (actual_k × actual_k, COLUMN-MAJOR).
    // We will need it in ROW-MAJOR for the normal equations; R is symmetric so
    // penalty_col[a + b*k] == penalty_col[b + a*k] — no explicit transpose needed.
    let penalty_col = bspline_penalty_matrix(argvals, nbasis, order, 2);
    // penalty_col may be sized (actual_k_pm)² — compute its actual side length.
    let actual_k_pm = (penalty_col.len() as f64).sqrt() as usize;
    // k used throughout: take the minimum in case of rare sizing mismatch.
    let k = actual_k.min(actual_k_pm);

    // Helper: Ψ_j(t_i) from column-major basis_flat.
    let psi_at = |i: usize, j: usize| -> f64 {
        if j < actual_k {
            basis_flat[i + j * m]
        } else {
            0.0
        }
    };

    // ── Initialisation ────────────────────────────────────────────────────────
    let t_range = (argvals[m - 1] - argvals[0]).max(1e-12);
    let beta0_init = data[0];
    // beta1 sign auto-detects data direction; avoids dividing by zero t_range.
    let beta1_init = (data[m - 1] - data[0]) / t_range;
    // If the data is flat, use a tiny positive value so the fit is still valid.
    let beta1_init = if beta1_init.abs() < 1e-12 {
        1e-6
    } else {
        beta1_init
    };
    let mut beta0 = beta0_init;
    let mut beta1 = beta1_init;
    let mut alpha = vec![0.0_f64; k]; // w = 0 everywhere → exp(w) = 1

    let big_p = 2 + k; // total parameter dimension

    // ── Helper: build W[i] and Iexp_psi[i*k + j] via cumulative trapezoid ────
    // argvals may not start at 0; we shift so integration starts at argvals[0].
    // W[i] = ∫_{argvals[0]}^{argvals[i]} exp(w(u)) du (trapezoid rule on data grid)
    // Iexp_psi[i*k+j] = ∫_{argvals[0]}^{argvals[i]} exp(w(u))·Ψ_j(u) du
    let build_integrals = |alpha: &[f64]| -> (Vec<f64>, Vec<f64>) {
        let mut exp_w = vec![0.0_f64; m];
        for i in 0..m {
            let w_i: f64 = (0..k).map(|j| alpha[j] * psi_at(i, j)).sum();
            // Clamp to [-30, 30] to avoid overflow (threat T-44-08).
            let w_clamped = w_i.clamp(-30.0, 30.0);
            exp_w[i] = w_clamped.exp();
        }
        // W[i]: cumulative trapezoid of exp_w.
        let mut w_int = vec![0.0_f64; m];
        for i in 1..m {
            let dt = argvals[i] - argvals[i - 1];
            w_int[i] = w_int[i - 1] + 0.5 * (exp_w[i - 1] + exp_w[i]) * dt;
        }
        // Iexp_psi: cumulative trapezoid of exp_w * Ψ_j for each j.
        let mut iexp_psi = vec![0.0_f64; m * k];
        for i in 1..m {
            let dt = argvals[i] - argvals[i - 1];
            for j in 0..k {
                let integrand_prev = exp_w[i - 1] * psi_at(i - 1, j);
                let integrand_cur = exp_w[i] * psi_at(i, j);
                iexp_psi[i * k + j] =
                    iexp_psi[(i - 1) * k + j] + 0.5 * (integrand_prev + integrand_cur) * dt;
            }
        }
        (w_int, iexp_psi)
    };

    // ── Gauss-Newton loop ─────────────────────────────────────────────────────
    let mut iterations = 0_usize;
    let mut converged = false;

    for iter in 0..max_iter {
        let (w_int, iexp_psi) = build_integrals(&alpha);

        // Residuals: r_i = data[i] - f(t_i).
        let mut resid = vec![0.0_f64; m];
        for i in 0..m {
            let f_i = beta0 + beta1 * w_int[i];
            resid[i] = data[i] - f_i;
        }

        // Build normal equations A = JᵀJ (row-major big_p × big_p) and g = Jᵀr.
        let mut a_mat = vec![0.0_f64; big_p * big_p];
        let mut g_vec = vec![0.0_f64; big_p];

        for i in 0..m {
            // Jacobian row for observation i: [1, W[i], β₁·Iexp_psi[i,0..k]]
            // col-0: ∂f/∂β₀ = 1
            // col-1: ∂f/∂β₁ = W[i]
            // col-(2+j): ∂f/∂α_j = β₁ · Iexp_psi[i*k+j]
            let j0 = 1.0_f64;
            let j1 = w_int[i];

            // Accumulate JᵀJ and Jᵀr.
            // Row/col indices: 0 = β₀, 1 = β₁, 2+j = α_j.

            // (0,0): j0*j0
            a_mat[0] += j0 * j0;
            // (0,1): j0*j1
            a_mat[1] += j0 * j1;
            a_mat[big_p] += j1 * j0; // (1,0) = symmetric
                                     // (1,1): j1*j1
            a_mat[big_p + 1] += j1 * j1;

            // g[0], g[1]
            g_vec[0] += j0 * resid[i];
            g_vec[1] += j1 * resid[i];

            for aj in 0..k {
                let jc = beta1 * iexp_psi[i * k + aj];
                // (0, 2+aj) and (2+aj, 0)
                a_mat[2 + aj] += j0 * jc;
                a_mat[(2 + aj) * big_p] += jc * j0;
                // (1, 2+aj) and (2+aj, 1)
                a_mat[big_p + (2 + aj)] += j1 * jc;
                a_mat[(2 + aj) * big_p + 1] += jc * j1;
                // (2+aj, 2+aj) and cross-terms
                for bj in 0..k {
                    let jd = beta1 * iexp_psi[i * k + bj];
                    a_mat[(2 + aj) * big_p + (2 + bj)] += jc * jd;
                }
                // g[2+aj]
                g_vec[2 + aj] += jc * resid[i];
            }
        }

        // Levenberg ridge on diagonal (threat T-44-09, T-44-10).
        for d in 0..big_p {
            a_mat[d * big_p + d] += 1e-6 * (1.0 + a_mat[d * big_p + d].abs());
        }

        // Add λ·R into the α–α sub-block.
        // penalty_col is COLUMN-MAJOR (k×k): penalty_col[a + b*k] = R[a,b].
        // a_mat is ROW-MAJOR (big_p×big_p): a_mat[(2+a)*big_p + (2+b)] += λ·R[a,b].
        for a in 0..k {
            for b in 0..k {
                let r_ab = if a < actual_k_pm && b < actual_k_pm {
                    penalty_col[a + b * actual_k_pm]
                } else {
                    0.0
                };
                a_mat[(2 + a) * big_p + (2 + b)] += lambda * r_ab;
            }
        }

        // Solve normal equations.
        let delta = crate::linalg::cholesky_solve(&a_mat, &g_vec, big_p)?;

        // Apply update.
        beta0 += delta[0];
        beta1 += delta[1];
        for j in 0..k {
            alpha[j] += delta[2 + j];
        }

        iterations = iter + 1;

        // Convergence check: ‖δ‖₂ < 1e-8.
        let delta_norm: f64 = delta.iter().map(|d| d * d).sum::<f64>().sqrt();
        if delta_norm < 1e-8 {
            converged = true;
            break;
        }
    }

    // ── Final fit ─────────────────────────────────────────────────────────────
    let (w_int_final, _) = build_integrals(&alpha);
    let fitted: Vec<f64> = (0..m).map(|i| beta0 + beta1 * w_int_final[i]).collect();

    Ok(SmoothMonotoneResult {
        fitted,
        beta0,
        beta1,
        w_coefficients: alpha,
        iterations,
        converged,
    })
}

// ─── Config Structs ─────────────────────────────────────────────────────────

/// Configuration for GCV-based smoothing parameter selection.
///
/// Collects all tuning parameters for [`smooth_basis_gcv_with_config`], with
/// sensible defaults obtained via [`SmoothBasisGcvConfig::default()`].
///
/// # Example
/// ```no_run
/// use fdars_core::smooth_basis::{SmoothBasisGcvConfig, BasisType};
///
/// let mut config = SmoothBasisGcvConfig::default();
/// config.nbasis = 20;
/// config.n_grid = 100;
/// ```
///
/// Construct via `SmoothBasisGcvConfig::default()`, then assign the fields you need (e.g. `let mut c = SmoothBasisGcvConfig::default(); c.field = …;`). This struct is `#[non_exhaustive]`, so external crates cannot build it with a struct literal — not even functional-update `..Default::default()` form.
#[non_exhaustive]
#[derive(Debug, Clone, PartialEq)]
pub struct SmoothBasisGcvConfig {
    /// Basis type (BSpline or Fourier).
    pub basis_type: BasisType,
    /// Number of basis functions (default: 15).
    pub nbasis: usize,
    /// Order of the roughness penalty differential operator (default: 2).
    pub lfd_order: usize,
    /// Range of log10(lambda) values to search (default: (-10.0, 2.0)).
    pub log_lambda_range: (f64, f64),
    /// Number of grid points in the lambda search (default: 50).
    pub n_grid: usize,
}

impl Default for SmoothBasisGcvConfig {
    fn default() -> Self {
        Self {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis: 15,
            lfd_order: 2,
            log_lambda_range: (-10.0, 2.0),
            n_grid: 50,
        }
    }
}

/// Perform basis-penalized smoothing with GCV-optimal lambda using a config struct.
///
/// This is the config-based alternative to [`smooth_basis_gcv`]. It takes data
/// parameters directly and reads all tuning parameters from the config.
///
/// # Arguments
/// * `data` — Functional data matrix (n × m)
/// * `argvals` — Evaluation points (length m)
/// * `config` — Tuning parameters
///
/// # Errors
///
/// Returns [`crate::FdarError::ComputationFailed`] if no valid smoothing result
/// is found for any lambda in the search grid.
#[must_use = "expensive computation whose result should not be discarded"]
pub fn smooth_basis_gcv_with_config(
    data: &FdMatrix,
    argvals: &[f64],
    config: &SmoothBasisGcvConfig,
) -> Result<SmoothBasisResult, crate::FdarError> {
    smooth_basis_gcv(
        data,
        argvals,
        &config.basis_type,
        config.nbasis,
        config.lfd_order,
        config.log_lambda_range,
        config.n_grid,
    )
    .ok_or_else(|| crate::FdarError::ComputationFailed {
        operation: "smooth_basis_gcv_with_config",
        detail: "no valid smoothing result found in GCV lambda search".to_string(),
    })
}

/// Configuration for cross-validation-based basis selection.
///
/// Collects all tuning parameters for [`basis_nbasis_cv_with_config`], with
/// sensible defaults obtained via [`BasisNbasisCvConfig::default()`].
///
/// # Example
/// ```no_run
/// use fdars_core::smooth_basis::{BasisNbasisCvConfig, BasisType, BasisCriterion};
///
/// let mut config = BasisNbasisCvConfig::default();
/// config.nbasis_range = (5, 25);
/// config.criterion = BasisCriterion::Aic;
/// ```
///
/// Construct via `BasisNbasisCvConfig::default()`, then assign the fields you need (e.g. `let mut c = BasisNbasisCvConfig::default(); c.field = …;`). This struct is `#[non_exhaustive]`, so external crates cannot build it with a struct literal — not even functional-update `..Default::default()` form.
#[non_exhaustive]
#[derive(Debug, Clone, PartialEq)]
pub struct BasisNbasisCvConfig {
    /// Basis type (default: BSpline with order 4).
    pub basis_type: BasisType,
    /// Range of nbasis values to try, inclusive (default: (5, 30)).
    pub nbasis_range: (usize, usize),
    /// Roughness penalty lambda (default: 1e-4).
    pub lambda: f64,
    /// Penalty order (default: 2).
    pub lfd_order: usize,
    /// Number of CV folds (default: 5). Only used when `criterion` is `Cv`.
    pub n_folds: usize,
    /// Selection criterion (default: `Gcv`).
    pub criterion: BasisCriterion,
}

impl Default for BasisNbasisCvConfig {
    fn default() -> Self {
        Self {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis_range: (5, 30),
            lambda: 1e-4,
            lfd_order: 2,
            n_folds: 5,
            criterion: BasisCriterion::Gcv,
        }
    }
}

/// Select the optimal number of basis functions using a config struct.
///
/// This is the config-based alternative to [`basis_nbasis_cv`]. It takes data
/// parameters directly and reads all tuning parameters from the config.
///
/// The `nbasis_range` tuple `(lo, hi)` is expanded to `lo..=hi` to form the
/// candidate set.
///
/// # Arguments
/// * `data` — Functional data matrix (n × m)
/// * `argvals` — Evaluation points (length m)
/// * `config` — Tuning parameters
///
/// # Errors
///
/// Returns [`crate::FdarError::ComputationFailed`] if no valid result is found
/// for any nbasis in the search range.
#[must_use = "expensive computation whose result should not be discarded"]
pub fn basis_nbasis_cv_with_config(
    data: &FdMatrix,
    argvals: &[f64],
    config: &BasisNbasisCvConfig,
) -> Result<BasisNbasisCvResult, crate::FdarError> {
    let nbasis_range: Vec<usize> = (config.nbasis_range.0..=config.nbasis_range.1).collect();
    basis_nbasis_cv(
        data,
        argvals,
        &nbasis_range,
        &config.basis_type,
        config.criterion,
        config.n_folds,
        config.lambda,
    )
    .ok_or_else(|| crate::FdarError::ComputationFailed {
        operation: "basis_nbasis_cv_with_config",
        detail: "no valid result found in nbasis CV search".to_string(),
    })
}

// ─── Internal Helpers ───────────────────────────────────────────────────────

/// Differentiate column-major basis matrix `lfd_order` times using gradient_uniform.
pub(crate) fn differentiate_basis_columns(
    basis: &[f64],
    n_quad: usize,
    nbasis: usize,
    h: f64,
    lfd_order: usize,
) -> Vec<f64> {
    let mut deriv = basis.to_vec();
    for _ in 0..lfd_order {
        let mut new_deriv = vec![0.0; n_quad * nbasis];
        for j in 0..nbasis {
            let col: Vec<f64> = (0..n_quad).map(|i| deriv[i + j * n_quad]).collect();
            let grad = crate::helpers::gradient_uniform(&col, h);
            for i in 0..n_quad {
                new_deriv[i + j * n_quad] = grad[i];
            }
        }
        deriv = new_deriv;
    }
    deriv
}

/// Integrate symmetric penalty: R[j,k] = ∫ D^m B_j · D^m B_k dt.
pub(crate) fn integrate_symmetric_penalty(
    deriv_basis: &[f64],
    weights: &[f64],
    k: usize,
    n_quad: usize,
) -> Vec<f64> {
    let mut penalty = vec![0.0; k * k];
    for j in 0..k {
        for l in j..k {
            let mut val = 0.0;
            for i in 0..n_quad {
                val += deriv_basis[i + j * n_quad] * deriv_basis[i + l * n_quad] * weights[i];
            }
            penalty[j + l * k] = val;
            penalty[l + j * k] = val;
        }
    }
    penalty
}

/// Evaluate basis functions on argvals, returning (flat column-major, actual_nbasis).
fn evaluate_basis(argvals: &[f64], basis_type: &BasisType, nbasis: usize) -> (Vec<f64>, usize) {
    let m = argvals.len();
    match basis_type {
        BasisType::Bspline { order } => {
            let nknots = nbasis.saturating_sub(*order).max(2);
            let basis = bspline_basis(argvals, nknots, *order);
            let actual = basis.len() / m;
            (basis, actual)
        }
        BasisType::Fourier { period } => {
            let basis = fourier_basis_with_period(argvals, nbasis, *period);
            (basis, nbasis)
        }
    }
}

/// Invert the penalized system matrix via Cholesky or SVD pseudoinverse.
fn invert_penalized_system(system: &DMatrix<f64>, k: usize) -> Option<DMatrix<f64>> {
    if let Some(chol) = system.clone().cholesky() {
        return Some(chol.inverse());
    }
    // SVD fallback
    let svd = nalgebra::SVD::new(system.clone(), true, true);
    let u = svd.u.as_ref()?;
    let v_t = svd.v_t.as_ref()?;
    let max_sv: f64 = svd.singular_values.iter().copied().fold(0.0_f64, f64::max);
    let eps = 1e-10 * max_sv;
    let mut inv = DMatrix::<f64>::zeros(k, k);
    for ii in 0..k {
        for jj in 0..k {
            let mut sum = 0.0;
            for s in 0..k.min(svd.singular_values.len()) {
                if svd.singular_values[s] > eps {
                    sum += v_t[(s, ii)] / svd.singular_values[s] * u[(jj, s)];
                }
            }
            inv[(ii, jj)] = sum;
        }
    }
    Some(inv)
}

/// Project all curves onto basis, returning (coefficients, fitted, total_rss).
fn project_all_curves(
    data: &FdMatrix,
    b_mat: &DMatrix<f64>,
    proj: &DMatrix<f64>,
    n: usize,
    m: usize,
    k: usize,
) -> (FdMatrix, FdMatrix, f64) {
    let mut all_coefs = FdMatrix::zeros(n, k);
    let mut all_fitted = FdMatrix::zeros(n, m);
    let mut total_rss = 0.0;

    for i in 0..n {
        let curve: Vec<f64> = (0..m).map(|j| data[(i, j)]).collect();
        let y_vec = nalgebra::DVector::from_vec(curve.clone());
        let coefs = proj * &y_vec;

        for j in 0..k {
            all_coefs[(i, j)] = coefs[j];
        }
        let fitted = b_mat * &coefs;
        for j in 0..m {
            all_fitted[(i, j)] = fitted[j];
            let resid = curve[j] - fitted[j];
            total_rss += resid * resid;
        }
    }

    (all_coefs, all_fitted, total_rss)
}

/// Compute GCV score.
fn compute_gcv(rss: f64, n_points: f64, edf: f64, m: usize) -> f64 {
    let gcv_denom = 1.0 - edf / m as f64;
    if gcv_denom.abs() > 1e-10 {
        (rss / n_points) / (gcv_denom * gcv_denom)
    } else {
        f64::INFINITY
    }
}

// ─── Nbasis Selection via CV ────────────────────────────────────────────────

/// Criterion for nbasis selection.
#[non_exhaustive]
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum BasisCriterion {
    /// Generalized cross-validation.
    Gcv,
    /// Leave-one-out cross-validation (k-fold).
    Cv,
    /// Akaike Information Criterion.
    Aic,
    /// Bayesian Information Criterion.
    Bic,
}

/// Result of nbasis selection.
#[derive(Debug, Clone, PartialEq)]
#[non_exhaustive]
pub struct BasisNbasisCvResult {
    /// Optimal number of basis functions.
    pub optimal_nbasis: usize,
    /// Score for each nbasis tested.
    pub scores: Vec<f64>,
    /// Range of nbasis values tested.
    pub nbasis_range: Vec<usize>,
    /// Criterion used.
    pub criterion: BasisCriterion,
}

/// Evaluate information criterion (GCV/AIC/BIC) for a range of nbasis values.
fn evaluate_nbasis_info_criterion(
    data: &FdMatrix,
    argvals: &[f64],
    nbasis_range: &[usize],
    basis_type: &BasisType,
    criterion: BasisCriterion,
    lambda: f64,
) -> Vec<f64> {
    let mut scores = Vec::with_capacity(nbasis_range.len());
    for &nb in nbasis_range {
        if nb < 2 {
            scores.push(f64::INFINITY);
            continue;
        }
        let penalty = match basis_type {
            BasisType::Bspline { order } => bspline_penalty_matrix(argvals, nb, *order, 2),
            BasisType::Fourier { period } => fourier_penalty_matrix(nb, *period, 2),
        };
        let fdpar = FdPar {
            basis_type: basis_type.clone(),
            nbasis: nb,
            lambda,
            lfd_order: 2,
            penalty_matrix: penalty,
        };
        match smooth_basis(data, argvals, &fdpar) {
            Ok(result) => {
                let score = match criterion {
                    BasisCriterion::Gcv => result.gcv,
                    BasisCriterion::Aic => result.aic,
                    BasisCriterion::Bic => result.bic,
                    BasisCriterion::Cv => unreachable!(),
                };
                scores.push(score);
            }
            Err(_) => scores.push(f64::INFINITY),
        }
    }
    scores
}

/// Evaluate nbasis via k-fold cross-validation of reconstruction error.
fn evaluate_nbasis_cv(
    data: &FdMatrix,
    argvals: &[f64],
    nbasis_range: &[usize],
    basis_type: &BasisType,
    lambda: f64,
    n_folds: usize,
) -> Vec<f64> {
    let (n, m) = data.shape();
    // Cross-validate over TIME POINTS, not curves. Leaving out curves is
    // ill-posed for per-curve basis fitting — each curve has its own
    // coefficients, so a held-out curve can only be scored against its own data,
    // which is an in-sample residual that decreases monotonically in `nbasis`
    // and always selects the largest candidate (GH #33). Holding out points and
    // predicting them from a fit on the remaining points gives a genuine
    // predictive score that penalizes overfitting and shows an interior minimum.
    let n_folds = n_folds.max(2).min(m);
    let point_folds = crate::cv::create_folds(m, n_folds, 42);
    let mut scores = Vec::with_capacity(nbasis_range.len());

    for &nb in nbasis_range {
        if nb < 2 {
            scores.push(f64::INFINITY);
            continue;
        }
        let penalty = match basis_type {
            BasisType::Bspline { order } => bspline_penalty_matrix(argvals, nb, *order, 2),
            BasisType::Fourier { period } => fourier_penalty_matrix(nb, *period, 2),
        };
        let (basis_flat, actual_k) = evaluate_basis(argvals, basis_type, nb);
        let b_full = DMatrix::from_column_slice(m, actual_k, &basis_flat);
        let r_mat = DMatrix::from_column_slice(actual_k, actual_k, &penalty);

        let mut total_se = 0.0;
        let mut count = 0usize;

        for fold in 0..n_folds {
            let (train_pts, test_pts) = crate::cv::fold_indices(&point_folds, fold);
            if train_pts.is_empty() || test_pts.is_empty() {
                continue;
            }
            // Basis rows at the train/test points, and the projection operator
            // built from the training points only (independent of any curve).
            let b_train = b_full.select_rows(train_pts.iter());
            let b_test = b_full.select_rows(test_pts.iter());
            let btb = b_train.transpose() * &b_train;
            let ridge_eps = 1e-10;
            let system: DMatrix<f64> =
                &btb + lambda * &r_mat + ridge_eps * DMatrix::<f64>::identity(actual_k, actual_k);
            let Some(system_inv) = invert_penalized_system(&system, actual_k) else {
                continue;
            };
            let proj = &system_inv * b_train.transpose(); // (k x |train|)

            for i in 0..n {
                let y_train = nalgebra::DVector::from_iterator(
                    train_pts.len(),
                    train_pts.iter().map(|&j| data[(i, j)]),
                );
                let coefs = &proj * &y_train;
                let pred = &b_test * &coefs; // predictions at held-out points
                for (t_idx, &j) in test_pts.iter().enumerate() {
                    let err = data[(i, j)] - pred[t_idx];
                    total_se += err * err;
                    count += 1;
                }
            }
        }

        if count > 0 {
            scores.push(total_se / count as f64);
        } else {
            scores.push(f64::INFINITY);
        }
    }
    scores
}

/// Select the optimal number of basis functions using multiple criteria
/// (R's `fdata2basis_cv`).
pub fn basis_nbasis_cv(
    data: &FdMatrix,
    argvals: &[f64],
    nbasis_range: &[usize],
    basis_type: &BasisType,
    criterion: BasisCriterion,
    n_folds: usize,
    lambda: f64,
) -> Option<BasisNbasisCvResult> {
    let (n, m) = data.shape();
    if n == 0 || m == 0 || argvals.len() != m || nbasis_range.is_empty() {
        return None;
    }

    let scores = match criterion {
        BasisCriterion::Gcv | BasisCriterion::Aic | BasisCriterion::Bic => {
            evaluate_nbasis_info_criterion(
                data,
                argvals,
                nbasis_range,
                basis_type,
                criterion,
                lambda,
            )
        }
        BasisCriterion::Cv => {
            evaluate_nbasis_cv(data, argvals, nbasis_range, basis_type, lambda, n_folds)
        }
    };

    let (best_idx, _) = scores
        .iter()
        .enumerate()
        .min_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal))?;

    Some(BasisNbasisCvResult {
        optimal_nbasis: nbasis_range[best_idx],
        scores,
        nbasis_range: nbasis_range.to_vec(),
        criterion,
    })
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::test_helpers::uniform_grid;
    use std::f64::consts::PI;

    #[test]
    fn test_bspline_penalty_matrix_symmetric() {
        let t = uniform_grid(101);
        let penalty = bspline_penalty_matrix(&t, 15, 4, 2);
        let _k = 15; // may differ from actual due to knot construction
        let actual_k = (penalty.len() as f64).sqrt() as usize;
        for i in 0..actual_k {
            for j in 0..actual_k {
                assert!(
                    (penalty[i + j * actual_k] - penalty[j + i * actual_k]).abs() < 1e-10,
                    "Penalty matrix not symmetric at ({}, {})",
                    i,
                    j
                );
            }
        }
    }

    #[test]
    fn test_bspline_penalty_matrix_positive_semidefinite() {
        let t = uniform_grid(101);
        let penalty = bspline_penalty_matrix(&t, 10, 4, 2);
        let k = (penalty.len() as f64).sqrt() as usize;
        // Diagonal elements should be non-negative
        for i in 0..k {
            assert!(
                penalty[i + i * k] >= -1e-10,
                "Diagonal element {} is negative: {}",
                i,
                penalty[i + i * k]
            );
        }
    }

    #[test]
    fn test_fourier_penalty_diagonal() {
        let penalty = fourier_penalty_matrix(7, 1.0, 2);
        // Should be diagonal
        for i in 0..7 {
            for j in 0..7 {
                if i != j {
                    assert!(
                        penalty[i + j * 7].abs() < 1e-10,
                        "Off-diagonal ({},{}) = {}",
                        i,
                        j,
                        penalty[i + j * 7]
                    );
                }
            }
        }
        // Constant term should have zero penalty
        assert!(penalty[0].abs() < 1e-10);
        // Higher frequency terms should have larger penalties
        assert!(penalty[1 + 7] > 0.0);
        assert!(penalty[3 + 3 * 7] > penalty[1 + 7]);
    }

    #[test]
    fn test_smooth_basis_bspline() {
        let m = 101;
        let n = 5;
        let t = uniform_grid(m);

        // Generate noisy sine curves
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin() + 0.1 * (i as f64 * 0.3 + j as f64 * 0.01);
            }
        }

        let nbasis = 15;
        let penalty = bspline_penalty_matrix(&t, nbasis, 4, 2);
        let _actual_k = (penalty.len() as f64).sqrt() as usize;

        let fdpar = FdPar {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis,
            lambda: 1e-4,
            lfd_order: 2,
            penalty_matrix: penalty,
        };

        let result = smooth_basis(&data, &t, &fdpar);
        assert!(result.is_ok(), "smooth_basis should succeed");

        let res = result.unwrap();
        assert_eq!(res.fitted.shape(), (n, m));
        assert_eq!(res.coefficients.nrows(), n);
        assert!(res.edf > 0.0, "EDF should be positive");
        assert!(res.gcv > 0.0, "GCV should be positive");
    }

    #[test]
    fn test_smooth_basis_fourier() {
        let m = 101;
        let n = 3;
        let t = uniform_grid(m);

        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin() + (4.0 * PI * t[j]).cos();
            }
        }

        let nbasis = 7;
        let period = 1.0;
        let penalty = fourier_penalty_matrix(nbasis, period, 2);

        let fdpar = FdPar {
            basis_type: BasisType::Fourier { period },
            nbasis,
            lambda: 1e-6,
            lfd_order: 2,
            penalty_matrix: penalty,
        };

        let result = smooth_basis(&data, &t, &fdpar);
        assert!(result.is_ok());

        let res = result.unwrap();
        // Fourier basis should fit periodic data well
        for j in 0..m {
            let expected = (2.0 * PI * t[j]).sin() + (4.0 * PI * t[j]).cos();
            assert!(
                (res.fitted[(0, j)] - expected).abs() < 0.1,
                "Fourier fit poor at j={}: got {}, expected {}",
                j,
                res.fitted[(0, j)],
                expected
            );
        }
    }

    #[test]
    fn test_smooth_basis_gcv_selects_reasonable_lambda() {
        let m = 101;
        let n = 5;
        let t = uniform_grid(m);

        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] =
                    (2.0 * PI * t[j]).sin() + 0.1 * ((i * 37 + j * 13) % 20) as f64 / 20.0;
            }
        }

        let basis_type = BasisType::Bspline { order: 4 };
        let result = smooth_basis_gcv(&data, &t, &basis_type, 15, 2, (-8.0, 4.0), 25);
        assert!(result.is_some(), "GCV search should succeed");
    }

    #[test]
    fn test_smooth_basis_aic_matches_brute_force_grid() {
        // smooth_basis_aic must return the SmoothBasisResult whose aic equals the
        // argmin of an explicit AIC grid built by calling smooth_basis at each
        // grid lambda and reading result.aic over the identical log-lambda grid.
        let m = 81;
        let n = 4;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] =
                    (2.0 * PI * t[j]).sin() + 0.15 * ((i * 41 + j * 7) % 23) as f64 / 23.0;
            }
        }

        let basis_type = BasisType::Bspline { order: 4 };
        let nbasis = 12;
        let lfd_order = 2;
        let range = (-8.0, 4.0);
        let n_grid = 25;

        // Brute-force the AIC grid using the same penalty + grid as the selector.
        let penalty = bspline_penalty_matrix(&t, nbasis, 4, lfd_order);
        let (lo, hi) = range;
        let mut brute_best_aic = f64::INFINITY;
        for k in 0..n_grid {
            let log_lam = lo + (hi - lo) * k as f64 / (n_grid - 1) as f64;
            let lam = 10.0_f64.powf(log_lam);
            let fdpar = FdPar {
                basis_type: basis_type.clone(),
                nbasis,
                lambda: lam,
                lfd_order,
                penalty_matrix: penalty.clone(),
            };
            if let Ok(result) = smooth_basis(&data, &t, &fdpar) {
                if result.aic < brute_best_aic {
                    brute_best_aic = result.aic;
                }
            }
        }

        let selected =
            smooth_basis_aic(&data, &t, &basis_type, nbasis, lfd_order, range, n_grid).unwrap();
        assert!(
            (selected.aic - brute_best_aic).abs() < 1e-9,
            "selected aic={}, brute-force min aic={}",
            selected.aic,
            brute_best_aic
        );
    }

    #[test]
    fn test_smooth_basis_aic_prefers_smoother_fit_than_smallest_lambda() {
        // On noisy data, the AIC-selected fit must have lower edf (be smoother)
        // than the smallest-lambda candidate from the same grid, which over-fits.
        let m = 81;
        let n = 4;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(n, m);
        // Smooth trend + deterministic pseudo-noise.
        for i in 0..n {
            for j in 0..m {
                let noise = ((((i * 91 + j * 53) % 101) as f64) / 101.0) - 0.5;
                data[(i, j)] = (2.0 * PI * t[j]).sin() + 0.6 * noise;
            }
        }

        let basis_type = BasisType::Bspline { order: 4 };
        let nbasis = 15;
        let lfd_order = 2;
        let range = (-8.0, 4.0);
        let n_grid = 25;

        // Smallest-lambda candidate (most flexible / highest edf).
        let penalty = bspline_penalty_matrix(&t, nbasis, 4, lfd_order);
        let smallest_lam = 10.0_f64.powf(range.0);
        let fdpar_small = FdPar {
            basis_type: basis_type.clone(),
            nbasis,
            lambda: smallest_lam,
            lfd_order,
            penalty_matrix: penalty.clone(),
        };
        let overfit = smooth_basis(&data, &t, &fdpar_small).unwrap();

        let selected =
            smooth_basis_aic(&data, &t, &basis_type, nbasis, lfd_order, range, n_grid).unwrap();

        assert!(
            selected.edf < overfit.edf,
            "AIC-selected edf ({}) should be smaller (smoother) than the smallest-lambda edf ({})",
            selected.edf,
            overfit.edf
        );
    }

    #[test]
    fn test_smooth_basis_large_lambda_reduces_edf() {
        let m = 101;
        let n = 3;
        let t = uniform_grid(m);

        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin();
            }
        }

        let nbasis = 15;
        let penalty = bspline_penalty_matrix(&t, nbasis, 4, 2);
        let _actual_k = (penalty.len() as f64).sqrt() as usize;

        let fdpar_small = FdPar {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis,
            lambda: 1e-8,
            lfd_order: 2,
            penalty_matrix: penalty.clone(),
        };
        let fdpar_large = FdPar {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis,
            lambda: 1e2,
            lfd_order: 2,
            penalty_matrix: penalty,
        };

        let res_small = smooth_basis(&data, &t, &fdpar_small).unwrap();
        let res_large = smooth_basis(&data, &t, &fdpar_large).unwrap();

        assert!(
            res_large.edf < res_small.edf,
            "Larger lambda should reduce EDF: {} vs {}",
            res_large.edf,
            res_small.edf
        );
    }

    // ============== basis_nbasis_cv tests ==============

    #[test]
    fn test_basis_nbasis_cv_gcv() {
        let m = 101;
        let n = 5;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] =
                    (2.0 * PI * t[j]).sin() + 0.1 * ((i * 37 + j * 13) % 20) as f64 / 20.0;
            }
        }

        let nbasis_range: Vec<usize> = (4..=20).step_by(2).collect();
        let result = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Gcv,
            5,
            1e-4,
        );
        assert!(result.is_some());
        let res = result.unwrap();
        assert!(nbasis_range.contains(&res.optimal_nbasis));
        assert_eq!(res.scores.len(), nbasis_range.len());
        assert_eq!(res.criterion, BasisCriterion::Gcv);
    }

    #[test]
    fn test_basis_nbasis_cv_aic_bic() {
        let m = 51;
        let n = 5;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin();
            }
        }

        let nbasis_range: Vec<usize> = vec![5, 7, 9, 11];
        let aic_result = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Aic,
            5,
            0.0,
        );
        let bic_result = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Bic,
            5,
            0.0,
        );
        assert!(aic_result.is_some());
        assert!(bic_result.is_some());
    }

    #[test]
    fn test_basis_nbasis_cv_kfold() {
        let m = 51;
        let n = 10;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin() + 0.05 * ((i * 7 + j * 3) % 10) as f64;
            }
        }

        let nbasis_range: Vec<usize> = vec![5, 7, 9];
        let result = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Cv,
            5,
            1e-4,
        );
        assert!(result.is_some());
        let res = result.unwrap();
        assert!(nbasis_range.contains(&res.optimal_nbasis));
        assert_eq!(res.criterion, BasisCriterion::Cv);
    }

    /// Regression for GH #33: the CV path scored held-out curves against their
    /// own data (no true hold-out), so scores fell monotonically in `n_basis`
    /// and it always selected the maximum candidate. With point-wise CV,
    /// overfitting the noise is penalized and the maximum is not chosen.
    #[test]
    fn test_basis_nbasis_cv_penalizes_overfitting() {
        let m = 120;
        let n = 6;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                // Smooth signal + deterministic pseudo-noise.
                let noise = 0.2 * (((i * 31 + j * 17) % 13) as f64 / 13.0 - 0.5);
                data[(i, j)] = (2.0 * PI * t[j]).sin() + noise;
            }
        }

        let nbasis_range: Vec<usize> = vec![5, 8, 12, 20, 30];
        let res = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Cv,
            5,
            1e-6,
        )
        .unwrap();

        assert_ne!(
            res.optimal_nbasis, 30,
            "CV must not always select the maximum n_basis (GH #33); scores={:?}",
            res.scores
        );
        let monotone_decreasing = res.scores.windows(2).all(|w| w[1] <= w[0] + 1e-12);
        assert!(
            !monotone_decreasing,
            "CV scores must not be monotone-decreasing in n_basis; scores={:?}",
            res.scores
        );
    }

    // ============== Comprehensive additional tests ==============

    // Helper: generate standard test data (sine + high-freq component)
    fn make_test_data(n: usize, m: usize) -> (FdMatrix, Vec<f64>) {
        let t: Vec<f64> = (0..m).map(|i| i as f64 / (m - 1) as f64).collect();
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin()
                    + 0.1 * (10.0 * t[j]).sin()
                    + 0.05 * ((i * 37 + j * 13) % 20) as f64 / 20.0;
            }
        }
        (data, t)
    }

    // Helper: create an FdPar for B-spline smoothing
    fn make_bspline_fdpar(argvals: &[f64], nbasis: usize, lambda: f64) -> FdPar {
        let penalty = bspline_penalty_matrix(argvals, nbasis, 4, 2);
        FdPar {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis,
            lambda,
            lfd_order: 2,
            penalty_matrix: penalty,
        }
    }

    // Helper: create an FdPar for Fourier smoothing
    fn make_fourier_fdpar(nbasis: usize, period: f64, lambda: f64) -> FdPar {
        let penalty = fourier_penalty_matrix(nbasis, period, 2);
        FdPar {
            basis_type: BasisType::Fourier { period },
            nbasis,
            lambda,
            lfd_order: 2,
            penalty_matrix: penalty,
        }
    }

    // ─── BasisType enum tests ───────────────────────────────────────────────

    #[test]
    fn test_basis_type_bspline_variant() {
        let bt = BasisType::Bspline { order: 4 };
        assert_eq!(bt, BasisType::Bspline { order: 4 });
        // Different orders are not equal
        assert_ne!(bt, BasisType::Bspline { order: 3 });
    }

    #[test]
    fn test_basis_type_fourier_variant() {
        let bt = BasisType::Fourier { period: 1.0 };
        assert_eq!(bt, BasisType::Fourier { period: 1.0 });
        assert_ne!(bt, BasisType::Fourier { period: 2.0 });
    }

    #[test]
    fn test_basis_type_cross_variant_inequality() {
        let bspline = BasisType::Bspline { order: 4 };
        let fourier = BasisType::Fourier { period: 1.0 };
        assert_ne!(bspline, fourier);
    }

    #[test]
    fn test_basis_type_clone_and_debug() {
        let bt = BasisType::Bspline { order: 4 };
        let cloned = bt.clone();
        assert_eq!(bt, cloned);
        let debug_str = format!("{:?}", bt);
        assert!(debug_str.contains("Bspline"));
        assert!(debug_str.contains("4"));
    }

    // ─── FdPar struct tests ─────────────────────────────────────────────────

    #[test]
    fn test_fdpar_construction_and_fields() {
        let penalty = vec![1.0, 0.0, 0.0, 1.0];
        let fdpar = FdPar {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis: 2,
            lambda: 0.01,
            lfd_order: 2,
            penalty_matrix: penalty.clone(),
        };
        assert_eq!(fdpar.nbasis, 2);
        assert!((fdpar.lambda - 0.01).abs() < 1e-15);
        assert_eq!(fdpar.lfd_order, 2);
        assert_eq!(fdpar.penalty_matrix.len(), 4);
    }

    #[test]
    fn test_fdpar_clone_and_debug() {
        let t = uniform_grid(50);
        let fdpar = make_bspline_fdpar(&t, 8, 1e-3);
        let cloned = fdpar.clone();
        assert_eq!(fdpar, cloned);
        let debug_str = format!("{:?}", fdpar);
        assert!(debug_str.contains("FdPar"));
    }

    // ─── BasisCriterion enum tests ──────────────────────────────────────────

    #[test]
    fn test_basis_criterion_variants() {
        assert_eq!(BasisCriterion::Gcv, BasisCriterion::Gcv);
        assert_eq!(BasisCriterion::Cv, BasisCriterion::Cv);
        assert_eq!(BasisCriterion::Aic, BasisCriterion::Aic);
        assert_eq!(BasisCriterion::Bic, BasisCriterion::Bic);
        assert_ne!(BasisCriterion::Gcv, BasisCriterion::Aic);
        assert_ne!(BasisCriterion::Cv, BasisCriterion::Bic);
    }

    #[test]
    fn test_basis_criterion_copy() {
        let c = BasisCriterion::Gcv;
        let copied = c; // Copy
        assert_eq!(c, copied);
    }

    #[test]
    fn test_basis_criterion_debug() {
        let debug_str = format!("{:?}", BasisCriterion::Bic);
        assert!(debug_str.contains("Bic"));
    }

    // ─── SmoothBasisResult tests ────────────────────────────────────────────

    #[test]
    fn test_smooth_basis_result_all_fields() {
        let (data, t) = make_test_data(3, 50);
        let fdpar = make_bspline_fdpar(&t, 10, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();

        // coefficients: n curves x k basis functions
        assert_eq!(res.coefficients.nrows(), 3);
        assert!(res.coefficients.ncols() > 0);
        assert_eq!(res.nbasis, res.coefficients.ncols());
        // fitted: n x m
        assert_eq!(res.fitted.shape(), (3, 50));
        // edf should be between 1 and nbasis
        assert!(res.edf > 0.0 && res.edf <= res.nbasis as f64);
        // gcv, aic, bic should be finite
        assert!(res.gcv.is_finite());
        assert!(res.aic.is_finite());
        assert!(res.bic.is_finite());
        // penalty_matrix should be k x k
        let k = res.nbasis;
        assert_eq!(res.penalty_matrix.len(), k * k);
    }

    #[test]
    fn test_smooth_basis_result_clone() {
        let (data, t) = make_test_data(2, 50);
        let fdpar = make_bspline_fdpar(&t, 8, 1e-3);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        let cloned = res.clone();
        assert_eq!(res, cloned);
    }

    // ─── smooth_basis: B-spline detailed tests ──────────────────────────────

    #[test]
    fn test_smooth_basis_bspline_coefficient_shape() {
        let (data, t) = make_test_data(4, 50);
        let nbasis = 12;
        let fdpar = make_bspline_fdpar(&t, nbasis, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        assert_eq!(res.coefficients.nrows(), 4);
        // actual nbasis may differ from requested due to knot construction
        assert!(res.coefficients.ncols() >= 2);
        assert_eq!(res.nbasis, res.coefficients.ncols());
    }

    #[test]
    fn test_smooth_basis_bspline_fitted_values_shape() {
        let m = 80;
        let n = 6;
        let (data, t) = make_test_data(n, m);
        let fdpar = make_bspline_fdpar(&t, 15, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        assert_eq!(res.fitted.shape(), (n, m));
    }

    #[test]
    fn test_smooth_basis_bspline_zero_lambda_interpolates() {
        // With lambda=0, the smoother should nearly interpolate the data
        let m = 30;
        let n = 2;
        let (data, t) = make_test_data(n, m);
        let fdpar = make_bspline_fdpar(&t, 15, 0.0);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();

        // Residuals should be very small (near interpolation)
        let mut max_resid = 0.0_f64;
        for i in 0..n {
            for j in 0..m {
                let resid = (data[(i, j)] - res.fitted[(i, j)]).abs();
                max_resid = max_resid.max(resid);
            }
        }
        assert!(
            max_resid < 0.5,
            "Zero-lambda B-spline should closely interpolate; max_resid = {}",
            max_resid
        );
    }

    #[test]
    fn test_smooth_basis_bspline_large_lambda_oversmooths() {
        // With very large lambda, the fit should be much smoother (lower variance)
        // than with small lambda
        let m = 50;
        let n = 1;
        let (data, t) = make_test_data(n, m);

        let fdpar_small = make_bspline_fdpar(&t, 15, 1e-6);
        let res_small = smooth_basis(&data, &t, &fdpar_small).unwrap();

        let fdpar_large = make_bspline_fdpar(&t, 15, 1e6);
        let res_large = smooth_basis(&data, &t, &fdpar_large).unwrap();

        let compute_variance = |fitted: &FdMatrix, row: usize, ncols: usize| -> f64 {
            let vals: Vec<f64> = (0..ncols).map(|j| fitted[(row, j)]).collect();
            let mean = vals.iter().sum::<f64>() / ncols as f64;
            vals.iter().map(|&v| (v - mean).powi(2)).sum::<f64>() / ncols as f64
        };

        let var_small = compute_variance(&res_small.fitted, 0, m);
        let var_large = compute_variance(&res_large.fitted, 0, m);
        assert!(
            var_large < var_small,
            "Large lambda should yield lower variance fit: var_large={}, var_small={}",
            var_large,
            var_small
        );
    }

    #[test]
    fn test_smooth_basis_bspline_penalty_effect_on_smoothness() {
        // Compare roughness of fits with small vs large lambda
        let m = 50;
        let n = 1;
        let (data, t) = make_test_data(n, m);

        let fdpar_small = make_bspline_fdpar(&t, 15, 1e-8);
        let fdpar_large = make_bspline_fdpar(&t, 15, 1.0);

        let res_small = smooth_basis(&data, &t, &fdpar_small).unwrap();
        let res_large = smooth_basis(&data, &t, &fdpar_large).unwrap();

        // Measure roughness as sum of squared second differences
        let roughness = |fitted: &FdMatrix, row: usize, ncols: usize| -> f64 {
            (1..ncols - 1)
                .map(|j| {
                    let d2 = fitted[(row, j + 1)] - 2.0 * fitted[(row, j)] + fitted[(row, j - 1)];
                    d2 * d2
                })
                .sum::<f64>()
        };

        let r_small = roughness(&res_small.fitted, 0, m);
        let r_large = roughness(&res_large.fitted, 0, m);
        assert!(
            r_large < r_small,
            "Larger lambda should produce smoother fit: roughness_large={}, roughness_small={}",
            r_large,
            r_small
        );
    }

    #[test]
    fn test_smooth_basis_bspline_single_curve() {
        let m = 50;
        let (data, t) = make_test_data(1, m);
        let fdpar = make_bspline_fdpar(&t, 10, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        assert_eq!(res.fitted.nrows(), 1);
        assert_eq!(res.fitted.ncols(), m);
        assert!(res.gcv.is_finite());
    }

    #[test]
    fn test_smooth_basis_bspline_many_curves() {
        let m = 50;
        let n = 20;
        let (data, t) = make_test_data(n, m);
        let fdpar = make_bspline_fdpar(&t, 10, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        assert_eq!(res.fitted.nrows(), n);
        assert_eq!(res.coefficients.nrows(), n);
    }

    #[test]
    fn test_smooth_basis_bspline_minimal_nbasis() {
        // nbasis = 2 is the minimum allowed
        let m = 50;
        let (data, t) = make_test_data(1, m);
        let fdpar = make_bspline_fdpar(&t, 2, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar);
        // Should succeed (or at least not panic); the fit may be poor
        assert!(res.is_ok());
    }

    #[test]
    fn test_smooth_basis_bspline_different_orders() {
        let m = 50;
        let (data, t) = make_test_data(2, m);
        // Order 3 (quadratic B-splines)
        let penalty3 = bspline_penalty_matrix(&t, 10, 3, 2);
        let fdpar3 = FdPar {
            basis_type: BasisType::Bspline { order: 3 },
            nbasis: 10,
            lambda: 1e-4,
            lfd_order: 2,
            penalty_matrix: penalty3,
        };
        let res3 = smooth_basis(&data, &t, &fdpar3);
        assert!(res3.is_ok());

        // Order 5 (quartic B-splines)
        let penalty5 = bspline_penalty_matrix(&t, 10, 5, 2);
        let fdpar5 = FdPar {
            basis_type: BasisType::Bspline { order: 5 },
            nbasis: 10,
            lambda: 1e-4,
            lfd_order: 2,
            penalty_matrix: penalty5,
        };
        let res5 = smooth_basis(&data, &t, &fdpar5);
        assert!(res5.is_ok());
    }

    // ─── smooth_basis: Fourier detailed tests ───────────────────────────────

    #[test]
    fn test_smooth_basis_fourier_coefficient_shape() {
        let m = 50;
        let n = 3;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin();
            }
        }
        let nbasis = 7;
        let fdpar = make_fourier_fdpar(nbasis, 1.0, 1e-6);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        assert_eq!(res.coefficients.nrows(), n);
        assert_eq!(res.coefficients.ncols(), nbasis);
        assert_eq!(res.nbasis, nbasis);
    }

    #[test]
    fn test_smooth_basis_fourier_fits_pure_sine() {
        // Fourier basis should perfectly fit a pure sine with enough basis fns
        let m = 100;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(1, m);
        for j in 0..m {
            data[(0, j)] = (2.0 * PI * t[j]).sin();
        }
        let fdpar = make_fourier_fdpar(5, 1.0, 1e-8);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();

        for j in 0..m {
            let expected = (2.0 * PI * t[j]).sin();
            assert!(
                (res.fitted[(0, j)] - expected).abs() < 0.05,
                "Fourier should fit pure sine; j={}, got={}, expected={}",
                j,
                res.fitted[(0, j)],
                expected
            );
        }
    }

    #[test]
    fn test_smooth_basis_fourier_different_periods() {
        let m = 50;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(1, m);
        for j in 0..m {
            data[(0, j)] = (2.0 * PI * t[j]).sin();
        }

        // Period = 1.0 (matches the data)
        let fdpar1 = make_fourier_fdpar(7, 1.0, 1e-6);
        let res1 = smooth_basis(&data, &t, &fdpar1).unwrap();

        // Period = 2.0 (mismatch, but should still produce a result)
        let fdpar2 = make_fourier_fdpar(7, 2.0, 1e-6);
        let res2 = smooth_basis(&data, &t, &fdpar2).unwrap();

        // Both should succeed and have valid shapes
        assert_eq!(res1.fitted.shape(), (1, m));
        assert_eq!(res2.fitted.shape(), (1, m));
    }

    #[test]
    fn test_smooth_basis_fourier_zero_lambda() {
        let m = 50;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(1, m);
        for j in 0..m {
            data[(0, j)] = (2.0 * PI * t[j]).sin() + (4.0 * PI * t[j]).cos();
        }
        let fdpar = make_fourier_fdpar(9, 1.0, 0.0);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        assert_eq!(res.fitted.shape(), (1, m));
        // EDF should be close to nbasis with zero penalty
        assert!(res.edf > 1.0);
    }

    #[test]
    fn test_smooth_basis_fourier_large_lambda() {
        let m = 50;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(1, m);
        for j in 0..m {
            data[(0, j)] = (2.0 * PI * t[j]).sin();
        }
        let fdpar = make_fourier_fdpar(9, 1.0, 1e6);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        // EDF should be very small with huge penalty
        assert!(
            res.edf < 5.0,
            "Large lambda should reduce EDF; edf={}",
            res.edf
        );
    }

    // ─── smooth_basis: Lambda comparison tests ──────────────────────────────

    #[test]
    fn test_smooth_basis_lambda_gradient_edf() {
        // EDF should monotonically decrease with increasing lambda
        let m = 50;
        let (data, t) = make_test_data(3, m);
        let lambdas = [1e-8, 1e-4, 1e-2, 1.0, 1e2];
        let mut prev_edf = f64::INFINITY;
        for &lam in &lambdas {
            let fdpar = make_bspline_fdpar(&t, 12, lam);
            let res = smooth_basis(&data, &t, &fdpar).unwrap();
            assert!(
                res.edf <= prev_edf + 0.01,
                "EDF should decrease: lambda={}, edf={}, prev_edf={}",
                lam,
                res.edf,
                prev_edf
            );
            prev_edf = res.edf;
        }
    }

    #[test]
    fn test_smooth_basis_lambda_gradient_rss() {
        // RSS should monotonically increase with increasing lambda
        let m = 50;
        let n = 2;
        let (data, t) = make_test_data(n, m);
        let lambdas = [0.0, 1e-6, 1e-2, 1.0, 1e4];
        let mut prev_rss = -1.0;
        for &lam in &lambdas {
            let fdpar = make_bspline_fdpar(&t, 12, lam);
            let res = smooth_basis(&data, &t, &fdpar).unwrap();
            let mut rss = 0.0;
            for i in 0..n {
                for j in 0..m {
                    rss += (data[(i, j)] - res.fitted[(i, j)]).powi(2);
                }
            }
            assert!(
                rss >= prev_rss - 1e-8,
                "RSS should increase: lambda={}, rss={}, prev_rss={}",
                lam,
                rss,
                prev_rss
            );
            prev_rss = rss;
        }
    }

    // ─── smooth_basis: Error cases ──────────────────────────────────────────

    #[test]
    fn test_smooth_basis_empty_data_rows() {
        let t = uniform_grid(50);
        let data = FdMatrix::zeros(0, 50);
        let fdpar = make_bspline_fdpar(&t, 10, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar);
        assert!(res.is_err());
    }

    #[test]
    fn test_smooth_basis_empty_data_cols() {
        let data = FdMatrix::zeros(5, 0);
        let fdpar = FdPar {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis: 10,
            lambda: 1e-4,
            lfd_order: 2,
            penalty_matrix: vec![0.0; 100],
        };
        let res = smooth_basis(&data, &[], &fdpar);
        assert!(res.is_err());
    }

    #[test]
    fn test_smooth_basis_mismatched_argvals() {
        let t = uniform_grid(50);
        let data = FdMatrix::zeros(3, 40); // m=40 but argvals has 50
        let fdpar = make_bspline_fdpar(&t, 10, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar);
        assert!(res.is_err());
    }

    #[test]
    fn test_smooth_basis_nbasis_too_small() {
        let t = uniform_grid(50);
        let data = FdMatrix::zeros(3, 50);
        // nbasis = 1, which is below minimum of 2
        let fdpar = FdPar {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis: 1,
            lambda: 1e-4,
            lfd_order: 2,
            penalty_matrix: vec![0.0; 1],
        };
        let res = smooth_basis(&data, &t, &fdpar);
        assert!(res.is_err());
    }

    #[test]
    fn test_smooth_basis_error_is_invalid_dimension() {
        let t = uniform_grid(50);
        let data = FdMatrix::zeros(0, 50);
        let fdpar = make_bspline_fdpar(&t, 10, 1e-4);
        let err = smooth_basis(&data, &t, &fdpar).unwrap_err();
        match err {
            crate::FdarError::InvalidDimension { .. } => {} // expected
            other => panic!("Expected InvalidDimension, got {:?}", other),
        }
    }

    // ─── Penalty matrix detailed tests ──────────────────────────────────────

    #[test]
    fn test_bspline_penalty_matrix_different_orders() {
        let t = uniform_grid(101);
        // Order 1 penalty (penalize derivatives)
        let p1 = bspline_penalty_matrix(&t, 10, 4, 1);
        // Order 2 penalty (penalize curvature)
        let p2 = bspline_penalty_matrix(&t, 10, 4, 2);
        // Both should be square and same size
        assert_eq!(p1.len(), p2.len());
        // But they should differ
        let diff: f64 = p1.iter().zip(p2.iter()).map(|(a, b)| (a - b).abs()).sum();
        assert!(
            diff > 1e-10,
            "Different lfd_orders should produce different penalties"
        );
    }

    #[test]
    fn test_bspline_penalty_matrix_edge_cases() {
        // Too few argvals
        let t = vec![0.0];
        let p = bspline_penalty_matrix(&t, 10, 4, 2);
        // Should return zero matrix
        assert!(p.iter().all(|&v| v == 0.0));

        // nbasis < 2
        let t2 = uniform_grid(50);
        let p2 = bspline_penalty_matrix(&t2, 1, 4, 2);
        assert!(p2.iter().all(|&v| v == 0.0));

        // lfd_order >= order
        let p3 = bspline_penalty_matrix(&t2, 10, 4, 4);
        assert!(p3.iter().all(|&v| v == 0.0));
    }

    #[test]
    fn test_bspline_penalty_nonnegative_diagonal() {
        let t = uniform_grid(101);
        for nbasis in [5, 10, 20] {
            let p = bspline_penalty_matrix(&t, nbasis, 4, 2);
            let k = (p.len() as f64).sqrt() as usize;
            for i in 0..k {
                assert!(
                    p[i + i * k] >= -1e-10,
                    "Diagonal ({},{}) negative for nbasis={}: {}",
                    i,
                    i,
                    nbasis,
                    p[i + i * k]
                );
            }
        }
    }

    #[test]
    fn test_fourier_penalty_increasing_with_frequency() {
        let penalty = fourier_penalty_matrix(11, 1.0, 2);
        let k = 11;
        // Constant term is zero
        assert!(penalty[0].abs() < 1e-15);
        // Pairs: (1,2) -> freq 1, (3,4) -> freq 2, etc.
        let mut prev_eigenval = 0.0;
        for freq in 1..=5 {
            let idx_sin = 2 * freq - 1;
            let eigenval = penalty[idx_sin + idx_sin * k];
            assert!(
                eigenval > prev_eigenval,
                "Higher frequency should have larger penalty: freq={}, eigenval={}, prev={}",
                freq,
                eigenval,
                prev_eigenval
            );
            prev_eigenval = eigenval;
            // cos and sin of same frequency should have same penalty
            let idx_cos = 2 * freq;
            if idx_cos < k {
                assert!(
                    (penalty[idx_cos + idx_cos * k] - eigenval).abs() < 1e-10,
                    "Sin and cos penalty should match at freq {}",
                    freq
                );
            }
        }
    }

    #[test]
    fn test_fourier_penalty_different_periods() {
        let p1 = fourier_penalty_matrix(7, 1.0, 2);
        let p2 = fourier_penalty_matrix(7, 2.0, 2);
        // Longer period -> smaller omega -> smaller penalty eigenvalues
        for i in 1..7 {
            assert!(
                p2[i + i * 7] < p1[i + i * 7] || (p1[i + i * 7] == 0.0 && p2[i + i * 7] == 0.0),
                "Longer period should have smaller penalties at i={}",
                i
            );
        }
    }

    #[test]
    fn test_fourier_penalty_first_order() {
        // lfd_order = 1: penalize first derivative
        let p = fourier_penalty_matrix(5, 1.0, 1);
        // Eigenvalues: (2*pi*freq)^2 for lfd_order=1
        let omega1 = 2.0 * PI;
        let expected1 = omega1.powi(2);
        assert!(
            (p[1 + 5] - expected1).abs() < 1e-6,
            "First-order penalty eigenval: got {}, expected {}",
            p[1 + 5],
            expected1
        );
    }

    #[test]
    fn test_fourier_penalty_zero_nbasis() {
        let p = fourier_penalty_matrix(0, 1.0, 2);
        assert!(p.is_empty());
    }

    #[test]
    fn test_fourier_penalty_nbasis_one() {
        let p = fourier_penalty_matrix(1, 1.0, 2);
        assert_eq!(p.len(), 1);
        assert!(p[0].abs() < 1e-15); // constant term has zero penalty
    }

    // ─── smooth_basis_gcv detailed tests ────────────────────────────────────

    #[test]
    fn test_smooth_basis_gcv_returns_valid_result() {
        let (data, t) = make_test_data(5, 50);
        let bt = BasisType::Bspline { order: 4 };
        let result = smooth_basis_gcv(&data, &t, &bt, 12, 2, (-6.0, 2.0), 20);
        assert!(result.is_some());
        let res = result.unwrap();
        assert_eq!(res.fitted.shape(), (5, 50));
        assert!(res.gcv.is_finite());
        assert!(res.edf > 0.0);
    }

    #[test]
    fn test_smooth_basis_gcv_fourier() {
        let m = 80;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(3, m);
        for i in 0..3 {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin() + 0.5 * (4.0 * PI * t[j]).cos();
            }
        }
        let bt = BasisType::Fourier { period: 1.0 };
        let result = smooth_basis_gcv(&data, &t, &bt, 9, 2, (-8.0, 4.0), 25);
        assert!(result.is_some());
        let res = result.unwrap();
        assert_eq!(res.fitted.nrows(), 3);
        assert_eq!(res.nbasis, 9);
    }

    #[test]
    fn test_smooth_basis_gcv_selects_finite_gcv() {
        let (data, t) = make_test_data(5, 60);
        let bt = BasisType::Bspline { order: 4 };
        let res = smooth_basis_gcv(&data, &t, &bt, 12, 2, (-6.0, 2.0), 15).unwrap();
        assert!(res.gcv.is_finite());
        assert!(res.gcv > 0.0);
    }

    #[test]
    fn test_smooth_basis_gcv_empty_data() {
        let data = FdMatrix::zeros(0, 50);
        let t = uniform_grid(50);
        let bt = BasisType::Bspline { order: 4 };
        let result = smooth_basis_gcv(&data, &t, &bt, 10, 2, (-6.0, 2.0), 10);
        // Should return None since smooth_basis will error for empty data
        assert!(result.is_none());
    }

    #[test]
    fn test_smooth_basis_gcv_empty_argvals() {
        let data = FdMatrix::zeros(5, 0);
        let bt = BasisType::Bspline { order: 4 };
        let result = smooth_basis_gcv(&data, &[], &bt, 10, 2, (-6.0, 2.0), 10);
        assert!(result.is_none());
    }

    #[test]
    fn test_smooth_basis_gcv_nbasis_too_small() {
        let (data, t) = make_test_data(5, 50);
        let bt = BasisType::Bspline { order: 4 };
        let result = smooth_basis_gcv(&data, &t, &bt, 1, 2, (-6.0, 2.0), 10);
        assert!(result.is_none());
    }

    #[test]
    fn test_smooth_basis_gcv_ngrid_too_small() {
        let (data, t) = make_test_data(5, 50);
        let bt = BasisType::Bspline { order: 4 };
        let result = smooth_basis_gcv(&data, &t, &bt, 10, 2, (-6.0, 2.0), 1);
        assert!(result.is_none());
    }

    #[test]
    fn test_smooth_basis_gcv_narrow_range() {
        let (data, t) = make_test_data(3, 50);
        let bt = BasisType::Bspline { order: 4 };
        // Very narrow search range
        let result = smooth_basis_gcv(&data, &t, &bt, 10, 2, (-3.0, -2.0), 5);
        assert!(result.is_some());
    }

    #[test]
    fn test_smooth_basis_gcv_wide_range() {
        let (data, t) = make_test_data(3, 50);
        let bt = BasisType::Bspline { order: 4 };
        // Very wide search range
        let result = smooth_basis_gcv(&data, &t, &bt, 10, 2, (-12.0, 8.0), 30);
        assert!(result.is_some());
    }

    // ─── basis_nbasis_cv detailed tests ─────────────────────────────────────

    #[test]
    fn test_basis_nbasis_cv_scores_length() {
        let (data, t) = make_test_data(5, 50);
        let nbasis_range: Vec<usize> = vec![4, 6, 8, 10, 12];
        let res = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Gcv,
            5,
            1e-4,
        )
        .unwrap();
        assert_eq!(res.scores.len(), 5);
        assert_eq!(res.nbasis_range.len(), 5);
        assert_eq!(res.nbasis_range, nbasis_range);
    }

    #[test]
    fn test_basis_nbasis_cv_optimal_within_range() {
        let (data, t) = make_test_data(8, 50);
        let nbasis_range: Vec<usize> = vec![5, 7, 9, 11, 13, 15];
        for criterion in [
            BasisCriterion::Gcv,
            BasisCriterion::Aic,
            BasisCriterion::Bic,
        ] {
            let res = basis_nbasis_cv(
                &data,
                &t,
                &nbasis_range,
                &BasisType::Bspline { order: 4 },
                criterion,
                5,
                1e-4,
            )
            .unwrap();
            assert!(
                nbasis_range.contains(&res.optimal_nbasis),
                "optimal_nbasis {} not in range for {:?}",
                res.optimal_nbasis,
                criterion
            );
        }
    }

    #[test]
    fn test_basis_nbasis_cv_fourier_gcv() {
        let m = 80;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(5, m);
        for i in 0..5 {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin()
                    + 0.3 * (4.0 * PI * t[j]).cos()
                    + 0.02 * ((i * 7 + j * 3) % 10) as f64;
            }
        }
        let nbasis_range: Vec<usize> = vec![5, 7, 9, 11];
        let res = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Fourier { period: 1.0 },
            BasisCriterion::Gcv,
            5,
            1e-4,
        )
        .unwrap();
        assert!(nbasis_range.contains(&res.optimal_nbasis));
    }

    #[test]
    fn test_basis_nbasis_cv_fourier_cv() {
        let m = 60;
        let t = uniform_grid(m);
        let n = 10;
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin() + 0.02 * ((i * 11 + j) % 15) as f64;
            }
        }
        let nbasis_range: Vec<usize> = vec![5, 7, 9];
        let res = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Fourier { period: 1.0 },
            BasisCriterion::Cv,
            5,
            1e-4,
        )
        .unwrap();
        assert!(nbasis_range.contains(&res.optimal_nbasis));
        assert_eq!(res.criterion, BasisCriterion::Cv);
    }

    #[test]
    fn test_basis_nbasis_cv_with_nbasis_below_minimum() {
        // Range includes nbasis = 1 which is invalid
        let (data, t) = make_test_data(5, 50);
        let nbasis_range: Vec<usize> = vec![1, 5, 10];
        let res = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Gcv,
            5,
            1e-4,
        )
        .unwrap();
        // Score for nbasis=1 should be infinity, so optimal should be 5 or 10
        assert!(
            res.optimal_nbasis >= 5,
            "Should skip invalid nbasis=1, got optimal={}",
            res.optimal_nbasis
        );
        assert!(res.scores[0].is_infinite());
    }

    #[test]
    fn test_basis_nbasis_cv_empty_range() {
        let (data, t) = make_test_data(5, 50);
        let nbasis_range: Vec<usize> = vec![];
        let result = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Gcv,
            5,
            1e-4,
        );
        assert!(result.is_none());
    }

    #[test]
    fn test_basis_nbasis_cv_empty_data() {
        let data = FdMatrix::zeros(0, 50);
        let t = uniform_grid(50);
        let nbasis_range: Vec<usize> = vec![5, 10];
        let result = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Gcv,
            5,
            1e-4,
        );
        assert!(result.is_none());
    }

    #[test]
    fn test_basis_nbasis_cv_mismatched_argvals() {
        let data = FdMatrix::zeros(5, 50);
        let t = uniform_grid(40); // mismatch
        let nbasis_range: Vec<usize> = vec![5, 10];
        let result = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Gcv,
            5,
            1e-4,
        );
        assert!(result.is_none());
    }

    #[test]
    fn test_basis_nbasis_cv_single_nbasis() {
        let (data, t) = make_test_data(5, 50);
        let nbasis_range: Vec<usize> = vec![10];
        let res = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Gcv,
            5,
            1e-4,
        )
        .unwrap();
        assert_eq!(res.optimal_nbasis, 10);
        assert_eq!(res.scores.len(), 1);
    }

    #[test]
    fn test_basis_nbasis_cv_bic_penalizes_more_than_aic() {
        // BIC penalizes complexity more heavily than AIC, so it should generally
        // select the same or fewer basis functions
        let (data, t) = make_test_data(5, 80);
        let nbasis_range: Vec<usize> = (4..=20).step_by(2).collect();

        let aic_res = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Aic,
            5,
            1e-4,
        )
        .unwrap();
        let bic_res = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Bic,
            5,
            1e-4,
        )
        .unwrap();
        // BIC should select at most as many basis functions as AIC
        // (not guaranteed in all cases, but typical behavior)
        assert!(
            bic_res.optimal_nbasis <= aic_res.optimal_nbasis + 4,
            "BIC selected {} vs AIC selected {} -- BIC should not select much more than AIC",
            bic_res.optimal_nbasis,
            aic_res.optimal_nbasis
        );
    }

    // ─── Fitted values quality tests ────────────────────────────────────────

    #[test]
    fn test_smooth_basis_fitted_close_to_data() {
        // With moderate penalty and enough basis functions, fitted should be close to data
        let m = 50;
        let n = 3;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin();
            }
        }
        let fdpar = make_bspline_fdpar(&t, 15, 1e-6);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();

        let mut max_err = 0.0_f64;
        for i in 0..n {
            for j in 0..m {
                let err = (data[(i, j)] - res.fitted[(i, j)]).abs();
                max_err = max_err.max(err);
            }
        }
        assert!(
            max_err < 0.1,
            "Fitted should be close to smooth data; max_err={}",
            max_err
        );
    }

    #[test]
    fn test_smooth_basis_constant_data() {
        // Constant data should be fit exactly
        let m = 50;
        let n = 2;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] = 3.15;
            }
        }
        let fdpar = make_bspline_fdpar(&t, 10, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        for i in 0..n {
            for j in 0..m {
                assert!(
                    (res.fitted[(i, j)] - 3.15).abs() < 0.01,
                    "Constant data should be fit well at ({},{}): got {}",
                    i,
                    j,
                    res.fitted[(i, j)]
                );
            }
        }
    }

    #[test]
    fn test_smooth_basis_linear_data() {
        // Linear data should be fit well with cubic B-splines
        let m = 50;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(1, m);
        for j in 0..m {
            data[(0, j)] = 2.0 * t[j] + 1.0;
        }
        let fdpar = make_bspline_fdpar(&t, 10, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        for j in 0..m {
            let expected = 2.0 * t[j] + 1.0;
            assert!(
                (res.fitted[(0, j)] - expected).abs() < 0.05,
                "Linear data should be fit well at j={}: got {}, expected {}",
                j,
                res.fitted[(0, j)],
                expected
            );
        }
    }

    // ─── EDF and diagnostic tests ───────────────────────────────────────────

    #[test]
    fn test_smooth_basis_edf_bounded() {
        let m = 50;
        let (data, t) = make_test_data(3, m);
        let fdpar = make_bspline_fdpar(&t, 12, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        // EDF should be between 1 and m (evaluation points)
        assert!(
            res.edf > 0.0 && res.edf <= m as f64,
            "EDF should be in (0, {}]; got {}",
            m,
            res.edf
        );
    }

    #[test]
    fn test_smooth_basis_gcv_aic_bic_all_finite() {
        let (data, t) = make_test_data(4, 60);
        let fdpar = make_bspline_fdpar(&t, 12, 1e-3);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        assert!(res.gcv.is_finite(), "GCV should be finite: {}", res.gcv);
        assert!(res.aic.is_finite(), "AIC should be finite: {}", res.aic);
        assert!(res.bic.is_finite(), "BIC should be finite: {}", res.bic);
    }

    // ─── Penalty matrix size consistency tests ──────────────────────────────

    #[test]
    fn test_smooth_basis_penalty_matrix_in_result() {
        let (data, t) = make_test_data(3, 50);
        let nbasis = 10;
        let fdpar = make_bspline_fdpar(&t, nbasis, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        let k = res.nbasis;
        assert_eq!(
            res.penalty_matrix.len(),
            k * k,
            "Penalty matrix should be k*k = {}*{} = {}; got {}",
            k,
            k,
            k * k,
            res.penalty_matrix.len()
        );
    }

    // ─── Regression: multiple identical curves ──────────────────────────────

    #[test]
    fn test_smooth_basis_identical_curves_same_coefficients() {
        let m = 50;
        let t = uniform_grid(m);
        let curve: Vec<f64> = (0..m).map(|j| (2.0 * PI * t[j]).sin()).collect();
        let n = 4;
        let mut data = FdMatrix::zeros(n, m);
        for i in 0..n {
            for j in 0..m {
                data[(i, j)] = curve[j];
            }
        }
        let fdpar = make_bspline_fdpar(&t, 10, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();

        // All curves should have the same coefficients
        let k = res.coefficients.ncols();
        for i in 1..n {
            for j in 0..k {
                assert!(
                    (res.coefficients[(i, j)] - res.coefficients[(0, j)]).abs() < 1e-10,
                    "Identical curves should have identical coefficients: curve {} col {} differs",
                    i,
                    j
                );
            }
        }
    }

    // ─── Cross-validation: different numbers of folds ───────────────────────

    #[test]
    fn test_basis_nbasis_cv_different_nfolds() {
        let (data, t) = make_test_data(12, 50);
        let nbasis_range: Vec<usize> = vec![5, 8, 11];
        for nfolds in [2, 3, 5, 10] {
            let res = basis_nbasis_cv(
                &data,
                &t,
                &nbasis_range,
                &BasisType::Bspline { order: 4 },
                BasisCriterion::Cv,
                nfolds,
                1e-4,
            );
            assert!(res.is_some(), "CV should succeed with nfolds={}", nfolds);
            let r = res.unwrap();
            assert!(nbasis_range.contains(&r.optimal_nbasis));
        }
    }

    // ─── Large nbasis / more basis than reasonable ──────────────────────────

    #[test]
    fn test_smooth_basis_many_basis_functions() {
        let m = 100;
        let (data, t) = make_test_data(2, m);
        // Many basis functions relative to data points
        let fdpar = make_bspline_fdpar(&t, 40, 1e-2);
        let res = smooth_basis(&data, &t, &fdpar);
        assert!(
            res.is_ok(),
            "Should handle many basis functions with penalty"
        );
    }

    // ─── evaluate_basis internal function (indirectly tested) ───────────────

    #[test]
    fn test_smooth_basis_bspline_vs_fourier_different_results() {
        let m = 50;
        let (data, t) = make_test_data(2, m);
        let fdpar_bs = make_bspline_fdpar(&t, 9, 1e-4);
        let fdpar_f = make_fourier_fdpar(9, 1.0, 1e-4);
        let res_bs = smooth_basis(&data, &t, &fdpar_bs).unwrap();
        let res_f = smooth_basis(&data, &t, &fdpar_f).unwrap();
        // Results should differ between the two basis types
        let diff: f64 = (0..m)
            .map(|j| (res_bs.fitted[(0, j)] - res_f.fitted[(0, j)]).abs())
            .sum();
        // They fit the same data, so some difference is expected but not huge
        assert!(
            diff > 1e-10,
            "B-spline and Fourier fits should differ for the same data"
        );
    }

    // ─── compute_gcv edge cases (indirectly tested) ─────────────────────────

    #[test]
    fn test_smooth_basis_gcv_positive_for_noisy_data() {
        let m = 50;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(1, m);
        for j in 0..m {
            // Noisy data
            data[(0, j)] = (2.0 * PI * t[j]).sin() + 0.5 * ((j * 37) % 20) as f64 / 20.0 - 0.25;
        }
        let fdpar = make_bspline_fdpar(&t, 10, 1e-3);
        let res = smooth_basis(&data, &t, &fdpar).unwrap();
        assert!(res.gcv > 0.0, "GCV should be positive for noisy data");
    }

    // ─── Penalty order (lfd_order) tests ────────────────────────────────────

    #[test]
    fn test_smooth_basis_different_lfd_orders() {
        let m = 50;
        let (data, t) = make_test_data(2, m);

        // lfd_order = 1 (penalize first derivative)
        let penalty1 = bspline_penalty_matrix(&t, 10, 4, 1);
        let fdpar1 = FdPar {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis: 10,
            lambda: 1e-2,
            lfd_order: 1,
            penalty_matrix: penalty1,
        };
        let res1 = smooth_basis(&data, &t, &fdpar1);
        assert!(res1.is_ok());

        // lfd_order = 2 (penalize second derivative)
        let penalty2 = bspline_penalty_matrix(&t, 10, 4, 2);
        let fdpar2 = FdPar {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis: 10,
            lambda: 1e-2,
            lfd_order: 2,
            penalty_matrix: penalty2,
        };
        let res2 = smooth_basis(&data, &t, &fdpar2);
        assert!(res2.is_ok());

        // Different penalty orders should produce different fitted values
        let r1 = res1.unwrap();
        let r2 = res2.unwrap();
        let diff: f64 = (0..m)
            .map(|j| (r1.fitted[(0, j)] - r2.fitted[(0, j)]).abs())
            .sum();
        assert!(
            diff > 1e-10,
            "Different lfd_orders should produce different fits"
        );
    }

    // ─── BasisNbasisCvResult field tests ────────────────────────────────────

    #[test]
    fn test_basis_nbasis_cv_result_fields() {
        let (data, t) = make_test_data(6, 50);
        let nbasis_range: Vec<usize> = vec![5, 7, 9, 11, 13];
        let res = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Aic,
            5,
            1e-4,
        )
        .unwrap();

        assert!(nbasis_range.contains(&res.optimal_nbasis));
        assert_eq!(res.scores.len(), nbasis_range.len());
        assert_eq!(res.nbasis_range, nbasis_range);
        assert_eq!(res.criterion, BasisCriterion::Aic);
        // optimal_nbasis should correspond to minimum score
        let min_score = res.scores.iter().copied().fold(f64::INFINITY, f64::min);
        let best_idx = res
            .scores
            .iter()
            .position(|&s| (s - min_score).abs() < 1e-15)
            .unwrap();
        assert_eq!(res.optimal_nbasis, nbasis_range[best_idx]);
    }

    #[test]
    fn test_basis_nbasis_cv_result_clone() {
        let (data, t) = make_test_data(5, 50);
        let nbasis_range: Vec<usize> = vec![5, 10];
        let res = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &BasisType::Bspline { order: 4 },
            BasisCriterion::Gcv,
            5,
            1e-4,
        )
        .unwrap();
        let cloned = res.clone();
        assert_eq!(res, cloned);
    }

    // ─── Non-uniform argvals ────────────────────────────────────────────────

    #[test]
    fn test_smooth_basis_nonuniform_argvals() {
        let m = 50;
        // Non-uniform grid: denser at the ends
        let t: Vec<f64> = (0..m)
            .map(|i| {
                let x = i as f64 / (m - 1) as f64;
                0.5 * (1.0 - (PI * x).cos())
            })
            .collect();
        let mut data = FdMatrix::zeros(2, m);
        for i in 0..2 {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin() + 0.1 * i as f64;
            }
        }
        let fdpar = make_bspline_fdpar(&t, 10, 1e-4);
        let res = smooth_basis(&data, &t, &fdpar);
        assert!(res.is_ok(), "Should handle non-uniform argvals");
        let r = res.unwrap();
        assert_eq!(r.fitted.shape(), (2, m));
    }

    // ─── Numerical stability with extreme lambda ────────────────────────────

    #[test]
    fn test_smooth_basis_very_small_lambda() {
        let m = 50;
        let (data, t) = make_test_data(2, m);
        let fdpar = make_bspline_fdpar(&t, 10, 1e-15);
        let res = smooth_basis(&data, &t, &fdpar);
        assert!(res.is_ok(), "Should handle very small lambda");
    }

    #[test]
    fn test_smooth_basis_very_large_lambda() {
        let m = 50;
        let (data, t) = make_test_data(2, m);
        let fdpar = make_bspline_fdpar(&t, 10, 1e10);
        let res = smooth_basis(&data, &t, &fdpar);
        assert!(res.is_ok(), "Should handle very large lambda");
    }

    // ─── Multiple curves consistency ────────────────────────────────────────

    #[test]
    fn test_smooth_basis_multi_curve_vs_single_curve() {
        // Smoothing multiple curves at once should give the same result as smoothing each individually
        let m = 50;
        let n = 3;
        let (data, t) = make_test_data(n, m);
        let fdpar = make_bspline_fdpar(&t, 10, 1e-3);

        // All at once
        let res_all = smooth_basis(&data, &t, &fdpar).unwrap();

        // One at a time
        for i in 0..n {
            let mut single = FdMatrix::zeros(1, m);
            for j in 0..m {
                single[(0, j)] = data[(i, j)];
            }
            let res_single = smooth_basis(&single, &t, &fdpar).unwrap();
            for j in 0..m {
                assert!(
                    (res_all.fitted[(i, j)] - res_single.fitted[(0, j)]).abs() < 1e-10,
                    "Multi-curve fit should match single-curve fit: curve {} point {}",
                    i,
                    j
                );
            }
        }
    }

    // ─── BasisCriterion comparison: all criteria produce finite scores ──────

    #[test]
    fn test_basis_nbasis_cv_all_criteria_finite_scores() {
        let (data, t) = make_test_data(10, 60);
        let nbasis_range: Vec<usize> = vec![5, 7, 9, 11];

        for criterion in [
            BasisCriterion::Gcv,
            BasisCriterion::Aic,
            BasisCriterion::Bic,
            BasisCriterion::Cv,
        ] {
            let res = basis_nbasis_cv(
                &data,
                &t,
                &nbasis_range,
                &BasisType::Bspline { order: 4 },
                criterion,
                5,
                1e-4,
            )
            .unwrap();
            // At least some scores should be finite (valid nbasis values)
            let finite_count = res.scores.iter().filter(|s| s.is_finite()).count();
            assert!(
                finite_count > 0,
                "At least one score should be finite for {:?}",
                criterion
            );
        }
    }

    // ─── SmoothBasisGcvConfig tests ────────────────────────────────────────

    #[test]
    fn test_smooth_basis_gcv_config_default() {
        let config = SmoothBasisGcvConfig::default();
        assert_eq!(config.basis_type, BasisType::Bspline { order: 4 });
        assert_eq!(config.nbasis, 15);
        assert_eq!(config.lfd_order, 2);
        assert_eq!(config.log_lambda_range, (-10.0, 2.0));
        assert_eq!(config.n_grid, 50);
    }

    #[test]
    fn test_smooth_basis_gcv_config_clone_eq() {
        let config = SmoothBasisGcvConfig {
            nbasis: 20,
            ..SmoothBasisGcvConfig::default()
        };
        let cloned = config.clone();
        assert_eq!(config, cloned);
    }

    #[test]
    fn test_smooth_basis_gcv_config_debug() {
        let config = SmoothBasisGcvConfig::default();
        let debug_str = format!("{:?}", config);
        assert!(debug_str.contains("SmoothBasisGcvConfig"));
        assert!(debug_str.contains("nbasis"));
    }

    #[test]
    fn test_smooth_basis_gcv_config_partial_override() {
        let config = SmoothBasisGcvConfig {
            basis_type: BasisType::Fourier { period: 2.0 },
            n_grid: 100,
            ..SmoothBasisGcvConfig::default()
        };
        assert_eq!(config.basis_type, BasisType::Fourier { period: 2.0 });
        assert_eq!(config.n_grid, 100);
        // defaults preserved
        assert_eq!(config.nbasis, 15);
        assert_eq!(config.lfd_order, 2);
    }

    #[test]
    fn test_smooth_basis_gcv_with_config_default() {
        let (data, t) = make_test_data(5, 101);
        let config = SmoothBasisGcvConfig::default();
        let result = smooth_basis_gcv_with_config(&data, &t, &config);
        assert!(result.is_ok(), "GCV with default config should succeed");
        let res = result.unwrap();
        assert_eq!(res.fitted.shape(), (5, 101));
        assert!(res.edf > 0.0);
        assert!(res.gcv.is_finite());
    }

    #[test]
    fn test_smooth_basis_gcv_with_config_custom() {
        let (data, t) = make_test_data(3, 50);
        let config = SmoothBasisGcvConfig {
            nbasis: 10,
            log_lambda_range: (-6.0, 0.0),
            n_grid: 15,
            ..SmoothBasisGcvConfig::default()
        };
        let result = smooth_basis_gcv_with_config(&data, &t, &config);
        assert!(result.is_ok());
    }

    #[test]
    fn test_smooth_basis_gcv_with_config_matches_direct() {
        let (data, t) = make_test_data(3, 50);
        let config = SmoothBasisGcvConfig {
            nbasis: 10,
            log_lambda_range: (-6.0, 0.0),
            n_grid: 20,
            ..SmoothBasisGcvConfig::default()
        };
        let with_config = smooth_basis_gcv_with_config(&data, &t, &config).unwrap();
        let direct = smooth_basis_gcv(
            &data,
            &t,
            &config.basis_type,
            config.nbasis,
            config.lfd_order,
            config.log_lambda_range,
            config.n_grid,
        )
        .unwrap();
        assert_eq!(with_config.gcv, direct.gcv);
        assert_eq!(with_config.edf, direct.edf);
        assert_eq!(with_config.nbasis, direct.nbasis);
    }

    #[test]
    fn test_smooth_basis_gcv_with_config_fourier() {
        let m = 100;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(2, m);
        for i in 0..2 {
            for j in 0..m {
                data[(i, j)] = (2.0 * PI * t[j]).sin() + (4.0 * PI * t[j]).cos();
            }
        }
        let config = SmoothBasisGcvConfig {
            basis_type: BasisType::Fourier { period: 1.0 },
            nbasis: 7,
            n_grid: 20,
            ..SmoothBasisGcvConfig::default()
        };
        let result = smooth_basis_gcv_with_config(&data, &t, &config);
        assert!(result.is_ok());
    }

    // ─── BasisNbasisCvConfig tests ─────────────────────────────────────────

    #[test]
    fn test_basis_nbasis_cv_config_default() {
        let config = BasisNbasisCvConfig::default();
        assert_eq!(config.basis_type, BasisType::Bspline { order: 4 });
        assert_eq!(config.nbasis_range, (5, 30));
        assert!((config.lambda - 1e-4).abs() < 1e-15);
        assert_eq!(config.lfd_order, 2);
        assert_eq!(config.n_folds, 5);
        assert_eq!(config.criterion, BasisCriterion::Gcv);
    }

    #[test]
    fn test_basis_nbasis_cv_config_clone_eq() {
        let config = BasisNbasisCvConfig {
            nbasis_range: (4, 15),
            ..BasisNbasisCvConfig::default()
        };
        let cloned = config.clone();
        assert_eq!(config, cloned);
    }

    #[test]
    fn test_basis_nbasis_cv_config_debug() {
        let config = BasisNbasisCvConfig::default();
        let debug_str = format!("{:?}", config);
        assert!(debug_str.contains("BasisNbasisCvConfig"));
        assert!(debug_str.contains("nbasis_range"));
    }

    #[test]
    fn test_basis_nbasis_cv_config_partial_override() {
        let config = BasisNbasisCvConfig {
            criterion: BasisCriterion::Aic,
            lambda: 1e-2,
            ..BasisNbasisCvConfig::default()
        };
        assert_eq!(config.criterion, BasisCriterion::Aic);
        assert!((config.lambda - 1e-2).abs() < 1e-15);
        // defaults preserved
        assert_eq!(config.nbasis_range, (5, 30));
        assert_eq!(config.n_folds, 5);
    }

    #[test]
    fn test_basis_nbasis_cv_with_config_default() {
        let (data, t) = make_test_data(5, 51);
        let config = BasisNbasisCvConfig {
            nbasis_range: (5, 12),
            ..BasisNbasisCvConfig::default()
        };
        let result = basis_nbasis_cv_with_config(&data, &t, &config);
        assert!(
            result.is_ok(),
            "nbasis CV with default config should succeed"
        );
        let res = result.unwrap();
        assert!(res.optimal_nbasis >= 5 && res.optimal_nbasis <= 12);
        assert_eq!(res.scores.len(), 8); // 5..=12 = 8 values
        assert_eq!(res.criterion, BasisCriterion::Gcv);
    }

    #[test]
    fn test_basis_nbasis_cv_with_config_aic() {
        let (data, t) = make_test_data(5, 51);
        let config = BasisNbasisCvConfig {
            nbasis_range: (5, 10),
            criterion: BasisCriterion::Aic,
            ..BasisNbasisCvConfig::default()
        };
        let result = basis_nbasis_cv_with_config(&data, &t, &config);
        assert!(result.is_ok());
        assert_eq!(result.unwrap().criterion, BasisCriterion::Aic);
    }

    #[test]
    fn test_basis_nbasis_cv_with_config_cv_folds() {
        let (data, t) = make_test_data(10, 51);
        let config = BasisNbasisCvConfig {
            nbasis_range: (5, 9),
            criterion: BasisCriterion::Cv,
            n_folds: 3,
            ..BasisNbasisCvConfig::default()
        };
        let result = basis_nbasis_cv_with_config(&data, &t, &config);
        assert!(result.is_ok());
        assert_eq!(result.unwrap().criterion, BasisCriterion::Cv);
    }

    #[test]
    fn test_basis_nbasis_cv_with_config_matches_direct() {
        let (data, t) = make_test_data(5, 51);
        let config = BasisNbasisCvConfig {
            nbasis_range: (5, 10),
            criterion: BasisCriterion::Bic,
            lambda: 1e-3,
            ..BasisNbasisCvConfig::default()
        };
        let with_config = basis_nbasis_cv_with_config(&data, &t, &config).unwrap();
        let nbasis_range: Vec<usize> = (5..=10).collect();
        let direct = basis_nbasis_cv(
            &data,
            &t,
            &nbasis_range,
            &config.basis_type,
            config.criterion,
            config.n_folds,
            config.lambda,
        )
        .unwrap();
        assert_eq!(with_config.optimal_nbasis, direct.optimal_nbasis);
        assert_eq!(with_config.scores, direct.scores);
        assert_eq!(with_config.nbasis_range, direct.nbasis_range);
    }

    #[test]
    fn test_basis_nbasis_cv_with_config_nbasis_range_expansion() {
        let (data, t) = make_test_data(5, 51);
        let config = BasisNbasisCvConfig {
            nbasis_range: (7, 7), // single value
            ..BasisNbasisCvConfig::default()
        };
        let result = basis_nbasis_cv_with_config(&data, &t, &config);
        assert!(result.is_ok());
        let res = result.unwrap();
        assert_eq!(res.optimal_nbasis, 7);
        assert_eq!(res.scores.len(), 1);
    }

    // ─── smooth_positive tests ──────────────────────────────────────────────

    /// Build a shared FdPar for positive-smoother tests (B-spline, nbasis=10, lambda=1e-3).
    fn make_positive_fdpar(argvals: &[f64]) -> FdPar {
        let nbasis = 10;
        let penalty = bspline_penalty_matrix(argvals, nbasis, 4, 2);
        FdPar {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis,
            lambda: 1e-3,
            lfd_order: 2,
            penalty_matrix: penalty,
        }
    }

    #[test]
    fn test_smooth_positive_is_positive() {
        // Positive signal: 2 + sin(2π t) stays in (1, 3) — strictly positive.
        let m = 41;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(1, m);
        for j in 0..m {
            // add a tiny deterministic wiggle so the curve is not constant
            let wiggle = 0.05 * ((j * 7) % 13) as f64 / 13.0;
            data[(0, j)] = 2.0 + (2.0 * PI * t[j]).sin() + wiggle;
        }

        let fdpar = make_positive_fdpar(&t);
        let result = smooth_positive(&data, &t, &fdpar);
        assert!(
            result.is_ok(),
            "smooth_positive should succeed on positive data"
        );

        let res = result.unwrap();
        assert_eq!(res.fitted.shape(), (1, m));
        assert_eq!(res.log_coefficients.nrows(), 1);

        for j in 0..m {
            let v = res.fitted[(0, j)];
            assert!(v > 0.0, "fitted value at j={} is not positive: {}", j, v);
            assert!(
                v.is_finite(),
                "fitted value at j={} is not finite: {}",
                j,
                v
            );
        }
    }

    #[test]
    fn test_smooth_positive_recovers_curve() {
        // With a small lambda, the smoother should closely track the true positive curve.
        let m = 41;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(1, m);
        let mut truth = vec![0.0f64; m];
        for j in 0..m {
            truth[j] = 2.0 + (2.0 * PI * t[j]).sin();
            data[(0, j)] = truth[j]; // no noise — test recovery directly
        }

        let nbasis = 10;
        let penalty = bspline_penalty_matrix(&t, nbasis, 4, 2);
        let fdpar = FdPar {
            basis_type: BasisType::Bspline { order: 4 },
            nbasis,
            lambda: 1e-6, // small lambda → close interpolation
            lfd_order: 2,
            penalty_matrix: penalty,
        };

        let res = smooth_positive(&data, &t, &fdpar).unwrap();
        let mae: f64 = (0..m)
            .map(|j| (res.fitted[(0, j)] - truth[j]).abs())
            .sum::<f64>()
            / m as f64;

        assert!(
            mae < 0.2,
            "Mean absolute error too large: {}; smooth_positive should recover positive curve",
            mae
        );
        // EDF and GCV should be reasonable.
        assert!(res.edf > 0.0, "EDF should be positive");
        assert!(res.gcv.is_finite(), "GCV should be finite");
    }

    #[test]
    fn test_smooth_positive_rejects_nonpositive() {
        // A data matrix with one zero element should be rejected with InvalidParameter.
        let m = 41;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(1, m);
        for j in 0..m {
            data[(0, j)] = 2.0 + (2.0 * PI * t[j]).sin();
        }
        // Inject a non-positive value at position (0, 10).
        data[(0, 10)] = 0.0;

        let fdpar = make_positive_fdpar(&t);
        let result = smooth_positive(&data, &t, &fdpar);
        assert!(
            result.is_err(),
            "smooth_positive must reject data with a zero element"
        );

        match result.unwrap_err() {
            crate::FdarError::InvalidParameter { parameter, .. } => {
                assert_eq!(parameter, "data");
            }
            other => panic!("Expected InvalidParameter, got {:?}", other),
        }
    }

    #[test]
    fn test_smooth_positive_rejects_negative() {
        // A data matrix with one negative element should also be rejected.
        let m = 41;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(1, m);
        for j in 0..m {
            data[(0, j)] = 2.0 + (2.0 * PI * t[j]).sin();
        }
        data[(0, 20)] = -0.5;

        let fdpar = make_positive_fdpar(&t);
        let result = smooth_positive(&data, &t, &fdpar);
        assert!(
            result.is_err(),
            "smooth_positive must reject data with a negative element"
        );

        match result.unwrap_err() {
            crate::FdarError::InvalidParameter { parameter, .. } => {
                assert_eq!(parameter, "data");
            }
            other => panic!("Expected InvalidParameter, got {:?}", other),
        }
    }

    // ============== smooth_monotone tests (Task 1: tracer) ==============

    #[test]
    fn test_smooth_monotone_is_monotone() {
        // Fit a monotone-increasing target: g(t) = t².
        let m = 41_usize;
        let t = uniform_grid(m); // [0, 1] uniform
        let y: Vec<f64> = t.iter().map(|&ti| ti * ti).collect();

        let result = smooth_monotone(&y, &t, 8, 4, 1e-3, 50)
            .expect("smooth_monotone should succeed on t² data");

        // Every fitted value must be finite.
        for (i, &v) in result.fitted.iter().enumerate() {
            assert!(v.is_finite(), "fitted[{}] = {} is not finite", i, v);
        }

        // Structural nondecreasing guarantee must hold regardless of convergence.
        for i in 1..m {
            assert!(
                result.fitted[i] >= result.fitted[i - 1] - 1e-9,
                "Monotonicity violated at i={}: fitted[{}]={} < fitted[{}]={}",
                i,
                i,
                result.fitted[i],
                i - 1,
                result.fitted[i - 1]
            );
        }

        // beta1 should be positive (increasing direction).
        assert!(
            result.beta1 > 0.0,
            "beta1 should be positive for increasing data, got {}",
            result.beta1
        );
    }

    // ============== smooth_monotone tests (Task 2: recovery + direction) ==============

    #[test]
    fn test_smooth_monotone_recovers_increasing() {
        // Logistic target: well-conditioned smooth increasing function on [0,1].
        // The logistic has high curvature near the inflection; the Gauss-Newton scheme
        // requires more iterations to shape the w(u) spline precisely.  We use
        // max_iter=100 so convergence can complete; the bounded-iterations oracle
        // (test_smooth_monotone_bounded_iterations) uses max_iter=50 separately.
        let m = 51_usize;
        let t = uniform_grid(m);
        let y: Vec<f64> = t
            .iter()
            .map(|&ti| 1.0 / (1.0 + (-8.0 * (ti - 0.5)).exp()))
            .collect();

        let result = smooth_monotone(&y, &t, 10, 4, 1e-4, 100)
            .expect("smooth_monotone should succeed on logistic data");

        // Recovery within tolerance.  Tolerance is 0.15 (not a tight 0.05) to allow
        // for boundary effects in cumulative-trapezoid integration and B-spline knot
        // placement at the grid edges — both standard sources of bias in v1.
        let mae: f64 = y
            .iter()
            .zip(result.fitted.iter())
            .map(|(&yi, &fi)| (yi - fi).abs())
            .sum::<f64>()
            / m as f64;
        assert!(
            mae < 0.15,
            "Mean absolute error {} too large for logistic recovery (tolerance 0.15, iterations={})",
            mae,
            result.iterations
        );

        // Monotonicity must still hold regardless of convergence.
        for i in 1..m {
            assert!(
                result.fitted[i] >= result.fitted[i - 1] - 1e-9,
                "Monotonicity violated at i={}: fitted[{}]={} < fitted[{}]={}",
                i,
                result.fitted[i],
                i,
                result.fitted[i - 1],
                i - 1
            );
        }
    }

    #[test]
    fn test_smooth_monotone_decreasing() {
        // Decreasing target: g(t) = 1 - t (simple linear decrease).
        let m = 41_usize;
        let t = uniform_grid(m);
        let y: Vec<f64> = t.iter().map(|&ti| 1.0 - ti).collect();

        let result = smooth_monotone(&y, &t, 8, 4, 1e-3, 50)
            .expect("smooth_monotone should succeed on decreasing data");

        // Direction auto-detected: beta1 must be negative.
        assert!(
            result.beta1 < 0.0,
            "beta1 should be negative for decreasing data, got {}",
            result.beta1
        );

        // Fit must be nonincreasing: fitted[i] <= fitted[i-1] + 1e-9.
        for i in 1..m {
            assert!(
                result.fitted[i] <= result.fitted[i - 1] + 1e-9,
                "Nonincreasing violated at i={}: fitted[{}]={} > fitted[{}]={}",
                i,
                result.fitted[i],
                i,
                result.fitted[i - 1],
                i - 1
            );
        }
    }

    #[test]
    fn test_smooth_monotone_bounded_iterations() {
        // Noisy monotone data (deterministic pseudo-noise so test is reproducible).
        let m = 51_usize;
        let t = uniform_grid(m);
        let y: Vec<f64> = t
            .iter()
            .enumerate()
            .map(|(i, &ti)| {
                let noise = 0.05 * ((i * 17 + 3) % 11) as f64 / 10.0 - 0.025;
                ti + noise // noisy increasing; not guaranteed sorted but trend is up
            })
            .collect();

        let result = smooth_monotone(&y, &t, 8, 4, 1e-2, 50)
            .expect("smooth_monotone should succeed on noisy data");

        // Bounded iteration count.
        assert!(
            result.iterations <= 50,
            "iterations ({}) must be <= max_iter (50)",
            result.iterations
        );

        // Structural monotonicity holds even if not converged.
        // Increasing data → beta1 > 0 → fitted is nondecreasing.
        // (beta1 < 0 would mean fitted is nonincreasing — check the appropriate direction.)
        if result.beta1 >= 0.0 {
            for i in 1..m {
                assert!(
                    result.fitted[i] >= result.fitted[i - 1] - 1e-9,
                    "Nondecreasing violated (beta1={}) at i={}",
                    result.beta1,
                    i
                );
            }
        } else {
            for i in 1..m {
                assert!(
                    result.fitted[i] <= result.fitted[i - 1] + 1e-9,
                    "Nonincreasing violated (beta1={}) at i={}",
                    result.beta1,
                    i
                );
            }
        }
    }

    // ============== smooth_monotone tests (Task 3: error paths) ==============

    #[test]
    fn test_smooth_monotone_errors_on_short_input() {
        // data.len() == 2 → need >= 3.
        let data = vec![0.0, 1.0];
        let argvals = vec![0.0, 1.0];
        let result = smooth_monotone(&data, &argvals, 4, 4, 1e-3, 50);
        assert!(
            result.is_err(),
            "smooth_monotone should fail for data.len() == 2"
        );
        match result.unwrap_err() {
            crate::FdarError::InvalidDimension { parameter, .. } => {
                assert!(
                    parameter.contains("data") || parameter.contains("argvals"),
                    "Expected data/argvals dimension error, got param={}",
                    parameter
                );
            }
            other => panic!("Expected InvalidDimension, got {:?}", other),
        }
    }

    #[test]
    fn test_smooth_monotone_errors_on_argvals_mismatch() {
        // argvals.len() != data.len()
        let m = 10_usize;
        let data: Vec<f64> = (0..m).map(|i| i as f64).collect();
        let argvals: Vec<f64> = (0..m - 1).map(|i| i as f64).collect();
        let result = smooth_monotone(&data, &argvals, 4, 4, 1e-3, 50);
        assert!(
            result.is_err(),
            "smooth_monotone should fail for argvals length mismatch"
        );
        match result.unwrap_err() {
            crate::FdarError::InvalidDimension { .. } => {}
            other => panic!("Expected InvalidDimension, got {:?}", other),
        }
    }

    #[test]
    fn test_smooth_monotone_errors_on_bad_params() {
        let m = 10_usize;
        let t = uniform_grid(m);
        let data: Vec<f64> = t.clone();

        // nbasis == 1 → InvalidParameter
        let result = smooth_monotone(&data, &t, 1, 4, 1e-3, 50);
        assert!(
            result.is_err(),
            "smooth_monotone should fail for nbasis == 1"
        );
        match result.unwrap_err() {
            crate::FdarError::InvalidParameter { parameter, .. } => {
                assert_eq!(parameter, "nbasis");
            }
            other => panic!("Expected InvalidParameter(nbasis), got {:?}", other),
        }

        // max_iter == 0 → InvalidParameter
        let result = smooth_monotone(&data, &t, 8, 4, 1e-3, 0);
        assert!(
            result.is_err(),
            "smooth_monotone should fail for max_iter == 0"
        );
        match result.unwrap_err() {
            crate::FdarError::InvalidParameter { parameter, .. } => {
                assert_eq!(parameter, "max_iter");
            }
            other => panic!("Expected InvalidParameter(max_iter), got {:?}", other),
        }
    }
}