//! Term construction: bridge from parsed formula terms to `TermCollectionSpec`.
//!
//! This module takes the AST produced by `inference::formula_dsl` and a loaded
//! dataset, resolves column references, infers knot counts and center strategies,
//! and produces a `TermCollectionSpec` ready for `build_term_collection_design`.

use std::collections::{BTreeMap, BTreeSet, HashMap};
use std::path::PathBuf;

use ndarray::{Array2, ArrayView1};

use crate::basis::{
    BSplineBasisSpec, BSplineBoundaryConditions, BSplineEndpointBoundaryCondition,
    BSplineIdentifiability, BSplineKnotSpec, CenterCountRequest, CenterStrategy, DuchonBasisSpec,
    DuchonNullspaceOrder, DuchonOperatorPenaltySpec, MaternBasisSpec, MaternIdentifiability,
    MaternNu, OneDimensionalBoundary, SpatialIdentifiability, SphereMethod, SphereWahbaKernel,
    SphericalSplineBasisSpec, SphericalSplineIdentifiability, ThinPlateBasisSpec,
    auto_spatial_center_strategy, default_num_centers, default_spatial_center_strategy,
    default_spherical_harmonic_degree, plan_spatial_basis,
};
use crate::inference::data::{DataError, EncodedDataset as Dataset};
use crate::inference::formula_dsl::{
    ParsedTerm, SmoothKind, option_bool, option_f64, option_f64_strict, option_usize,
    option_usize_any, option_usize_any_strict, option_usize_strict, strip_quotes,
};
use crate::inference::model::ColumnKindTag;
use crate::resource::ResourcePolicy;
use crate::smooth::{
    BySmoothKind, ByVarKind, ByVariableSpec, FactorSmoothFlavour, FactorSmoothSpec,
    LinearCoefficientGeometry, LinearTermSpec, RandomEffectTermSpec, ShapeConstraint,
    SmoothBasisSpec, SmoothTermSpec, TensorBSplineIdentifiability, TensorBSplineSpec,
    TermCollectionSpec,
};
use crate::types::ColIdx;

fn default_matern_length_scale(ds: &Dataset, cols: &[usize]) -> f64 {
    let mut diameter2 = 0.0_f64;
    for &col in cols {
        let column = ds.values.column(col);
        let mut lo = f64::INFINITY;
        let mut hi = f64::NEG_INFINITY;
        for &value in column.iter().filter(|v| v.is_finite()) {
            lo = lo.min(value);
            hi = hi.max(value);
        }
        if lo.is_finite() && hi.is_finite() && hi > lo {
            let span = hi - lo;
            diameter2 += span * span;
        }
    }
    let diameter = diameter2.sqrt();
    if diameter.is_finite() && diameter > 0.0 {
        (0.15 * diameter).max(1e-6)
    } else {
        1.0
    }
}

// ---------------------------------------------------------------------------
// Typed errors
// ---------------------------------------------------------------------------

/// Typed errors emitted by term-builder helpers. `Display` reproduces the exact
/// pre-refactor `format!(...)` text byte-for-byte, so callers that string-match
/// on the message (tests, log assertions) keep working unchanged. Public-API
/// functions still return `Result<_, String>` and use `.to_string()` shims at
/// their boundary to stay compatible with callers in protected modules.
#[derive(Clone, Debug)]
pub enum TermBuilderError {
    /// Column-resolution / column-kind lookup failures whose context is purely
    /// internal (column-kind table out-of-sync, alias map missing an entry,
    /// etc.). User-facing "this formula references a column that doesn't
    /// exist" diagnostics use the dedicated `ColumnNotFound` variant so the
    /// FFI boundary can lift the structured payload into a Python
    /// `ColumnNotFoundError` without parsing prose.
    MissingColumn { reason: String },
    /// A formula referenced a column that is not present in the input data.
    /// Mirrors `DataError::ColumnNotFound` field-for-field so the conversion
    /// across module boundaries is a pure data move (no re-derivation, no
    /// string re-parsing). Public callers see byte-identical `Display`
    /// output to the legacy `missing_column_message` text.
    ColumnNotFound {
        name: String,
        role: Option<String>,
        available: Vec<String>,
        similar: Vec<String>,
        tsv_hint: bool,
    },
    /// User-specified configuration is internally inconsistent (e.g. too few
    /// variables for a smooth type, conflicting size options, requested basis
    /// dimension below the polynomial nullspace).
    IncompatibleConfig { reason: String },
    /// Option parsing failure: malformed numeric expression, unknown option
    /// key, out-of-range integer, list-length mismatch, etc.
    InvalidOption { reason: String },
    /// User requested a feature that is intentionally not supported (unknown
    /// smooth type / method / kernel / identifiability, non-zero anchor,
    /// internal-only token, etc.).
    UnsupportedFeature { reason: String },
    /// Input data is degenerate for the requested term (constant column,
    /// non-finite categorical entries, ...).
    DegenerateData { reason: String },
    /// Term-collection-stage formula error — a node that the caller was
    /// supposed to resolve upstream reached the builder.
    MalformedFormula { reason: String },
}

impl std::fmt::Display for TermBuilderError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            TermBuilderError::MissingColumn { reason }
            | TermBuilderError::IncompatibleConfig { reason }
            | TermBuilderError::InvalidOption { reason }
            | TermBuilderError::UnsupportedFeature { reason }
            | TermBuilderError::DegenerateData { reason }
            | TermBuilderError::MalformedFormula { reason } => f.write_str(reason),
            // Delegate to the canonical `DataError::ColumnNotFound` formatter
            // so a single source of truth defines the human text. The
            // intermediate `DataError` constructed here owns its strings only
            // for the duration of the Display call — no allocation cost
            // beyond the original payload that this variant already holds.
            TermBuilderError::ColumnNotFound {
                name,
                role,
                available,
                similar,
                tsv_hint,
            } => {
                let canonical = DataError::ColumnNotFound {
                    name: name.clone(),
                    role: role.clone(),
                    available: available.clone(),
                    similar: similar.clone(),
                    tsv_hint: *tsv_hint,
                };
                std::fmt::Display::fmt(&canonical, f)
            }
        }
    }
}

impl From<TermBuilderError> for String {
    fn from(err: TermBuilderError) -> String {
        err.to_string()
    }
}

/// Catchall lift for the term-builder's internal `Result<_, String>` helpers
/// (numeric expression parsing, option lookup, boundary-condition parsing,
/// ...) that flow into `build_termspec` via `?`. Maps to
/// `IncompatibleConfig`, which is the most appropriate generic bucket for
/// option/config-style failures — leaf sites that emit structured payloads
/// (`From<DataError>` for column-not-found) bypass this fallback.
impl From<String> for TermBuilderError {
    fn from(reason: String) -> Self {
        Self::IncompatibleConfig { reason }
    }
}

/// Typed lift from data-layer errors. `DataError::ColumnNotFound` becomes
/// `TermBuilderError::ColumnNotFound` field-for-field — no stringification,
/// no information loss — so the FFI boundary downstream can dispatch on
/// the typed variant. Other `DataError` variants degrade into
/// `MissingColumn` since they describe column-resolution-time failures
/// without a dedicated structured destination.
impl From<DataError> for TermBuilderError {
    fn from(err: DataError) -> Self {
        match err {
            DataError::ColumnNotFound {
                name,
                role,
                available,
                similar,
                tsv_hint,
            } => Self::ColumnNotFound {
                name,
                role,
                available,
                similar,
                tsv_hint,
            },
            DataError::SchemaMismatch { reason }
            | DataError::ParseError { reason }
            | DataError::EncodingFailure { reason }
            | DataError::EmptyInput { reason }
            | DataError::InvalidValue { reason } => Self::MissingColumn { reason },
        }
    }
}

// Constructor helpers — keep error-site code compact and consistent.
impl TermBuilderError {
    #[inline]
    fn missing_column(reason: impl Into<String>) -> Self {
        TermBuilderError::MissingColumn {
            reason: reason.into(),
        }
    }
    #[inline]
    fn incompatible_config(reason: impl Into<String>) -> Self {
        TermBuilderError::IncompatibleConfig {
            reason: reason.into(),
        }
    }
    #[inline]
    fn invalid_option(reason: impl Into<String>) -> Self {
        TermBuilderError::InvalidOption {
            reason: reason.into(),
        }
    }
    #[inline]
    fn unsupported_feature(reason: impl Into<String>) -> Self {
        TermBuilderError::UnsupportedFeature {
            reason: reason.into(),
        }
    }
    #[inline]
    fn degenerate_data(reason: impl Into<String>) -> Self {
        TermBuilderError::DegenerateData {
            reason: reason.into(),
        }
    }
    #[inline]
    fn malformed_formula(reason: impl Into<String>) -> Self {
        TermBuilderError::MalformedFormula {
            reason: reason.into(),
        }
    }
}

// ---------------------------------------------------------------------------
// Column resolution
// ---------------------------------------------------------------------------

/// Resolve a bare column name to its index, returning a typed
/// `DataError::ColumnNotFound` on miss so the FFI boundary can surface a
/// structured `gamfit.ColumnNotFoundError(column=…, available=…)` rather
/// than rely on string-classification of human prose. Internal callers that
/// still flow `Result<_, String>` get byte-identical text via
/// `From<DataError> for String`.
pub fn resolve_col(col_map: &HashMap<String, usize>, name: &str) -> Result<usize, DataError> {
    col_map
        .get(name)
        .copied()
        .ok_or_else(|| DataError::column_not_found(col_map, name, None))
}

/// Like `resolve_col` but tags the missing-column payload with a role label
/// (`"response"`, `"entry"`, `"exit"`, `"event"`, `"z"`, `"id"`, …) so the
/// boundary-side Python exception can disambiguate which formula slot held
/// the bad reference.
pub fn resolve_role_col(
    col_map: &HashMap<String, usize>,
    name: &str,
    role: &str,
) -> Result<usize, DataError> {
    col_map
        .get(name)
        .copied()
        .ok_or_else(|| DataError::column_not_found(col_map, name, Some(role)))
}

fn encoded_levels_for_column(ds: &Dataset, col: ColIdx) -> Vec<(u64, String)> {
    let mut seen = BTreeSet::<u64>::new();
    for value in ds.values.column(col.get()) {
        if value.is_finite() {
            seen.insert(value.to_bits());
        }
    }
    let schema_levels = ds
        .schema
        .columns
        .get(col.get())
        .map(|column| column.levels.as_slice())
        .unwrap_or(&[]);
    seen.into_iter()
        .enumerate()
        .map(|(idx, bits)| {
            let fallback = format!("level{}", idx + 1);
            let label = schema_levels.get(idx).cloned().unwrap_or(fallback);
            (bits, label)
        })
        .collect()
}

pub fn column_map_with_alias(
    col_map: &HashMap<String, usize>,
    alias: &str,
    target_column: &str,
) -> HashMap<String, usize> {
    let mut aliased = col_map.clone();
    if let Some(idx) = col_map.get(target_column).copied() {
        aliased.entry(alias.to_string()).or_insert(idx);
    }
    aliased
}

// ---------------------------------------------------------------------------
// ParsedTerm[] + Dataset → TermCollectionSpec
// ---------------------------------------------------------------------------

pub fn build_termspec(
    terms: &[ParsedTerm],
    ds: &Dataset,
    col_map: &HashMap<String, usize>,
    inference_notes: &mut Vec<String>,
    policy: &ResourcePolicy,
) -> Result<TermCollectionSpec, TermBuilderError> {
    let mut linear_terms = Vec::<LinearTermSpec>::new();
    let mut random_terms = Vec::<RandomEffectTermSpec>::new();
    let mut smooth_terms = Vec::<SmoothTermSpec>::new();
    let smooth_coordinate_count = terms
        .iter()
        .map(|term| match term {
            ParsedTerm::Smooth { vars, .. } => vars.len(),
            _ => 0,
        })
        .sum::<usize>();

    for t in terms {
        match t {
            ParsedTerm::Linear {
                name,
                explicit,
                coefficient_min,
                coefficient_max,
            } => {
                let col = resolve_col(col_map, name)?;
                let auto_kind = ds.column_kinds.get(col).copied().ok_or_else(|| {
                    TermBuilderError::missing_column(format!(
                        "internal column-kind lookup failed for '{name}'"
                    ))
                    .to_string()
                })?;
                if *explicit {
                    linear_terms.push(LinearTermSpec {
                        name: name.clone(),
                        feature_col: col,
                        feature_cols: vec![col],
                        // Parametric linear terms are unpenalized by default
                        // (MLE, matching mgcv/glm); see #749.
                        double_penalty: false,
                        coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
                        coefficient_min: *coefficient_min,
                        coefficient_max: *coefficient_max,
                    });
                } else {
                    match auto_kind {
                        ColumnKindTag::Continuous | ColumnKindTag::Binary => {
                            linear_terms.push(LinearTermSpec {
                                name: name.clone(),
                                feature_col: col,
                                feature_cols: vec![col],
                                // Unpenalized parametric effect by default (#749).
                                double_penalty: false,
                                coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
                                coefficient_min: *coefficient_min,
                                coefficient_max: *coefficient_max,
                            });
                        }
                        ColumnKindTag::Categorical => {
                            if coefficient_min.is_some() || coefficient_max.is_some() {
                                return Err(TermBuilderError::incompatible_config(format!(
                                    "coefficient constraints are not supported for categorical auto-random-effect term '{name}'; use group({name}) or an unconstrained numeric term"
                                )));
                            }
                            random_terms.push(RandomEffectTermSpec {
                                name: name.clone(),
                                feature_col: col,
                                drop_first_level: false,
                                penalized: true,
                                frozen_levels: None,
                            });
                        }
                    }
                }
            }
            ParsedTerm::BoundedLinear {
                name,
                min,
                max,
                prior,
            } => {
                let col = resolve_col(col_map, name)?;
                let auto_kind = ds.column_kinds.get(col).copied().ok_or_else(|| {
                    TermBuilderError::missing_column(format!(
                        "internal column-kind lookup failed for '{name}'"
                    ))
                    .to_string()
                })?;
                if !matches!(auto_kind, ColumnKindTag::Continuous | ColumnKindTag::Binary) {
                    return Err(TermBuilderError::incompatible_config(format!(
                        "bounded() currently supports only numeric columns, got categorical '{name}'"
                    )));
                }
                linear_terms.push(LinearTermSpec {
                    name: name.clone(),
                    feature_col: col,
                    feature_cols: vec![col],
                    double_penalty: false,
                    coefficient_geometry: LinearCoefficientGeometry::Bounded {
                        min: *min,
                        max: *max,
                        prior: prior.clone(),
                    },
                    coefficient_min: None,
                    coefficient_max: None,
                });
            }
            ParsedTerm::RandomEffect { name } => {
                let col = resolve_col(col_map, name)?;
                random_terms.push(RandomEffectTermSpec {
                    name: name.clone(),
                    feature_col: col,
                    drop_first_level: false,
                    penalized: true,
                    frozen_levels: None,
                });
            }
            ParsedTerm::Smooth {
                label,
                vars,
                kind,
                options,
            } => {
                let mut smooth_vars = vars.clone();
                let by_name = options.get("by").cloned();
                let bs_name = options
                    .get("bs")
                    .or_else(|| options.get("type"))
                    .map(|v| v.to_ascii_lowercase());
                let is_sz = matches!(
                    bs_name.as_deref(),
                    Some("sz") | Some("sum-to-zero") | Some("sum_to_zero")
                );
                // For a sum-to-zero factor smooth, identify the categorical
                // factor AMONG the supplied variables rather than assuming a
                // fixed position. mgcv's `s(x, fac, bs="sz")` is
                // continuous-first (the factor acts like a `by=` producing
                // sum-to-zero per-level deviations), exactly like `bs="fs"`,
                // whose builder already detects the factor by column kind. We
                // mirror that here so `s(x, g, bs="sz")` and `s(g, x, ...)` both
                // resolve `g` as the factor — one consistent convention across
                // factor-smooth families.
                let sz_factor_var: Option<String> = if is_sz {
                    if vars.len() < 2 {
                        return Err(format!(
                            "bs=sz smooth '{}' expects a categorical factor and one or more smooth variables",
                            label
                        )
                        .into());
                    }
                    let mut factor_positions = vars.iter().enumerate().filter(|(_, v)| {
                        resolve_col(col_map, v).is_ok_and(|c| {
                            matches!(ds.column_kinds.get(c), Some(ColumnKindTag::Categorical))
                        })
                    });
                    let factor_pos = factor_positions.next().map(|(idx, _)| idx);
                    if factor_positions.next().is_some() {
                        return Err(format!(
                            "bs=sz smooth '{}' expects exactly one categorical factor among its variables; found more than one",
                            label
                        )
                        .into());
                    }
                    let Some(factor_pos) = factor_pos else {
                        return Err(format!(
                            "bs=sz smooth '{}' requires one categorical factor variable; none of its variables is categorical",
                            label
                        )
                        .into());
                    };
                    smooth_vars = vars
                        .iter()
                        .enumerate()
                        .filter(|(idx, _)| *idx != factor_pos)
                        .map(|(_, v)| v.clone())
                        .collect();
                    Some(vars[factor_pos].clone())
                } else {
                    None
                };
                let cols = smooth_vars
                    .iter()
                    .map(|v| resolve_col(col_map, v))
                    .collect::<Result<Vec<_>, _>>()?;
                let mut inner_options = options.clone();
                inner_options.remove("by");
                if is_sz {
                    inner_options.remove("bs");
                    inner_options.remove("type");
                    // mgcv's `bs="sz"` is a SINGLE-penalty smooth: the marginal
                    // wiggliness penalty is replicated across levels and the
                    // marginal's polynomial null space (the per-level linear
                    // trend, once centred) is left UNPENALISED, exactly as the
                    // default thin-plate marginal leaves its null space free.
                    // gam's 1-D smooth otherwise defaults to a double penalty
                    // (an added null-space shrinkage ridge, mgcv `select=TRUE`
                    // semantics). Replicated across the sz deviation blocks that
                    // extra ridge shrinks every level's linear-trend deviation
                    // toward zero — signal that the truth carries (e.g. the
                    // linear projection of sin(2*pi*x) is non-zero) — so REML
                    // over-shrinks the per-level deviations and truth recovery
                    // degrades relative to mgcv's free null space. Force the
                    // inner marginal to a single penalty so the per-level null
                    // space stays free, matching mgcv's sz construction. An
                    // explicit user `double_penalty=` still wins.
                    inner_options
                        .entry("double_penalty".to_string())
                        .or_insert_with(|| "false".to_string());
                }
                // Pop the shape constraint before `build_smooth_basis` runs so
                // it never reaches the per-kind `validate_known_options`
                // allow-lists (the constraint is a property of the smooth term,
                // not of any one basis kind). Basis-incompatible requests still
                // fail loudly downstream via `shape_supports_basis`.
                let shape = match inner_options.remove("shape") {
                    None => ShapeConstraint::None,
                    Some(raw) => crate::terms::smooth::parse_shape_constraint(&raw)
                        .map_err(TermBuilderError::invalid_option)?,
                };
                let inner_basis = build_smooth_basis(
                    *kind,
                    &smooth_vars,
                    &cols,
                    &inner_options,
                    ds,
                    inference_notes,
                    policy,
                    smooth_coordinate_count,
                )?;
                if let Some(factor_var) = sz_factor_var {
                    // `factor_var` was already confirmed categorical when it was
                    // selected above; resolve its column for the level scan.
                    let by_col = resolve_col(col_map, &factor_var)?;
                    let mut levels: Vec<u64> = ds
                        .values
                        .column(by_col)
                        .iter()
                        .map(|v| v.to_bits())
                        .collect();
                    levels.sort_unstable();
                    levels.dedup();
                    smooth_terms.push(SmoothTermSpec {
                        name: label.clone(),
                        basis: SmoothBasisSpec::FactorSumToZero {
                            inner: Box::new(inner_basis),
                            by_col,
                            levels,
                        },
                        shape,
                        joint_null_rotation: None,
                    });
                } else if let Some(by_name) = by_name {
                    let by_col = resolve_col(col_map, &by_name)?;
                    match ds.column_kinds.get(by_col).copied().ok_or_else(|| {
                        format!("internal column-kind lookup failed for by variable '{by_name}'")
                    })? {
                        ColumnKindTag::Categorical => {
                            let levels = encoded_levels_for_column(ds, ColIdx::new(by_col));
                            // Add an unpenalized treatment-coded fixed main effect for the factor, unless present already.
                            if !random_terms
                                .iter()
                                .any(|rt| rt.name == by_name && !rt.penalized)
                            {
                                random_terms.push(RandomEffectTermSpec {
                                    name: by_name.clone(),
                                    feature_col: by_col,
                                    drop_first_level: true,
                                    penalized: false,
                                    frozen_levels: None,
                                });
                            }
                            for (level_bits, level_label) in levels {
                                smooth_terms.push(SmoothTermSpec {
                                    name: format!("{}:{}", label, level_bits),
                                    basis: SmoothBasisSpec::ByVariable {
                                        inner: Box::new(inner_basis.clone()),
                                        by_col,
                                        kind: BySmoothKind::Level { level_bits },
                                        by: ByVariableSpec::Level {
                                            value_bits: level_bits,
                                            label: level_label,
                                        },
                                    },
                                    shape,
                                    joint_null_rotation: None,
                                });
                            }
                        }
                        ColumnKindTag::Binary | ColumnKindTag::Continuous => {
                            smooth_terms.push(SmoothTermSpec {
                                name: label.clone(),
                                basis: SmoothBasisSpec::ByVariable {
                                    inner: Box::new(inner_basis),
                                    by_col,
                                    kind: BySmoothKind::Numeric,
                                    by: ByVariableSpec::Numeric,
                                },
                                shape,
                                joint_null_rotation: None,
                            });
                        }
                    }
                } else {
                    smooth_terms.push(SmoothTermSpec {
                        name: label.clone(),
                        basis: inner_basis,
                        shape,
                        joint_null_rotation: None,
                    });
                }
            }
            ParsedTerm::LinkWiggle { .. }
            | ParsedTerm::TimeWiggle { .. }
            | ParsedTerm::LinkConfig { .. }
            | ParsedTerm::SurvivalConfig { .. } => {
                // Consumed at formula level, not design terms.
            }
            ParsedTerm::LogSlopeSurface { .. } => {
                return Err(TermBuilderError::malformed_formula(
                    "logslope(...) declarations must be resolved by the marginal-slope formula path before building a term spec",
                ));
            }
            ParsedTerm::Interaction { vars } => {
                // Validate every atom resolves to a numeric column. Linear
                // `:` interactions produce a design column equal to the
                // elementwise product of the operands; categorical operands
                // would need a factor-aware expansion (one column per cell)
                // that is not yet wired here.
                let mut cols = Vec::<usize>::with_capacity(vars.len());
                for var in vars {
                    let col = resolve_col(col_map, var)?;
                    let kind = ds.column_kinds.get(col).copied().ok_or_else(|| {
                        TermBuilderError::missing_column(format!(
                            "internal column-kind lookup failed for '{var}'"
                        ))
                        .to_string()
                    })?;
                    match kind {
                        ColumnKindTag::Continuous | ColumnKindTag::Binary => cols.push(col),
                        ColumnKindTag::Categorical => {
                            return Err(TermBuilderError::incompatible_config(format!(
                                "interaction `{}` references categorical column `{var}`; \
                                 linear `:` interactions currently support only numeric \
                                 (continuous/binary) operands. Use te(...) for smooth \
                                 interactions or pre-encode the factor as numeric columns.",
                                vars.join(":")
                            )));
                        }
                    }
                }
                let label = vars.join(":");
                linear_terms.push(LinearTermSpec {
                    name: label,
                    feature_col: cols[0],
                    feature_cols: cols,
                    // Parametric `:` interaction column is unpenalized by
                    // default, same as any other linear term (#749).
                    double_penalty: false,
                    coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
                    coefficient_min: None,
                    coefficient_max: None,
                });
                inference_notes.push(format!(
                    "wired linear interaction `{}` as product of numeric columns",
                    vars.join(":")
                ));
            }
        }
    }

    Ok(TermCollectionSpec {
        linear_terms,
        random_effect_terms: random_terms,
        smooth_terms,
    })
}

fn split_list_option(raw: &str) -> Vec<String> {
    let t = raw.trim();
    let inner = t
        .strip_prefix('[')
        .and_then(|u| u.strip_suffix(']'))
        .unwrap_or(t);
    inner
        .split(',')
        .map(|v| v.trim().to_string())
        .filter(|v| !v.is_empty())
        .collect()
}

fn parse_numeric_expr(raw: &str) -> Result<f64, String> {
    let mut acc = 1.0f64;
    let normalized = raw.replace(' ', "");
    if normalized.eq_ignore_ascii_case("none") {
        return Err("None is not numeric".to_string());
    }
    for factor in normalized.split('*') {
        if factor.is_empty() {
            return Err(format!("invalid numeric expression '{raw}'"));
        }
        let value = if factor.eq_ignore_ascii_case("pi") || factor == "π" {
            std::f64::consts::PI
        } else if factor.eq_ignore_ascii_case("tau") || factor == "τ" {
            std::f64::consts::TAU
        } else if let Some(prefix) = factor
            .strip_suffix("pi")
            .or_else(|| factor.strip_suffix("π"))
        {
            let coefficient = if prefix.is_empty() {
                1.0
            } else {
                prefix
                    .parse::<f64>()
                    .map_err(|err| format!("invalid numeric expression '{raw}': {err}"))?
            };
            coefficient * std::f64::consts::PI
        } else if let Some(prefix) = factor
            .strip_suffix("tau")
            .or_else(|| factor.strip_suffix("τ"))
        {
            let coefficient = if prefix.is_empty() {
                1.0
            } else {
                prefix
                    .parse::<f64>()
                    .map_err(|err| format!("invalid numeric expression '{raw}': {err}"))?
            };
            coefficient * std::f64::consts::TAU
        } else {
            factor
                .parse::<f64>()
                .map_err(|err| format!("invalid numeric expression '{raw}': {err}"))?
        };
        acc *= value;
    }
    Ok(acc)
}

/// Read an endpoint/period option as a numeric *expression* (`2*pi`, `tau`,
/// `0.5*tau`, `6.283185307179586`, ...) — the same grammar that `period=` and
/// `origin=` already accept via [`parse_numeric_expr`].
///
/// Returns `Ok(None)` when the key is absent, `Ok(Some(v))` when it parses, and
/// a hard `Err` when the key is *present but unparseable*. The crucial contrast
/// is with the lenient [`option_f64`], which collapses an unparseable value to
/// `None` and lets the caller silently substitute the data range — wrapping a
/// cyclic smooth at the wrong period with no diagnostic (the #815 failure mode).
fn option_numeric_expr(
    options: &BTreeMap<String, String>,
    key: &str,
) -> Result<Option<f64>, String> {
    match options.get(key) {
        None => Ok(None),
        Some(raw) => parse_numeric_expr(raw)
            .map(Some)
            .map_err(|err| format!("option `{key}={raw}` is not a valid numeric value: {err}")),
    }
}

fn parse_periods_option(
    options: &BTreeMap<String, String>,
    dim: usize,
) -> Result<Option<Vec<Option<f64>>>, String> {
    let Some(raw) = options.get("period") else {
        return Ok(None);
    };
    let values = split_list_option(raw);
    let mut periods = vec![None; dim];
    if values.len() == 1 && dim == 1 {
        periods[0] = Some(parse_numeric_expr(&values[0])?);
    } else {
        if values.len() != dim {
            return Err(format!(
                "period list length {} must match smooth dimension {}",
                values.len(),
                dim
            ));
        }
        for (i, v) in values.iter().enumerate() {
            if v.eq_ignore_ascii_case("none") {
                continue;
            }
            periods[i] = Some(parse_numeric_expr(v)?);
        }
    }
    Ok(Some(periods))
}

fn parse_periodic_axes_option(
    options: &BTreeMap<String, String>,
    dim: usize,
) -> Result<Option<Vec<Option<f64>>>, String> {
    let Some(raw_axes) = options.get("periodic") else {
        return Ok(None);
    };
    let mut periods = parse_periods_option(options, dim)?.unwrap_or_else(|| vec![None; dim]);
    let axes = split_list_option(raw_axes);
    if axes.is_empty() {
        return Ok(Some(periods));
    }
    for a in axes {
        let axis = a
            .parse::<usize>()
            .map_err(|err| format!("invalid periodic axis '{a}': {err}"))?;
        if axis >= dim {
            return Err(format!(
                "periodic axis {axis} out of range for {dim}D smooth"
            ));
        }
        if periods[axis].is_none() {
            return Err(format!(
                "periodic axis {axis} requires period[{axis}] to be finite"
            ));
        }
    }
    // Axes not listed are non-periodic even if period list has a finite placeholder.
    let listed: std::collections::BTreeSet<usize> = split_list_option(raw_axes)
        .into_iter()
        .filter_map(|a| a.parse::<usize>().ok())
        .collect();
    for i in 0..dim {
        if !listed.contains(&i) {
            periods[i] = None;
        }
    }
    Ok(Some(periods))
}

// ---------------------------------------------------------------------------
// Smooth basis spec construction
// ---------------------------------------------------------------------------

fn parse_option_list(raw: &str) -> Vec<String> {
    let trimmed = raw.trim();
    // Accept both the Python/JSON list form `[a, b]` and mgcv's R vector form
    // `c(a, b)` (and a bare `(a, b)`) as the bracketed wrapper around a
    // comma-separated option list. mgcv writes per-margin options as
    // `bs=c('tp','tp')` / `m=c(2,2)`, so the `c(...)` form must round-trip
    // through the same splitter the `[...]` form uses.
    let inner = trimmed
        .strip_prefix('[')
        .and_then(|v| v.strip_suffix(']'))
        .or_else(|| {
            trimmed
                .strip_prefix("c(")
                .or_else(|| trimmed.strip_prefix("C("))
                .or_else(|| trimmed.strip_prefix('('))
                .and_then(|v| v.strip_suffix(')'))
        })
        .unwrap_or(trimmed);
    inner
        .split(',')
        .map(|v| {
            v.trim()
                .trim_matches('"')
                .trim_matches('\'')
                .to_ascii_lowercase()
        })
        .filter(|v| !v.is_empty())
        .collect()
}

fn parse_periodic_axes(
    options: &BTreeMap<String, String>,
    dim: usize,
) -> Result<Vec<bool>, String> {
    let mut axes = vec![false; dim];
    if let Some(raw) = options.get("periodic").or_else(|| options.get("cyclic")) {
        let lowered = raw.trim().to_ascii_lowercase();
        match lowered.as_str() {
            "true" | "yes" | "y" => {
                axes.fill(true);
                return Ok(axes);
            }
            "false" | "no" | "n" => return Ok(axes),
            _ => {}
        }
        for axis_raw in parse_option_list(raw) {
            let axis = axis_raw
                .parse::<usize>()
                .map_err(|err| format!("invalid periodic axis '{axis_raw}': {err}"))?;
            if axis >= dim {
                return Err(format!(
                    "periodic axis {axis} out of range for {dim}D smooth"
                ));
            }
            axes[axis] = true;
        }
    }
    if let Some(raw) = options.get("boundary").or_else(|| options.get("bc")) {
        let boundary = parse_option_list(raw);
        if boundary.len() == dim {
            for (axis, value) in boundary.iter().enumerate() {
                if matches!(value.as_str(), "periodic" | "cyclic" | "cc") {
                    axes[axis] = true;
                }
            }
        } else if dim == 1
            && matches!(
                boundary.first().map(String::as_str),
                Some("periodic" | "cyclic" | "cc")
            )
        {
            axes[0] = true;
        }
    }
    Ok(axes)
}

fn parse_optional_numeric_list(
    options: &BTreeMap<String, String>,
    keys: &[&str],
    dim: usize,
) -> Result<Vec<Option<f64>>, String> {
    let Some(raw) = keys.iter().find_map(|key| options.get(*key)) else {
        return Ok(vec![None; dim]);
    };
    let values = split_list_option(raw);
    let mut out = vec![None; dim];
    if values.len() == 1 && dim == 1 {
        if !values[0].eq_ignore_ascii_case("none") {
            out[0] = Some(parse_numeric_expr(&values[0])?);
        }
        return Ok(out);
    }
    if values.len() != dim {
        return Err(format!(
            "numeric option list length {} must match smooth dimension {}",
            values.len(),
            dim
        ));
    }
    for (i, value) in values.iter().enumerate() {
        if !value.eq_ignore_ascii_case("none") {
            out[i] = Some(parse_numeric_expr(value)?);
        }
    }
    Ok(out)
}

fn parse_periods(
    options: &BTreeMap<String, String>,
    periodic_axes: &[bool],
) -> Result<Vec<Option<f64>>, String> {
    let periods =
        parse_optional_numeric_list(options, &["period", "periods"], periodic_axes.len())?;
    for (axis, (periodic, period)) in periodic_axes.iter().zip(periods.iter()).enumerate() {
        if *periodic
            && let Some(value) = period
            && (!value.is_finite() || *value <= 0.0)
        {
            return Err(format!(
                "period for periodic axis {axis} must be finite and positive, got {value}"
            ));
        }
    }
    Ok(periods)
}

fn parse_period_origins(
    options: &BTreeMap<String, String>,
    periodic_axes: &[bool],
) -> Result<Vec<Option<f64>>, String> {
    parse_optional_numeric_list(
        options,
        &[
            "origin",
            "origins",
            "period_origin",
            "period-origin",
            "domain_origin",
        ],
        periodic_axes.len(),
    )
}

/// Parse a per-axis periodic flag list for tensor smooths. Accepts three forms:
/// - `periodic=true` / `periodic=false` (scalar applied to every axis),
/// - `periodic=[true, false, ...]` (one flag per axis, length `dim`), and
/// - `periodic=[0, 2, ...]` (axis indices that are periodic; others are not).
///
/// `boundary=[..., "periodic"/"cyclic"/"cc", ...]` may also flip individual
/// axes on; non-matching tokens leave the existing flag unchanged.
fn parse_tensor_periodic_axes(
    options: &BTreeMap<String, String>,
    dim: usize,
) -> Result<Vec<bool>, String> {
    let mut axes = vec![false; dim];
    if let Some(raw) = options.get("periodic").or_else(|| options.get("cyclic")) {
        let lowered = raw.trim().to_ascii_lowercase();
        match lowered.as_str() {
            "true" | "yes" | "y" => {
                axes.fill(true);
            }
            "false" | "no" | "n" => {
                // Already false; allow `boundary=` below to flip axes if set.
            }
            _ => {
                let entries = parse_option_list(raw);
                let all_bool = !entries.is_empty()
                    && entries.iter().all(|v| {
                        matches!(
                            v.as_str(),
                            "true" | "yes" | "y" | "false" | "no" | "n" | "none"
                        )
                    });
                if all_bool {
                    if entries.len() != dim {
                        return Err(format!(
                            "periodic list length {} must match smooth dimension {}",
                            entries.len(),
                            dim
                        ));
                    }
                    for (i, v) in entries.iter().enumerate() {
                        axes[i] = matches!(v.as_str(), "true" | "yes" | "y");
                    }
                } else {
                    for axis_raw in entries {
                        let axis = axis_raw
                            .parse::<usize>()
                            .map_err(|err| format!("invalid periodic axis '{axis_raw}': {err}"))?;
                        if axis >= dim {
                            return Err(format!(
                                "periodic axis {axis} out of range for {dim}D smooth"
                            ));
                        }
                        axes[axis] = true;
                    }
                }
            }
        }
    }
    if let Some(raw) = options.get("boundary").or_else(|| options.get("bc")) {
        let boundary = parse_option_list(raw);
        if boundary.len() == dim {
            for (axis, value) in boundary.iter().enumerate() {
                if matches!(value.as_str(), "periodic" | "cyclic" | "cc") {
                    axes[axis] = true;
                }
            }
        }
    }
    Ok(axes)
}

/// Parse a per-margin basis dimension list (`k=<scalar>` or `k=[k0, k1, ...]`).
/// A scalar is broadcast across all axes; `None` returns the heuristic from the
/// data column.
fn parse_tensor_k_list(
    options: &BTreeMap<String, String>,
    cols: &[usize],
    ds: &Dataset,
) -> Result<(Vec<usize>, bool), String> {
    let raw = options
        .get("k")
        .or_else(|| options.get("basis_dim"))
        .or_else(|| options.get("basis-dim"))
        .or_else(|| options.get("basisdim"));
    let Some(raw) = raw else {
        let inferred = cols
            .iter()
            .map(|&c| heuristic_knots_for_column(ds.values.column(c)))
            .collect();
        return Ok((inferred, true));
    };
    let entries = split_list_option(raw);
    if entries.len() == 1 {
        let k: usize = entries[0]
            .parse()
            .map_err(|err| format!("invalid tensor k '{}': {err}", entries[0]))?;
        return Ok((vec![k; cols.len()], false));
    }
    if entries.len() != cols.len() {
        return Err(format!(
            "tensor k list length {} must match smooth dimension {}",
            entries.len(),
            cols.len()
        ));
    }
    let mut out = Vec::with_capacity(entries.len());
    for entry in entries {
        let k: usize = entry
            .parse()
            .map_err(|err| format!("invalid tensor k '{entry}': {err}"))?;
        out.push(k);
    }
    Ok((out, false))
}

/// Parse the `identifiability=` option for tensor-product smooths. Mirrors the
/// vocabulary of the Matern/Duchon parsers so the formula DSL is consistent.
///
/// `kind` selects the default identifiability when no explicit
/// `identifiability=` option is supplied: `te(...)` ([`SmoothKind::Te`]) keeps
/// the full-tensor sum-to-zero default, while `ti(...)` ([`SmoothKind::Ti`])
/// defaults to per-margin sum-to-zero so the marginal main effects are excluded
/// (the mgcv tensor-interaction semantics). An explicit option always wins.
fn parse_tensor_identifiability(
    options: &BTreeMap<String, String>,
    kind: SmoothKind,
) -> Result<TensorBSplineIdentifiability, String> {
    let Some(raw) = options.get("identifiability").map(String::as_str) else {
        return Ok(match kind {
            SmoothKind::Ti => TensorBSplineIdentifiability::MarginalSumToZero,
            _ => TensorBSplineIdentifiability::default(),
        });
    };
    match raw.trim().to_ascii_lowercase().as_str() {
        "none" => Ok(TensorBSplineIdentifiability::None),
        "sum_tozero" | "sum-to-zero" | "center_sum_tozero" | "center-sum-to-zero" | "centered"
        | "sumtozero" => Ok(TensorBSplineIdentifiability::SumToZero),
        "marginal_sum_tozero" | "marginal-sum-to-zero" | "marginal_sumtozero"
        | "marginalsumtozero" | "interaction" => {
            Ok(TensorBSplineIdentifiability::MarginalSumToZero)
        }
        other => Err(TermBuilderError::unsupported_feature(format!(
            "invalid tensor identifiability '{other}'; expected one of: none, sum_tozero, marginal_sum_tozero"
        ))
        .to_string()),
    }
}

fn bspline_boundary_declares_periodic_axis(options: &BTreeMap<String, String>) -> bool {
    options
        .get("boundary")
        .or_else(|| options.get("bc"))
        .map(|raw| {
            parse_option_list(raw)
                .into_iter()
                .any(|value| matches!(value.as_str(), "periodic" | "cyclic" | "cc"))
        })
        .unwrap_or(false)
}

/// Canonical-name lookup for the `bs=`/`type=` smooth selector.
///
/// User-facing names — including mgcv-compatible spellings whose semantics
/// match an existing gamfit smooth exactly — collapse to the engine-internal
/// canonical names used by the dispatch in [`build_smooth_basis`]. Adding a
/// new exactly-equivalent alias is a one-line entry here; the match arms
/// below remain the single dispatch site.
///
/// Aliases listed here MUST be true semantic equivalents of the canonical
/// target, not approximations. mgcv names whose semantics differ from any
/// gamfit smooth (e.g. `bs="cs"` shrinkage thin-plate, `bs="ad"` adaptive)
/// are intentionally NOT mapped here — they should reach the unsupported-type
/// path so users get a real diagnostic instead of a silent semantic
/// substitution.
///
/// Unrecognised inputs pass through unchanged so the dispatch can produce its
/// usual "unsupported smooth type" error, preserving the existing diagnostic
/// surface for genuine typos.
pub(crate) fn canonicalize_smooth_type(raw: &str) -> &str {
    match raw {
        // Thin-plate spline. mgcv `bs="tp"` is the default thin-plate
        // regression spline — exact semantic equivalent of gamfit's `"tps"`.
        "tp" => "tps",
        // Gaussian process / Matérn. mgcv `bs="gp"` defaults to a Matérn
        // covariance kernel with REML smoothing parameter selection, which
        // matches gamfit's `"matern"` exactly (same kernel-Gram identity,
        // same REML route).
        "gp" => "matern",
        other => other,
    }
}

/// Is `margin_bs` a per-margin basis name that the tensor builder realizes as a
/// penalized 1-D B-spline margin?
///
/// gam's tensor product is built from penalized B-spline marginals. mgcv's
/// thin-plate (`tp`/`tps`), P-spline (`ps`), B-spline (`bs`), cubic-regression
/// (`cr`/`cs`), and cyclic (`cc`/`cp`/`cyclic`) marginals are all penalized
/// splines spanning the same per-axis smoothing space, so a B-spline margin
/// reproduces the same tensor smoothing class. Margin kinds with fundamentally
/// different structure (adaptive, random-effect, sphere) are NOT accepted as
/// tensor margins.
pub(crate) fn tensor_margin_bs_is_supported(margin_bs: &str) -> bool {
    matches!(
        canonicalize_smooth_type(margin_bs),
        "tps" | "ps" | "bs" | "bspline" | "cr" | "cs" | "cc" | "cp" | "cyclic"
    )
}

/// Does the smooth request a periodic/cyclic axis via its options?
///
/// Mirrors the boundary-condition reading used by the periodic-aware dispatch
/// branches. Factored out so the type resolver and `build_smooth_basis` agree
/// on a single notion of "periodic requested".
pub(crate) fn smooth_options_declare_periodic(options: &BTreeMap<String, String>) -> bool {
    options.contains_key("periodic")
        || options.contains_key("cyclic")
        || options
            .get("boundary")
            .or_else(|| options.get("bc"))
            .map(|boundary| {
                boundary.to_ascii_lowercase().contains("periodic")
                    || boundary.to_ascii_lowercase().contains("cyclic")
            })
            .unwrap_or(false)
}

/// Resolve the canonical engine-internal smooth-type name for a term.
///
/// Reads the user-facing `type=`/`bs=` selector and collapses mgcv-compatible
/// aliases (`tp`→`tps`, `gp`→`matern`) via [`canonicalize_smooth_type`], or
/// derives the default from the smooth kind/arity when no selector is given.
/// This is the single source of truth for the dispatch in
/// [`build_smooth_basis`]; other call sites (e.g. predictor-specific basis
/// policy) use it so the classification never drifts from the dispatch.
/// Is the raw `bs=`/`type=` selector a vector literal (`c('tp','tp')`,
/// `['tp','tp']`, `(tp, tp)`) rather than a scalar smooth-type name?
///
/// mgcv's tensor smooths take a *per-margin* basis vector
/// (`te(x1, x2, bs=c('tp','tp'))`). Such a value is not a scalar canonical
/// type and must not be fed through [`canonicalize_smooth_type`] — it has to be
/// recognized as a tensor request and split into per-margin types. A scalar
/// selector (`bs="tp"`) is left untouched.
pub(crate) fn bs_selector_is_vector(raw: &str) -> bool {
    let trimmed = raw.trim();
    let bracketed = (trimmed.starts_with('[') && trimmed.ends_with(']'))
        || (trimmed.starts_with("c(") || trimmed.starts_with("C(")) && trimmed.ends_with(')')
        || (trimmed.starts_with('(') && trimmed.ends_with(')'));
    bracketed && !parse_option_list(trimmed).is_empty()
}

pub(crate) fn resolve_smooth_type_name(
    kind: SmoothKind,
    n_cols: usize,
    options: &BTreeMap<String, String>,
) -> String {
    let selector = options.get("type").or_else(|| options.get("bs"));
    // A per-margin basis vector is a tensor request, never a scalar type. Route
    // it to the tensor builder, which reads the per-margin types out of the
    // same `bs=` option. (A vector on a non-tensor smooth is ill-formed and
    // falls through to the scalar path below so the existing diagnostic fires.)
    if let Some(raw) = selector
        && bs_selector_is_vector(raw)
        && matches!(kind, SmoothKind::Te | SmoothKind::Ti)
    {
        return "tensor".to_string();
    }
    selector
        .map(|s| canonicalize_smooth_type(&s.to_ascii_lowercase()).to_string())
        .unwrap_or_else(|| match kind {
            SmoothKind::Te | SmoothKind::Ti => "tensor".to_string(),
            SmoothKind::S if n_cols == 1 => "bspline".to_string(),
            // Mixed periodic Euclidean radial kernels are not separable on the
            // cylinder. Use a tensor product with a cyclic margin so s(theta,h)
            // honors seam continuity while preserving the formula-level s(...).
            SmoothKind::S if smooth_options_declare_periodic(options) => "tensor".to_string(),
            SmoothKind::S => "tps".to_string(),
        })
}

/// Does this canonical smooth type size its basis through the generous spatial
/// center heuristic ([`crate::terms::basis::default_num_centers`])?
///
/// Only the radial spatial bases (thin-plate, Matérn/GP, Duchon) route their
/// default basis dimension through `plan_spatial_basis(.., Default, ..)`. The
/// B-spline, cyclic, tensor, and factor-smooth bases use their own modest
/// knot-based defaults, so they are unaffected by — and must not be perturbed
/// by — secondary-predictor basis-parsimony adjustments (#501).
pub(crate) fn smooth_type_uses_spatial_center_heuristic(canonical_type: &str) -> bool {
    matches!(canonical_type, "tps" | "matern" | "duchon")
}

pub fn build_smooth_basis(
    kind: SmoothKind,
    vars: &[String],
    cols: &[usize],
    options: &BTreeMap<String, String>,
    ds: &Dataset,
    inference_notes: &mut Vec<String>,
    policy: &ResourcePolicy,
    smooth_coordinate_count: usize,
) -> Result<SmoothBasisSpec, String> {
    // Fail fast on degenerate input columns: a smooth over a column that takes
    // only one finite value can only ever fit the response mean — the design
    // matrix is rank-1, and the user almost certainly didn't mean to model a
    // constant predictor as a smooth. Without this guard, `smooth(x)` and
    // `matern(x)` silently fit the mean of `y` regardless of `x`, and the
    // user has no way to tell from looking at the predictions (they're all
    // the same number). Duchon already errors loudly via the basis layer
    // ("smooth basis collapses onto the parametric block"); this lift makes
    // the same diagnosis explicit and uniform across smooth families.
    for (var, &col) in vars.iter().zip(cols.iter()) {
        if matches!(ds.column_kinds.get(col), Some(ColumnKindTag::Categorical)) {
            continue;
        }
        if unique_count_column(ds.values.column(col)) <= 1 {
            return Err(TermBuilderError::degenerate_data(format!(
                "smooth term over '{var}' has only one unique value in the training data \
                 — a smooth on a constant column is degenerate and would only fit the response mean. \
                 Remove `{var}` from the smooth, drop the term, or check the data."
            ))
            .to_string());
        }
    }
    if let Some(by_name) = options.get("by").cloned() {
        let by_col = options
            .get("__by_col")
            .and_then(|raw| raw.parse::<usize>().ok())
            .or_else(|| vars.iter().position(|v| v == &by_name).map(|idx| cols[idx]))
            .ok_or_else(|| format!("unknown by= column '{by_name}'"))?;
        let mut inner_options = options.clone();
        inner_options.remove("by");
        inner_options.remove("__by_col");
        inner_options.remove("id");
        let inner = build_smooth_basis(
            kind,
            vars,
            cols,
            &inner_options,
            ds,
            inference_notes,
            policy,
            smooth_coordinate_count,
        )?;
        let by_kind = match ds.column_kinds.get(by_col).copied() {
            Some(ColumnKindTag::Categorical) => ByVarKind::Factor {
                feature_col: by_col,
                ordered: option_bool(options, "ordered").unwrap_or(false),
                frozen_levels: None,
            },
            Some(ColumnKindTag::Continuous | ColumnKindTag::Binary) => ByVarKind::Numeric {
                feature_col: by_col,
            },
            None => {
                return Err(format!(
                    "internal column-kind lookup failed for by='{by_name}'"
                ));
            }
        };
        return Ok(SmoothBasisSpec::BySmooth {
            smooth: Box::new(inner),
            by_kind,
        });
    }

    let smooth_double_penalty = option_bool(options, "double_penalty").unwrap_or(true);
    let type_opt = resolve_smooth_type_name(kind, cols.len(), options);

    if matches!(type_opt.as_str(), "fs" | "sz" | "re") {
        validate_known_options(
            type_opt.as_str(),
            options,
            &[
                "type",
                "bs",
                "k",
                "basis_dim",
                "basis-dim",
                "basisdim",
                "knots",
                "knot_placement",
                "knot-placement",
                "knotplacement",
                "degree",
                "penalty_order",
                "m",
                "double_penalty",
                "ordered",
            ],
        )?;
        if cols.len() != 2 {
            return Err(format!(
                "{} factor-smooth currently expects exactly two variables (one numeric, one categorical)",
                type_opt
            ));
        }
        let kinds = cols
            .iter()
            .map(|&c| ds.column_kinds.get(c).copied())
            .collect::<Vec<_>>();
        let (cont_idx, group_idx) = if type_opt == "re" {
            // mgcv random-slope examples are often s(g, x, bs="re").
            match (kinds[0], kinds[1]) {
                (Some(ColumnKindTag::Categorical), _) => (1usize, 0usize),
                (_, Some(ColumnKindTag::Categorical)) => (0usize, 1usize),
                _ => (1usize, 0usize),
            }
        } else {
            match (kinds[0], kinds[1]) {
                (_, Some(ColumnKindTag::Categorical)) => (0usize, 1usize),
                (Some(ColumnKindTag::Categorical), _) => (1usize, 0usize),
                _ => {
                    return Err(format!(
                        "{} factor-smooth requires one categorical factor variable",
                        type_opt
                    ));
                }
            }
        };
        let c = cols[cont_idx];
        let (minv, maxv) = col_minmax(ds.values.column(c))?;
        let degree = if type_opt == "re" {
            1
        } else {
            option_usize(options, "degree").unwrap_or(3)
        };
        // For a factor smooth every group's curve is fit from THAT group's rows
        // alone, so the marginal's flexibility must respect the least-resolved
        // group, not the pooled column. The pooled heuristic can hand the marginal
        // a basis that saturates (or exceeds) a small group's sample — e.g. the
        // sleepstudy panel has 8 training days per subject, and a default cubic
        // basis of 8 functions interpolates each subject's 8 points, leaving no
        // room for the wiggliness penalty to collapse the curve toward the
        // per-subject line. The factor smooth then fits within-group noise and
        // extrapolates badly (held-out forecast worse than the population mean).
        //
        // Cap the marginal basis below the minimum per-group covariate resolution
        // so the penalty always retains residual degrees of freedom to shrink each
        // group's curvature toward its linear null space (the random-slope
        // estimand). Groups with ample data (e.g. 40 points each) keep the full
        // pooled flexibility; only small-sample groups are protected. The cap is
        // skipped for the explicit `re` random-effect form, whose degree-1 marginal
        // carries no curvature to over-fit.
        let pooled_internal = heuristic_knots_for_column(ds.values.column(c));
        let default_internal = if type_opt == "re" {
            pooled_internal
        } else {
            let min_group_resolution =
                min_per_group_unique_count(ds.values.column(c), ds.values.column(cols[group_idx]));
            // Per-group basis dim = degree + 1 + internal. Hold it well below the
            // smallest group's resolution (leave at least two residual points per
            // group) so the smooth cannot interpolate that group and the
            // wiggliness penalty retains the room to collapse each curve toward
            // its linear null space. Never drop below `degree + 2`, which keeps
            // exactly the linear span plus a single curvature direction — the
            // minimal smoother that can still bend if the data demand it.
            let basis_cap = min_group_resolution.saturating_sub(2).max(degree + 2);
            let internal_cap = basis_cap.saturating_sub(degree + 1);
            pooled_internal.min(internal_cap.max(1))
        };
        let (n_knots, _) = parse_ps_internal_knots(options, degree, default_internal)?;
        let marginal = BSplineBasisSpec {
            degree,
            penalty_order: option_usize(options, "penalty_order").unwrap_or(if degree > 1 {
                2
            } else {
                1
            }),
            knotspec: resolve_nonperiodic_bspline_knotspec(
                options,
                ds.values.column(c),
                (minv, maxv),
                degree,
                n_knots,
            )?,
            // mgcv's `bs="fs"` is a random-effect-style smooth: EVERY per-level
            // coefficient, including the marginal null space, is penalized so
            // unobserved groups can be predicted — so `fs` keeps the null-space
            // (double) penalty. mgcv's `bs="sz"` is a pure across-level
            // *deviation* smooth that, under the default `select=FALSE`, leaves
            // the per-level null space UNPENALIZED; carrying the double penalty
            // there shrinks the genuine deviation signal and over-smooths the
            // recovered curves relative to mgcv (gam#700). `re` carries its own
            // identity ridge below and ignores this flag. Honour an explicit
            // user `double_penalty=` either way.
            double_penalty: option_bool(options, "double_penalty")
                .unwrap_or(type_opt.as_str() != "sz"),
            identifiability: BSplineIdentifiability::None,
            boundary_conditions: Default::default(),
            boundary: OneDimensionalBoundary::Open,
        };
        let flavour = match type_opt.as_str() {
            "fs" => FactorSmoothFlavour::Fs {
                m_null_penalty_orders: vec![option_usize(options, "m").unwrap_or(2)],
            },
            "sz" => FactorSmoothFlavour::Sz,
            "re" => FactorSmoothFlavour::Re,
            // Outer `matches!` already restricts to fs/sz/re.
            other => {
                return Err(format!(
                    "internal: factor-smooth flavour dispatch reached unexpected type `{}`",
                    other
                ));
            }
        };
        return Ok(SmoothBasisSpec::FactorSmooth {
            spec: FactorSmoothSpec {
                continuous_cols: vec![c],
                group_col: cols[group_idx],
                marginal,
                flavour,
                group_frozen_levels: None,
            },
        });
    }

    match type_opt.as_str() {
        "cyclic" | "cc" | "cp" | "cyclic-ps" => {
            validate_known_options(
                "cyclic",
                options,
                &[
                    "type",
                    "bs",
                    "by",
                    "k",
                    "basis_dim",
                    "basis-dim",
                    "basisdim",
                    "degree",
                    "penalty_order",
                    "period",
                    "periods",
                    "period_start",
                    "period_end",
                    "start",
                    "end",
                    "origin",
                    "origins",
                    "period_origin",
                    "period-origin",
                    "domain_origin",
                    "double_penalty",
                    "id",
                    "__by_col",
                    "identifiability",
                ],
            )?;
            if cols.len() != 1 {
                return Err(format!(
                    "periodic smooth expects one variable, got {}",
                    cols.len()
                ));
            }
            let c = cols[0];
            let (minv, maxv) = col_minmax(ds.values.column(c))?;
            let degree = option_usize(options, "degree").unwrap_or(3);
            let mut default_internal = heuristic_knots_for_column(ds.values.column(c));
            if ds.values.nrows() <= 32 && smooth_coordinate_count >= 5 {
                default_internal = default_internal.min(1);
            }
            let default_basis = default_internal + degree + 1;
            let num_basis = option_usize_any(options, &["k", "basis_dim", "basis-dim", "basisdim"])
                .unwrap_or(default_basis);
            if num_basis < degree + 1 {
                return Err(format!(
                    "periodic smooth: k={} too small for degree {}; expected k >= {}",
                    num_basis,
                    degree,
                    degree + 1
                ));
            }
            // The cyclic arm is periodic on its single axis by construction, so
            // resolve the period exactly the way the `s()`/`ps` arm does: honour
            // `period=`/`periods=` first (with `origin=` setting the domain
            // start), and fall back to the `period_start`/`period_end` endpoint
            // form only when `period=` is absent. Previously this arm jumped
            // straight to `parse_periodic_domain_1d`, so a `period=<v>`
            // declaration was silently dropped and the smooth wrapped at the
            // data range (#816). All three helpers route through
            // `parse_numeric_expr`, so `period=2*pi` and `period_end=2*pi` parse
            // identically (#815).
            let periodic_axes = [true];
            let periods = parse_periods(options, &periodic_axes)?;
            let origins = parse_period_origins(options, &periodic_axes)?;
            let (domain_start, period) = if let Some(p) = periods[0] {
                (origins[0].unwrap_or(minv), p)
            } else {
                parse_periodic_domain_1d(options, minv, maxv)?
            };
            Ok(SmoothBasisSpec::BSpline1D {
                feature_col: c,
                spec: BSplineBasisSpec {
                    degree,
                    penalty_order: option_usize(options, "penalty_order").unwrap_or(2),
                    knotspec: BSplineKnotSpec::PeriodicUniform {
                        data_range: (domain_start, domain_start + period),
                        num_basis,
                    },
                    double_penalty: smooth_double_penalty,
                    identifiability: BSplineIdentifiability::default(),
                    boundary_conditions: Default::default(),
                    boundary: OneDimensionalBoundary::Cyclic {
                        start: domain_start,
                        end: domain_start + period,
                    },
                },
            })
        }
        "bspline" | "ps" | "p-spline" => {
            validate_known_options(
                "bspline",
                options,
                &[
                    "type",
                    "bs",
                    "by",
                    "k",
                    "basis_dim",
                    "basis-dim",
                    "basisdim",
                    "knots",
                    "knot_placement",
                    "knot-placement",
                    "knotplacement",
                    "degree",
                    "penalty_order",
                    "boundary",
                    "bc",
                    "boundary_conditions",
                    "bc_left",
                    "bc_right",
                    "left_bc",
                    "right_bc",
                    "start_bc",
                    "end_bc",
                    "side",
                    "anchor",
                    "anchor_value",
                    "value",
                    "anchor_left",
                    "left_anchor",
                    "anchor_right",
                    "right_anchor",
                    "periodic",
                    "period",
                    "periods",
                    "period_start",
                    "period_end",
                    "origin",
                    "double_penalty",
                    "by",
                    "id",
                    "__by_col",
                    "identifiability",
                    "by",
                ],
            )?;
            if cols.len() != 1 {
                return Err(TermBuilderError::incompatible_config(format!(
                    "bspline smooth expects one variable, got {}",
                    cols.len()
                ))
                .to_string());
            }
            let c = cols[0];
            let (minv, maxv) = col_minmax(ds.values.column(c))?;
            let degree = option_usize(options, "degree").unwrap_or(3);
            let default_internal = heuristic_knots_for_column(ds.values.column(c));
            let (mut n_knots, inferred) =
                parse_ps_internal_knots(options, degree, default_internal)?;
            if inferred && ds.values.nrows() <= 32 && smooth_coordinate_count >= 5 {
                n_knots = n_knots.min(1);
            }
            if inferred {
                let unique = unique_count_column(ds.values.column(c));
                let ceiling = ((unique as f64).cbrt() as usize).max(20);
                inference_notes.push(format!(
                    "Automatically set {} internal knots for smooth '{}' from {} unique values (rule: clamp(unique/4, 4..max(20, cbrt(unique))) = clamp(unique/4, 4..{})). Override with knots=... or k=....",
                    n_knots,
                    vars.join(","),
                    unique,
                    ceiling,
                ));
            }
            let periodic_axes = parse_periodic_axes(options, 1).map_err(|e| e.to_string())?;
            let boundary_conditions =
                if periodic_axes[0] && bspline_boundary_declares_periodic_axis(options) {
                    BSplineBoundaryConditions::default()
                } else {
                    parse_bspline_boundary_conditions(options).map_err(|e| e.to_string())?
                };
            let periods = parse_periods(options, &periodic_axes).map_err(|e| e.to_string())?;
            let origins =
                parse_period_origins(options, &periodic_axes).map_err(|e| e.to_string())?;
            let (knotspec, boundary) = if periodic_axes[0] {
                if !boundary_conditions.is_free() {
                    return Err(TermBuilderError::incompatible_config(
                        "periodic B-splines cannot also declare endpoint boundary conditions",
                    )
                    .to_string());
                }
                {
                    let (domain_start, p_value) = if periods[0].is_some() {
                        (origins[0].unwrap_or(minv), periods[0].unwrap())
                    } else {
                        parse_periodic_domain_1d(options, minv, maxv).map_err(|e| e.to_string())?
                    };
                    let domain_end = domain_start + p_value;
                    (
                        BSplineKnotSpec::PeriodicUniform {
                            data_range: (domain_start, domain_end),
                            num_basis: n_knots + degree + 1,
                        },
                        OneDimensionalBoundary::Cyclic {
                            start: domain_start,
                            end: domain_end,
                        },
                    )
                }
            } else {
                (
                    resolve_nonperiodic_bspline_knotspec(
                        options,
                        ds.values.column(c),
                        (minv, maxv),
                        degree,
                        n_knots,
                    )?,
                    parse_cyclic_boundary(options, minv, maxv)?,
                )
            };
            Ok(SmoothBasisSpec::BSpline1D {
                feature_col: c,
                spec: BSplineBasisSpec {
                    degree,
                    penalty_order: option_usize(options, "penalty_order").unwrap_or(2),
                    knotspec,
                    double_penalty: smooth_double_penalty,
                    identifiability: BSplineIdentifiability::default(),
                    boundary,
                    boundary_conditions,
                },
            })
        }
        "tps" | "thinplate" | "thin-plate" => {
            validate_known_options(
                "thinplate",
                options,
                &[
                    SECONDARY_CENTER_CAP_OPTION,
                    "type",
                    "bs",
                    "by",
                    "length_scale",
                    "centers",
                    "k",
                    "basis_dim",
                    "basis-dim",
                    "basisdim",
                    "knots",
                    "include_intercept",
                    "double_penalty",
                    "by",
                    "id",
                    "__by_col",
                    "identifiability",
                    "by",
                    "scale_dims",
                ],
            )?;
            let plan = plan_spatial_basis(
                ds.values.nrows(),
                cols.len(),
                CenterCountRequest::Default,
                DuchonNullspaceOrder::Linear,
                option_bool(options, "scale_dims").unwrap_or(false),
                policy,
            )
            .map_err(|e| e.to_string())?;
            let centers = parse_countwith_basis_alias(
                options,
                "centers",
                cap_default_spatial_centers(options, plan.centers),
            )?;
            let center_strategy = if has_explicit_countwith_basis_alias(options, "centers") {
                spatial_center_strategy_for_dimension(centers, cols.len())
            } else {
                auto_spatial_center_strategy(centers, cols.len())
            };
            Ok(SmoothBasisSpec::ThinPlate {
                feature_cols: cols.to_vec(),
                spec: ThinPlateBasisSpec {
                    center_strategy,
                    periodic: parse_periodic_axes_option(options, cols.len())?,
                    // Sentinel: leave at 0.0 when the user didn't pass an
                    // explicit length_scale so `auto_init_length_scale_in_place`
                    // can replace it with a data-derived initialization. The
                    // old hard-coded 1.0 was the documented basin (see
                    // smooth.rs `auto_init_length_scale_in_place`) that the
                    // spatial optimizer could not escape, leaving TPS terms
                    // initialized off the data scale.
                    length_scale: option_f64(options, "length_scale").unwrap_or(0.0),
                    double_penalty: smooth_double_penalty,
                    identifiability: parse_spatial_identifiability(options)
                        .map_err(|e| e.to_string())?,
                    radial_reparam: None,
                },
                input_scales: None,
            })
        }
        "sphere" | "s2" | "sos" => {
            validate_known_options(
                "sphere",
                options,
                &[
                    "type",
                    "bs",
                    "by",
                    "centers",
                    "k",
                    "basis_dim",
                    "basis-dim",
                    "basisdim",
                    "knots",
                    "penalty_order",
                    "m",
                    "double_penalty",
                    "id",
                    "__by_col",
                    "kernel",
                    "method",
                    "radians",
                    "units",
                    "degree",
                    "l",
                    "max_degree",
                    "max-degree",
                ],
            )?;
            if cols.len() != 2 {
                return Err(format!(
                    "sphere smooth expects exactly two variables (lat, lon), got {}",
                    cols.len()
                ));
            }
            let radians = option_bool(options, "radians").unwrap_or_else(|| {
                options
                    .get("units")
                    .map(|u| u.eq_ignore_ascii_case("radian") || u.eq_ignore_ascii_case("radians"))
                    .unwrap_or(false)
            });
            let kernel = options
                .get("kernel")
                .or_else(|| options.get("method"))
                .map(|raw| strip_quotes(raw).trim().to_ascii_lowercase())
                .unwrap_or_else(|| "sobolev".to_string());
            let (method, wahba_kernel) = match kernel.as_str() {
                "sobolev" | "wahba" => (SphereMethod::Wahba, SphereWahbaKernel::Sobolev),
                "pseudo" | "mgcv" | "sos" => (SphereMethod::Wahba, SphereWahbaKernel::Pseudo),
                "harmonic" | "spherical_harmonic" | "spherical-harmonic" => {
                    (SphereMethod::Harmonic, SphereWahbaKernel::Sobolev)
                }
                other => {
                    return Err(format!(
                        "unsupported sphere kernel '{other}'; expected sobolev, pseudo, or harmonic"
                    ));
                }
            };
            let max_degree = if matches!(method, SphereMethod::Harmonic) {
                let degree =
                    option_usize_any(options, &["degree", "l", "max_degree", "max-degree"])
                        .or_else(|| option_usize(options, "centers"))
                        .or_else(|| {
                            option_usize_any(options, &["k", "basis_dim", "basis-dim", "basisdim"])
                                .and_then(|k| (1..=128).find(|&l| l * (l + 2) >= k))
                        })
                        .unwrap_or_else(|| default_spherical_harmonic_degree(ds.values.nrows()));
                if degree == 0 {
                    return Err("sphere smooth requires degree/max_degree >= 1".to_string());
                }
                if degree > 32 {
                    return Err(format!(
                        "sphere smooth max_degree={} is too large for the dense harmonic engine (limit 32)",
                        degree
                    ));
                }
                Some(degree)
            } else {
                None
            };
            let center_strategy = if matches!(method, SphereMethod::Wahba) {
                let centers = parse_countwith_basis_alias(
                    options,
                    "centers",
                    default_num_centers(ds.values.nrows(), cols.len()),
                )?;
                CenterStrategy::FarthestPoint {
                    num_centers: centers,
                }
            } else {
                CenterStrategy::FarthestPoint { num_centers: 0 }
            };
            Ok(SmoothBasisSpec::Sphere {
                feature_cols: cols.to_vec(),
                spec: SphericalSplineBasisSpec {
                    center_strategy,
                    penalty_order: option_usize(options, "penalty_order")
                        .or_else(|| option_usize(options, "m"))
                        .unwrap_or(2),
                    double_penalty: smooth_double_penalty,
                    radians,
                    method,
                    max_degree,
                    wahba_kernel,
                    identifiability: SphericalSplineIdentifiability::CenterSumToZero,
                },
            })
        }
        "matern" => {
            // Catch typos like `lengt_scale=` / `nyu=` / `centerz=` before
            // they get silently ignored and the user wonders why their
            // option had no effect. The matern() term accepts exactly
            // these options.
            validate_known_options(
                "matern",
                options,
                &[
                    SECONDARY_CENTER_CAP_OPTION,
                    "type",
                    "bs",
                    "by",
                    "nu",
                    "length_scale",
                    "centers",
                    "k",
                    "basis_dim",
                    "basis-dim",
                    "basisdim",
                    "knots",
                    "include_intercept",
                    "double_penalty",
                    "by",
                    "id",
                    "__by_col",
                    "identifiability",
                    "by",
                    "scale_dims",
                ],
            )?;
            let plan = plan_spatial_basis(
                ds.values.nrows(),
                cols.len(),
                CenterCountRequest::Default,
                DuchonNullspaceOrder::Zero,
                option_bool(options, "scale_dims").unwrap_or(false),
                policy,
            )
            .map_err(|e| e.to_string())?;
            let centers = parse_countwith_basis_alias(
                options,
                "centers",
                cap_default_spatial_centers(options, plan.centers),
            )?;
            let center_strategy = if has_explicit_countwith_basis_alias(options, "centers") {
                spatial_center_strategy_for_dimension(centers, cols.len())
            } else {
                auto_spatial_center_strategy(centers, cols.len())
            };
            let nu = parse_matern_nu(options.get("nu").map(String::as_str).unwrap_or("5/2"))?;
            // The exponential (ν = 1/2) Matérn kernel has a singular Laplacian
            // at zero in d ≥ 2, so the operator-collocation penalty machinery
            // hits a non-invertible matrix during fit. Surface the cause
            // up-front instead of letting the user see the generic
            // "Matrix conditioning issue detected" wrapper from PIRLS.
            if matches!(nu, MaternNu::Half) && cols.len() >= 2 {
                return Err(TermBuilderError::unsupported_feature(format!(
                    "matern() with nu=1/2 is not supported for d>=2 (got {} covariates): \
                     the exponential kernel's Laplacian is singular at center collisions, \
                     which makes the operator-collocation penalty non-invertible. \
                     Choose nu>=3/2 (e.g. nu=3/2 or the default nu=5/2) for multi-dimensional smooths.",
                    cols.len()
                ))
                .to_string());
            }
            let aniso_log_scales = if option_bool(options, "scale_dims").unwrap_or(false) {
                Some(vec![0.0; cols.len()])
            } else {
                None
            };
            Ok(SmoothBasisSpec::Matern {
                feature_cols: cols.to_vec(),
                spec: MaternBasisSpec {
                    center_strategy,
                    periodic: parse_periodic_axes_option(options, cols.len())?,
                    length_scale: option_f64(options, "length_scale")
                        .unwrap_or_else(|| default_matern_length_scale(ds, cols)),
                    nu,
                    include_intercept: option_bool(options, "include_intercept").unwrap_or(false),
                    // Matérn is a STRICTLY POSITIVE-DEFINITE kernel: its penalty
                    // (the kernel Gram `K_CC`) is full-rank with NO structural
                    // polynomial nullspace, so the double-penalty nullspace ridge
                    // (mgcv `select=TRUE`, which shrinks a smooth's UNPENALIZED
                    // trend directions) has nothing to shrink. The default `true`
                    // (`smooth_double_penalty`) therefore emitted a SPURIOUS
                    // `DoublePenaltyNullspace` whose presence flipped with the
                    // κ-relative spectral rank test (`build_nullspace_shrinkage_penalty`):
                    // at a κ where the projected `K_CC` had a near-zero eigenvalue
                    // (a conditioning artifact, not a real null direction) it
                    // appeared, else not — flipping the learned-penalty count 6↔7
                    // across the kappa optimizer's design rebuilds and breaking the
                    // joint-setup rho contract (the BMS spatial setup counts only
                    // genuinely-learned penalties; see the
                    // `spatial_joint_setup_counts_only_learned_penalties_in_rho`
                    // invariant), which raised "joint hyper rho dimension mismatch"
                    // (gam#787). Default OFF for the strictly-PD Matérn kernel; an
                    // explicit user `double_penalty=` still wins. (Duchon/thin-plate
                    // is handled in its own arm and KEEPS its native nullspace
                    // shrinkage — it genuinely has a polynomial nullspace, #754.)
                    double_penalty: option_bool(options, "double_penalty").unwrap_or(false),
                    identifiability: parse_matern_identifiability(options)
                        .map_err(|e| e.to_string())?,
                    aniso_log_scales,
                },
                input_scales: None,
            })
        }
        "duchon" => {
            validate_known_options(
                "duchon",
                options,
                &[
                    SECONDARY_CENTER_CAP_OPTION,
                    "type",
                    "bs",
                    "by",
                    "length_scale",
                    "centers",
                    "k",
                    "basis_dim",
                    "basis-dim",
                    "basisdim",
                    "knots",
                    "power",
                    "p",
                    "nullspace_order",
                    "order",
                    "identifiability",
                    "by",
                    "periodic",
                    "cyclic",
                    "period",
                    "period_start",
                    "period_end",
                    "scale_dims",
                    "double_penalty",
                    "by",
                    "id",
                    "__by_col",
                ],
            )?;
            if options.contains_key("double_penalty") {
                return Err(TermBuilderError::incompatible_config(format!(
                    "Duchon smooth '{}' does not support double_penalty; the Duchon smoother already ships its native reproducing-norm penalty plus a null-space shrinkage ridge.",
                    vars.join(", ")
                ))
                .to_string());
            }
            let requested_nullspace_order = parse_duchon_order(options)?;
            let length_scale = option_f64_strict(options, "length_scale")?;
            // Resolve `(nullspace_order, power)`. The default (magic) path is a
            // structural amplitude/slope/curvature smoother: an affine (`Linear`)
            // polynomial nullspace and spectral power `s = (d - 1)/2`, giving the
            // cubic kernel `r^3` in 1D. There is no nullspace-order escalation —
            // the structural cubic smoother is well-defined for every dimension.
            //
            // Explicit `power=...` honors the user's value verbatim against their
            // requested nullspace order; the kernel validator emits a precise
            // diagnostic for any inadmissible combination. In the scale-free
            // (non-hybrid) regime fractional powers are admitted and threaded as
            // `f64`. The hybrid Duchon-Matérn kernel (`length_scale=Some`) is
            // restricted to integer powers.
            let (nullspace_order, power) = match parse_duchon_power_policy(options)? {
                DuchonPowerPolicy::Explicit(req_power) => {
                    if length_scale.is_some() && req_power.fract() != 0.0 {
                        return Err(TermBuilderError::incompatible_config(format!(
                            "hybrid Duchon-Matern smooth '{}' (length_scale=...) requires an integer power, got power={}; \
                             drop length_scale to use the scale-free structural kernel with a fractional power.",
                            vars.join(", "),
                            req_power,
                        ))
                        .to_string());
                    }
                    (requested_nullspace_order, req_power)
                }
                DuchonPowerPolicy::CubicStructuralDefault => {
                    // Magic cubic rule (REQUEST-LAYER default): no explicit power ⇒
                    // affine null space + fractional spectral power s = (d-1)/2, i.e.
                    // the Duchon kernel φ(r)=r³ in every dimension. An EXPLICIT
                    // `power=0` is handled above and is honored as the s=0 Duchon
                    // kernel (r²·log r ≡ the thin-plate kernel in even d) — the magic
                    // default lives here, not in the basis builder.
                    match length_scale {
                        None => crate::basis::duchon_cubic_default(cols.len()),
                        Some(_) => {
                            // The hybrid Matérn-blended kernel (`length_scale=Some`)
                            // requires an INTEGER spectral power `s` (the partial-
                            // fraction split `1/(ρ^{2p}(κ²+ρ²)^s)` is only defined for
                            // integer `s`). The fractional cubic default `s=(d-1)/2` is
                            // a half-integer for even `d`, and the basis builder's
                            // `power_as_usize` maps a NON-integer to `0` (not its
                            // floor) — so for even `d ≥ 4` the realized kernel has
                            // `2(p+s) = 2p = 4 ≤ d`, which is non-finite at the origin
                            // and crashes the fit (historically a non-finite
                            // eigendecomposition; now a fit-time validation error).
                            //
                            // Rather than emit the fractional cubic and let it truncate
                            // into an inadmissible kernel, resolve the SMALLEST
                            // admissible integer `(nullspace, s)` at the requested
                            // nullspace order, honoring the collocation order of the
                            // default operator penalties (mass + tension ⇒ D1). This
                            // recovers the canonical thin-plate smoothness order
                            // `m = p + s = ⌊d/2⌋ + 1` for the hybrid kernel and agrees
                            // with the fractional cubic default for odd `d` (where the
                            // collocation floor already forces `s = (d-1)/2`).
                            let max_op = crate::basis::duchon_max_active_operator_derivative_order(
                                &DuchonOperatorPenaltySpec::default(),
                            );
                            let (ns, s) = crate::basis::resolve_duchon_orders(
                                cols.len(),
                                requested_nullspace_order,
                                max_op,
                                length_scale,
                            );
                            (ns, s as f64)
                        }
                    }
                }
            };
            let plan = plan_spatial_basis(
                ds.values.nrows(),
                cols.len(),
                CenterCountRequest::Default,
                nullspace_order,
                option_bool(options, "scale_dims").unwrap_or(false),
                policy,
            )
            .map_err(|e| e.to_string())?;
            let centers_explicit = has_explicit_countwith_basis_alias(options, "centers");
            let requested_centers = parse_countwith_basis_alias(
                options,
                "centers",
                cap_default_spatial_centers(options, plan.centers),
            )?;
            let polynomial_cols = match nullspace_order {
                DuchonNullspaceOrder::Zero => 1,
                DuchonNullspaceOrder::Linear => cols.len() + 1,
                DuchonNullspaceOrder::Degree(degree) => {
                    crate::basis::duchon_nullspace_dimension(cols.len(), degree)
                }
            };
            if requested_centers <= polynomial_cols {
                return Err(TermBuilderError::incompatible_config(format!(
                    "Duchon smooth '{}' requested basis dimension {} but order={:?} in {}D needs {} polynomial null-space columns; choose centers/k > {}",
                    vars.join(", "),
                    requested_centers,
                    nullspace_order,
                    cols.len(),
                    polynomial_cols,
                    polynomial_cols,
                ))
                .to_string());
            }
            let mut centers = requested_centers;
            if !centers_explicit && ds.values.nrows() <= 32 && smooth_coordinate_count >= 5 {
                centers = centers.max(polynomial_cols + 4);
            }
            let center_strategy = if centers_explicit {
                spatial_center_strategy_for_dimension(centers, cols.len())
            } else {
                auto_spatial_center_strategy(centers, cols.len())
            };
            let aniso_log_scales = if option_bool(options, "scale_dims").unwrap_or(false) {
                Some(vec![0.0; cols.len()])
            } else {
                None
            };
            // The default is the full Hilbert scale (curvature `Primary` + trend
            // ridge + mass + tension); REML deselects what the data don't support.
            let operator_penalties = DuchonOperatorPenaltySpec::default();
            Ok(SmoothBasisSpec::Duchon {
                feature_cols: cols.to_vec(),
                spec: DuchonBasisSpec {
                    center_strategy,
                    periodic: parse_periodic_axes_option(options, cols.len())?,
                    length_scale,
                    power,
                    nullspace_order,
                    identifiability: parse_spatial_identifiability(options)
                        .map_err(|e| e.to_string())?,
                    aniso_log_scales,
                    operator_penalties,
                    boundary: if cols.len() == 1 {
                        let c = cols[0];
                        let (minv, maxv) = col_minmax(ds.values.column(c))?;
                        parse_cyclic_boundary(options, minv, maxv)?
                    } else {
                        OneDimensionalBoundary::Open
                    },
                },
                input_scales: None,
            })
        }
        "tensor" | "te" | "ti" | "t2" => {
            validate_known_options(
                "tensor",
                options,
                &[
                    "type",
                    "bs",
                    "by",
                    "k",
                    "basis_dim",
                    "basis-dim",
                    "basisdim",
                    "knot_placement",
                    "knot-placement",
                    "knotplacement",
                    "degree",
                    "penalty_order",
                    "double_penalty",
                    "periodic",
                    "cyclic",
                    "period",
                    "periods",
                    "period_start",
                    "period_end",
                    "origin",
                    "origins",
                    "period_origin",
                    "period-origin",
                    "domain_origin",
                    "boundary",
                    "bc",
                    "identifiability",
                    "id",
                    "__by_col",
                ],
            )?;
            if cols.len() < 2 {
                return Err(TermBuilderError::incompatible_config(format!(
                    "tensor smooth expects at least 2 variables, got {}",
                    cols.len()
                ))
                .to_string());
            }
            let dim = cols.len();
            // Per-margin basis vector (`bs=c('tp','tp')` / `bs=['ps','cr']`):
            // validate each requested margin is a penalized-spline basis that
            // the tensor product realizes as a 1-D B-spline margin. mgcv's
            // `tp`/`ps`/`cr`/`bs`/`cc` margins are all penalized splines over
            // the same per-axis function space, so a B-spline margin recovers
            // the same tensor smoothing space; genuinely different margin kinds
            // (e.g. adaptive `ad`, random `re`) are rejected loudly rather than
            // silently substituted.
            if let Some(raw) = options.get("bs").or_else(|| options.get("type"))
                && bs_selector_is_vector(raw)
            {
                let per_margin = parse_option_list(raw);
                if per_margin.len() != dim {
                    return Err(TermBuilderError::invalid_option(format!(
                        "tensor smooth per-margin bs vector has {} entries but the smooth has {} margins",
                        per_margin.len(),
                        dim
                    ))
                    .to_string());
                }
                for (axis, margin_bs) in per_margin.iter().enumerate() {
                    if !tensor_margin_bs_is_supported(margin_bs) {
                        return Err(TermBuilderError::unsupported_feature(format!(
                            "tensor smooth margin {axis} basis '{margin_bs}' is not a supported penalized-spline margin; \
                             tensor margins accept tp/tps/ps/bs/cr/cc"
                        ))
                        .to_string());
                    }
                }
            }
            let periodic_axes = parse_tensor_periodic_axes(options, dim)?;
            let periods_opt = parse_periods(options, &periodic_axes)?;
            let origins_opt = parse_period_origins(options, &periodic_axes)?;
            let degree = option_usize(options, "degree").unwrap_or(3);
            let penalty_order =
                option_usize(options, "penalty_order").unwrap_or(if degree > 1 { 2 } else { 1 });
            let (mut k_list, k_inferred) = parse_tensor_k_list(options, cols, ds)?;
            if ds.values.nrows() <= 32 && smooth_coordinate_count >= 5 {
                for k in &mut k_list {
                    *k = (*k).min(degree + 2);
                }
            }
            if k_inferred {
                inference_notes.push(format!(
                    "Automatically set per-margin basis sizes {:?} for tensor smooth '{}' \
                     (unique/4 rule per column, clamped to [4, cbrt(unique).max(20)]). \
                     Override with k=<int> or k=[k0,k1,...].",
                    k_list,
                    vars.join(",")
                ));
            }
            let mut margins: Vec<BSplineBasisSpec> = Vec::with_capacity(dim);
            let mut emitted_periods: Vec<Option<f64>> = Vec::with_capacity(dim);
            for axis in 0..dim {
                let c = cols[axis];
                let (data_min, data_max) = col_minmax(ds.values.column(c))?;
                let k_axis = k_list[axis];
                let min_k = degree + 1;
                if k_axis < min_k {
                    return Err(TermBuilderError::invalid_option(format!(
                        "tensor smooth: k[{axis}]={k_axis} too small for degree {degree}; expected k >= {min_k}"
                    ))
                    .to_string());
                }
                let (knotspec, boundary, axis_period) = if periodic_axes[axis] {
                    let period_value = periods_opt[axis].ok_or_else(|| {
                        format!(
                            "tensor smooth axis {axis} is periodic but no period was supplied; \
                             pass period=<value> (scalar) or period=[..., <value>, ...]"
                        )
                    })?;
                    if !period_value.is_finite() || period_value <= 0.0 {
                        return Err(format!(
                            "tensor smooth axis {axis}: period must be a positive finite value, got {period_value}"
                        ));
                    }
                    let domain_start = origins_opt[axis].unwrap_or(data_min);
                    let domain_end = domain_start + period_value;
                    (
                        BSplineKnotSpec::PeriodicUniform {
                            data_range: (domain_start, domain_end),
                            num_basis: k_axis,
                        },
                        OneDimensionalBoundary::Cyclic {
                            start: domain_start,
                            end: domain_end,
                        },
                        Some(period_value),
                    )
                } else {
                    let num_internal_knots = k_axis.saturating_sub(degree + 1).max(1);
                    let knotspec = match parse_knot_placement(options)? {
                        crate::basis::BSplineKnotPlacement::Uniform => BSplineKnotSpec::Generate {
                            data_range: (data_min, data_max),
                            num_internal_knots,
                        },
                        crate::basis::BSplineKnotPlacement::Quantile => {
                            crate::basis::auto_knot_vector_1d_quantile(
                                ds.values.column(c),
                                num_internal_knots,
                                degree,
                            )
                            .map_err(|e| e.to_string())?;
                            BSplineKnotSpec::Automatic {
                                num_internal_knots: Some(num_internal_knots),
                                placement: crate::basis::BSplineKnotPlacement::Quantile,
                            }
                        }
                    };
                    (knotspec, OneDimensionalBoundary::Open, None)
                };
                margins.push(BSplineBasisSpec {
                    degree,
                    penalty_order,
                    knotspec,
                    double_penalty: false,
                    identifiability: BSplineIdentifiability::None,
                    boundary,
                    boundary_conditions: BSplineBoundaryConditions::default(),
                });
                emitted_periods.push(axis_period);
            }
            let any_periodic = emitted_periods.iter().any(|p| p.is_some());
            let periods_vec = if any_periodic {
                emitted_periods
            } else {
                Vec::new()
            };
            // Tensor smooths (`te`/`ti`/`t2`) must match mgcv's DEFAULT
            // `select = FALSE`: the joint null space of the per-margin
            // penalties — the bilinear, low-order interaction directions that
            // no marginal roughness operator can see — is left UNPENALIZED.
            // mgcv only adds a null-space shrinkage penalty there under the
            // opt-in `select = TRUE` (which gam exposes as `double_penalty`).
            //
            // The general smooth default (`smooth_double_penalty`, true) is
            // calibrated for 1-D `s()` terms; carrying it into tensors silently
            // shrinks the genuinely-present bilinear interaction signal, so
            // REML places positive weight on the extra ridge and systematically
            // OVER-SMOOTHS the recovered surface relative to mgcv's plain
            // `te`/`ti` (gam#700/#701/#702/#703). Default tensors to no extra
            // null-space penalty; an explicit user `double_penalty=`/`select=`
            // still wins.
            let tensor_double_penalty = option_bool(options, "double_penalty").unwrap_or(false);
            Ok(SmoothBasisSpec::TensorBSpline {
                feature_cols: cols.to_vec(),
                spec: TensorBSplineSpec {
                    marginalspecs: margins,
                    periods: periods_vec,
                    double_penalty: tensor_double_penalty,
                    identifiability: parse_tensor_identifiability(options, kind)?,
                },
            })
        }
        "pca" => {
            validate_known_options(
                "pca",
                options,
                &[
                    "type",
                    "bs",
                    "by",
                    "k",
                    "basis_dim",
                    "basis-dim",
                    "basisdim",
                    "lazy_path",
                    "path",
                    "pca_basis_path",
                    "chunk_size",
                    "smooth_penalty",
                    "centered",
                    "double_penalty",
                    "id",
                    "__by_col",
                ],
            )?;
            let path = options
                .get("lazy_path")
                .or_else(|| options.get("pca_basis_path"))
                .or_else(|| options.get("path"))
                .map(|raw| PathBuf::from(strip_quotes(raw)));
            let Some(path) = path else {
                return Err(TermBuilderError::incompatible_config(
                    "pca smooth requires lazy_path=... on the formula path",
                )
                .to_string());
            };
            let k = option_usize_any(options, &["k", "basis_dim", "basis-dim", "basisdim"])
                .unwrap_or(0);
            let chunk_size = option_usize(options, "chunk_size").unwrap_or(4096);
            Ok(SmoothBasisSpec::Pca {
                feature_cols: cols.to_vec(),
                basis_matrix: Array2::<f64>::zeros((cols.len(), k)),
                centered: option_bool(options, "centered").unwrap_or(true),
                smooth_penalty: option_f64(options, "smooth_penalty").unwrap_or(1.0),
                center_mean: None,
                pca_basis_path: Some(path),
                chunk_size,
            })
        }
        other => Err(TermBuilderError::unsupported_feature(format!(
            "unsupported smooth type '{other}'"
        ))
        .to_string()),
    }
}

/// Initialise per-axis anisotropic log-scales on eligible spatial smooth specs.
pub fn enable_scale_dimensions(spec: &mut TermCollectionSpec) {
    for smooth in spec.smooth_terms.iter_mut() {
        match &mut smooth.basis {
            SmoothBasisSpec::Matern {
                feature_cols,
                spec: matern,
                ..
            } => {
                if matern.aniso_log_scales.is_none() {
                    let d = feature_cols.len();
                    matern.aniso_log_scales = Some(vec![0.0; d]);
                }
            }
            SmoothBasisSpec::Duchon {
                feature_cols,
                spec: duchon,
                ..
            } => {
                if duchon.aniso_log_scales.is_none() {
                    let d = feature_cols.len();
                    duchon.aniso_log_scales = Some(vec![0.0; d]);
                }
            }
            _ => {}
        }
    }
}

// ---------------------------------------------------------------------------
// Data-aware helpers
// ---------------------------------------------------------------------------

pub fn spatial_center_strategy_for_dimension(num_centers: usize, d: usize) -> CenterStrategy {
    if d == 1 {
        // In one dimension, an explicit `k` is a resolution request rather than
        // a request for quantile-midpoint centers. Use the same maximin geometry
        // as the automatic 1D path so Duchon REML sees a well-resolved native
        // kernel block instead of compensating for endpoint under-resolution by
        // over-smoothing easy low-noise signals (#504).
        CenterStrategy::FarthestPoint { num_centers }
    } else {
        default_spatial_center_strategy(num_centers, d)
    }
}

pub fn col_minmax(col: ArrayView1<'_, f64>) -> Result<(f64, f64), String> {
    let min = col.iter().fold(f64::INFINITY, |a, &b| a.min(b));
    let max = col.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));
    if !min.is_finite() || !max.is_finite() {
        return Err(TermBuilderError::degenerate_data(
            "non-finite data encountered while inferring knot range",
        )
        .to_string());
    }
    if (max - min).abs() < 1e-12 {
        Ok((min, min + 1e-6))
    } else {
        Ok((min, max))
    }
}

pub fn unique_count_column(col: ArrayView1<'_, f64>) -> usize {
    use std::collections::HashSet;
    let mut set = HashSet::<u64>::with_capacity(col.len());
    for &v in col {
        let norm = if v == 0.0 { 0.0 } else { v };
        set.insert(norm.to_bits());
    }
    set.len().max(1)
}

/// Smallest number of distinct covariate values seen within any single group
/// of `group_col`. For a factor smooth this is the resolution that bounds the
/// marginal basis: a group with `m` distinct covariate values can only inform
/// `m` basis coefficients, so a marginal richer than that interpolates the
/// group instead of estimating a penalized trend. Bits are compared exactly so
/// integer-valued covariates (days, dose levels) collapse to their true count.
fn min_per_group_unique_count(
    feature_col: ArrayView1<'_, f64>,
    group_col: ArrayView1<'_, f64>,
) -> usize {
    use std::collections::{HashMap, HashSet};
    let mut per_group: HashMap<u64, HashSet<u64>> = HashMap::new();
    for (xi, gi) in feature_col.iter().zip(group_col.iter()) {
        let xnorm = if *xi == 0.0 { 0.0 } else { *xi };
        let gnorm = if *gi == 0.0 { 0.0 } else { *gi };
        per_group
            .entry(gnorm.to_bits())
            .or_default()
            .insert(xnorm.to_bits());
    }
    per_group
        .values()
        .map(|s| s.len())
        .min()
        .unwrap_or(1)
        .max(1)
}

/// Per-column knot count from the unique-value count, with the same n^(1/3)
/// ceiling growth as `heuristic_knots` so per-column smooths can support more
/// detail at biobank scale. The 4-knot floor stays put because we still need
/// enough basis functions to fit a non-trivial smooth at all.
pub fn heuristic_knots_for_column(col: ArrayView1<'_, f64>) -> usize {
    let unique = unique_count_column(col);
    let ceiling = ((unique as f64).cbrt() as usize).max(20);
    (unique / 4).clamp(4, ceiling)
}

pub fn heuristic_centers(n: usize, d: usize) -> usize {
    default_num_centers(n, d)
}

// ---------------------------------------------------------------------------
// Smooth option parsers
// ---------------------------------------------------------------------------

fn parse_endpoint_side(
    value: &str,
    context: &str,
) -> Result<BSplineEndpointBoundaryCondition, String> {
    match value.trim().to_ascii_lowercase().as_str() {
        "" | "none" | "open" | "unconstrained" => Ok(BSplineEndpointBoundaryCondition::Free),
        "clamped" | "clamp" | "zero_derivative" | "zero-derivative" => {
            Ok(BSplineEndpointBoundaryCondition::Clamped)
        }
        "anchored" | "anchor" | "zero" | "zero_value" | "zero-value" => {
            Ok(BSplineEndpointBoundaryCondition::Anchored { value: 0.0 })
        }
        other => Err(format!(
            "unsupported {context} boundary condition '{other}'; expected none, clamped, or anchored"
        )),
    }
}

fn boundary_anchor_value(
    options: &BTreeMap<String, String>,
    side: &str,
    fallback: Option<f64>,
) -> Option<f64> {
    [
        format!("anchor_{side}"),
        format!("{side}_anchor"),
        format!("anchor-value-{side}"),
    ]
    .iter()
    .find_map(|key| option_f64(options, key))
    .or(fallback)
}

fn apply_anchor_value(
    cond: BSplineEndpointBoundaryCondition,
    value: Option<f64>,
) -> BSplineEndpointBoundaryCondition {
    match cond {
        BSplineEndpointBoundaryCondition::Anchored { .. } => {
            BSplineEndpointBoundaryCondition::Anchored {
                value: value.unwrap_or(0.0),
            }
        }
        other => other,
    }
}

fn parse_bspline_boundary_conditions(
    options: &BTreeMap<String, String>,
) -> Result<BSplineBoundaryConditions, String> {
    let fallback_anchor = option_f64(options, "anchor")
        .or_else(|| option_f64(options, "anchor_value"))
        .or_else(|| option_f64(options, "value"));
    let global_boundary_conditions = options
        .get("boundary_conditions")
        .or_else(|| options.get("bc"));
    let mut boundary_conditions = BSplineBoundaryConditions::default();

    if let Some(raw_boundary_conditions) = global_boundary_conditions {
        let cond = parse_endpoint_side(raw_boundary_conditions, "boundary_conditions")?;
        let side = options
            .get("side")
            .map(|s| s.trim().to_ascii_lowercase())
            .unwrap_or_else(|| "both".to_string());
        match side.as_str() {
            "both" | "all" | "endpoints" => {
                boundary_conditions.left = cond;
                boundary_conditions.right = cond;
            }
            "left" | "start" | "lower" => boundary_conditions.left = cond,
            "right" | "end" | "upper" => boundary_conditions.right = cond,
            other => {
                return Err(format!(
                    "unsupported B-spline boundary side '{other}'; expected left, right, or both"
                ));
            }
        }
    }

    if let Some(raw) = options
        .get("bc_left")
        .or_else(|| options.get("left_bc"))
        .or_else(|| options.get("bc_start"))
        .or_else(|| options.get("start_bc"))
    {
        boundary_conditions.left = parse_endpoint_side(raw, "left endpoint")?;
    }
    if let Some(raw) = options
        .get("bc_right")
        .or_else(|| options.get("right_bc"))
        .or_else(|| options.get("bc_end"))
        .or_else(|| options.get("end_bc"))
    {
        boundary_conditions.right = parse_endpoint_side(raw, "right endpoint")?;
    }

    boundary_conditions.left = apply_anchor_value(
        boundary_conditions.left,
        boundary_anchor_value(options, "left", fallback_anchor),
    );
    boundary_conditions.right = apply_anchor_value(
        boundary_conditions.right,
        boundary_anchor_value(options, "right", fallback_anchor),
    );

    // Non-zero anchors require an affine offset term that the current basis
    // builder does not synthesize (see `build_bspline_basis_1d` in
    // src/terms/basis.rs). Surface the rejection at parse time with the side
    // and value in the diagnostic, instead of letting the value-only error
    // emerge deep inside the basis builder where the user has no context
    // about which anchor key (`anchor`, `left_anchor`, `right_anchor`, …)
    // routed into which endpoint.
    reject_nonzero_anchor("left", boundary_conditions.left)?;
    reject_nonzero_anchor("right", boundary_conditions.right)?;

    Ok(boundary_conditions)
}

fn reject_nonzero_anchor(side: &str, cond: BSplineEndpointBoundaryCondition) -> Result<(), String> {
    if let BSplineEndpointBoundaryCondition::Anchored { value } = cond {
        if value.abs() > 1e-12 {
            return Err(format!(
                "non-zero {side} anchor {value} requires an affine offset term that is not yet supported; only anchored value 0 is accepted at parse time"
            ));
        }
    }
    Ok(())
}

pub fn parse_ps_internal_knots(
    options: &BTreeMap<String, String>,
    degree: usize,
    default_internal_knots: usize,
) -> Result<(usize, bool), String> {
    const MIN_EXPRESSIVE_INTERNAL_KNOTS: usize = 2;
    // Strict variants: reject `k=-1`, `k=1.5`, `knots=-2` etc. with a
    // focused error instead of silently dropping the value and using the
    // default. Lenient `option_usize` / `option_usize_any` silently swallow
    // unparseable values, which leaves the user thinking they configured
    // something when they did not.
    // A list-valued `knots=[...]` carries explicit internal positions, not a
    // count; it is consumed by `parse_explicit_internal_knots`. Treat it as
    // "count not specified" here so the strict integer parse does not reject
    // the bracketed value (the Provided path ignores the returned count).
    let knots_internal = if knots_option_is_list(options) {
        None
    } else {
        option_usize_strict(options, "knots")?
    };
    let basis_dim = option_usize_any_strict(options, &["k", "basis_dim", "basis-dim", "basisdim"])?;
    if knots_internal.is_some() && basis_dim.is_some() {
        return Err(TermBuilderError::incompatible_config(
            "ps/bspline smooth: specify either knots=<internal_knots> or k=<basis_dim> (not both)",
        )
        .to_string());
    }
    if let Some(k) = basis_dim {
        let min_k = degree + 1;
        if k < min_k {
            return Err(TermBuilderError::invalid_option(format!(
                "ps/bspline smooth: k={} too small for degree {}; expected k >= {}",
                k, degree, min_k
            ))
            .to_string());
        }
        Ok(((k - min_k).max(MIN_EXPRESSIVE_INTERNAL_KNOTS), false))
    } else {
        Ok((
            knots_internal.unwrap_or(default_internal_knots),
            knots_internal.is_none(),
        ))
    }
}

/// True when the `knots` option value is a *list* literal (`[...]`, `c(...)`,
/// or `(...)`) rather than a scalar count. mgcv's `knots=` accepts both: a
/// single integer is an internal-knot count, while a vector is explicit
/// internal knot positions. We disambiguate purely on the wrapper syntax so a
/// bare `knots=5` keeps its historical count meaning.
fn knots_option_is_list(options: &BTreeMap<String, String>) -> bool {
    options
        .get("knots")
        .map(|raw| {
            let t = raw.trim();
            t.starts_with('[') || t.starts_with("c(") || t.starts_with("C(") || t.starts_with('(')
        })
        .unwrap_or(false)
}

/// Parse `knots=[k0, k1, ...]` (or `c(...)` / `(...)`) into explicit internal
/// knot positions. Returns `Ok(None)` when `knots` is absent or a scalar count
/// (handled by [`parse_ps_internal_knots`]); `Ok(Some(positions))` when it is a
/// non-empty numeric list; and an error for an empty or unparseable list.
fn parse_explicit_internal_knots(
    options: &BTreeMap<String, String>,
) -> Result<Option<Vec<f64>>, String> {
    if !knots_option_is_list(options) {
        return Ok(None);
    }
    let raw = options
        .get("knots")
        .expect("knots_option_is_list implies the key is present");
    // Unwrap the `c(...)` / `(...)` R-vector forms down to a bare `[...]`-style
    // body so `split_list_option` (which only peels `[...]`) sees a clean list.
    let trimmed = raw.trim();
    let body = trimmed
        .strip_prefix("c(")
        .or_else(|| trimmed.strip_prefix("C("))
        .or_else(|| trimmed.strip_prefix('('))
        .and_then(|v| v.strip_suffix(')'))
        .map(|inner| format!("[{inner}]"))
        .unwrap_or_else(|| trimmed.to_string());
    let tokens = split_list_option(&body);
    if tokens.is_empty() {
        return Err(TermBuilderError::invalid_option(format!(
            "knots={raw} is an empty list; supply at least one internal knot position \
             (e.g. knots=[0.2, 0.5, 0.8]) or a scalar count (e.g. knots=8)"
        ))
        .to_string());
    }
    let mut positions = Vec::with_capacity(tokens.len());
    for tok in &tokens {
        let value = parse_numeric_expr(tok).map_err(|err| {
            TermBuilderError::invalid_option(format!(
                "knots list entry '{tok}' is not a numeric position: {err}"
            ))
            .to_string()
        })?;
        positions.push(value);
    }
    Ok(Some(positions))
}

/// Resolve the `knot_placement=` option for an automatically generated knot
/// vector. Accepts `"uniform"` (the default, equal spacing on the data range)
/// and `"quantile"` (interior knots at empirical data quantiles, better for
/// skewed covariates). Unknown values are rejected so typos do not silently
/// fall back to uniform.
fn parse_knot_placement(
    options: &BTreeMap<String, String>,
) -> Result<crate::basis::BSplineKnotPlacement, String> {
    use crate::basis::BSplineKnotPlacement;
    match options
        .get("knot_placement")
        .or_else(|| options.get("knot-placement"))
        .or_else(|| options.get("knotplacement"))
    {
        None => Ok(BSplineKnotPlacement::Uniform),
        Some(raw) => match raw
            .trim()
            .trim_matches('"')
            .trim_matches('\'')
            .to_ascii_lowercase()
            .as_str()
        {
            "uniform" | "even" | "equal" => Ok(BSplineKnotPlacement::Uniform),
            "quantile" | "quantiles" | "data" | "empirical" => Ok(BSplineKnotPlacement::Quantile),
            other => Err(TermBuilderError::invalid_option(format!(
                "knot_placement={other} is not recognised; expected \"uniform\" or \"quantile\""
            ))
            .to_string()),
        },
    }
}

/// Build the non-periodic 1D B-spline knot spec for the `ps`/`bspline` and
/// factor-smooth marginal paths, honoring (in priority order):
///   1. `knots=[...]` explicit internal positions  → [`BSplineKnotSpec::Provided`]
///   2. `knot_placement="quantile"`                 → [`BSplineKnotSpec::Automatic`]
///   3. uniform generation                          → [`BSplineKnotSpec::Generate`]
///
/// `data` is the covariate column (used to clamp explicit positions to the
/// observed range and to drive quantile placement); `n_knots` is the resolved
/// internal-knot count from [`parse_ps_internal_knots`] used for the automatic
/// strategies.
fn resolve_nonperiodic_bspline_knotspec(
    options: &BTreeMap<String, String>,
    data: ArrayView1<'_, f64>,
    data_range: (f64, f64),
    degree: usize,
    n_knots: usize,
) -> Result<BSplineKnotSpec, String> {
    use crate::basis::{BSplineKnotPlacement, clamped_knot_vector_from_internal_positions};
    if let Some(positions) = parse_explicit_internal_knots(options)? {
        if option_usize_any_strict(options, &["k", "basis_dim", "basis-dim", "basisdim"])?.is_some()
        {
            return Err(TermBuilderError::incompatible_config(
                "ps/bspline smooth: specify either explicit knots=[...] positions or \
                 k=<basis_dim> (not both); the basis size is fixed by the knot vector",
            )
            .to_string());
        }
        let knots = clamped_knot_vector_from_internal_positions(data_range, &positions, degree)
            .map_err(|e| e.to_string())?;
        return Ok(BSplineKnotSpec::Provided(knots));
    }
    match parse_knot_placement(options)? {
        BSplineKnotPlacement::Uniform => Ok(BSplineKnotSpec::Generate {
            data_range,
            num_internal_knots: n_knots,
        }),
        BSplineKnotPlacement::Quantile => {
            // Validate the column up-front so an unfittable request surfaces a
            // user-correctable error at parse time rather than deep in basis
            // construction. The same data drives the eventual quantile knots.
            crate::basis::auto_knot_vector_1d_quantile(data, n_knots, degree)
                .map_err(|e| e.to_string())?;
            Ok(BSplineKnotSpec::Automatic {
                num_internal_knots: Some(n_knots),
                placement: BSplineKnotPlacement::Quantile,
            })
        }
    }
}

/// Reject unknown option keys with a focused error that names the term and
/// the offending key, plus suggests near-matches from the known-key list.
/// Without this, typos like `lengt_scale=0.1` or `nyu=5/2` are silently
/// dropped, the term uses the default, and the user has no idea why their
/// option had no effect.
pub fn validate_known_options(
    term_name: &str,
    options: &BTreeMap<String, String>,
    known: &[&str],
) -> Result<(), String> {
    let known_set: std::collections::BTreeSet<&&str> = known.iter().collect();
    for key in options.keys() {
        if !known_set.contains(&key.as_str()) {
            // Suggest near-matches (substring or shared prefix ≥ 3).
            let key_l = key.to_ascii_lowercase();
            let mut suggestions: Vec<&str> = known
                .iter()
                .filter(|k| {
                    let kl = k.to_ascii_lowercase();
                    kl.contains(&key_l) || key_l.contains(&kl) || {
                        let n = kl
                            .chars()
                            .zip(key_l.chars())
                            .take_while(|(a, b)| a == b)
                            .count();
                        n >= 3
                    }
                })
                .copied()
                .collect();
            suggestions.sort_unstable();
            suggestions.dedup();
            let hint = if suggestions.is_empty() {
                String::new()
            } else {
                format!(" — did you mean one of [{}]?", suggestions.join(", "))
            };
            return Err(TermBuilderError::invalid_option(format!(
                "{term_name}() does not accept option `{key}`{hint} Known options: [{}]",
                {
                    let mut sorted = known.to_vec();
                    sorted.sort_unstable();
                    sorted.join(", ")
                }
            ))
            .to_string());
        }
    }
    Ok(())
}

/// Private (engine-injected) option that caps the *default* spatial center
/// count for a secondary (distributional) predictor's smooth — see
/// `solver::workflow::apply_secondary_predictor_basis_parsimony` and #501.
///
/// It is deliberately NOT one of the user-facing count aliases recognised by
/// [`has_explicit_countwith_basis_alias`], so it never flips the spatial basis
/// onto the explicit (hard) center-placement strategy: the cap lowers the
/// *default* count while the `Auto` strategy is retained, so the count is still
/// softly reduced when the data can't support it.
pub(crate) const SECONDARY_CENTER_CAP_OPTION: &str = "__secondary_center_cap";

/// Apply the secondary-predictor center cap to a *default* spatial center
/// count. A no-op when the cap option is absent (the common case) or when the
/// user supplied an explicit count (then `default_count` is ignored downstream
/// by [`parse_countwith_basis_alias`] anyway).
pub(crate) fn cap_default_spatial_centers(
    options: &BTreeMap<String, String>,
    default_count: usize,
) -> usize {
    match option_usize(options, SECONDARY_CENTER_CAP_OPTION) {
        Some(cap) => default_count.min(cap),
        None => default_count,
    }
}

pub fn parse_countwith_basis_alias(
    options: &BTreeMap<String, String>,
    primarykey: &str,
    default_count: usize,
) -> Result<usize, String> {
    // Strict: reject unparseable values (e.g. `centers=many`, `centers=-1`,
    // `centers=1.5`) instead of silently dropping them and falling through
    // to the default. Without this the user gets the auto-inferred count
    // silently and never realizes their explicit option was ignored.
    let primary = option_usize_strict(options, primarykey)?;
    let basis_dim = option_usize_any_strict(
        options,
        &["k", "basis_dim", "basis-dim", "basisdim", "knots"],
    )?;
    if primary.is_some() && basis_dim.is_some() {
        return Err(TermBuilderError::incompatible_config(format!(
            "specify either {}=<count> or k=<basis_dim> (not both)",
            primarykey
        ))
        .to_string());
    }
    Ok(primary.or(basis_dim).unwrap_or(default_count))
}

pub(crate) fn has_explicit_countwith_basis_alias(
    options: &BTreeMap<String, String>,
    primarykey: &str,
) -> bool {
    options.contains_key(primarykey)
        || ["k", "basis_dim", "basis-dim", "basisdim", "knots"]
            .iter()
            .any(|alias| options.contains_key(*alias))
}

pub fn parse_cyclic_boundary(
    options: &BTreeMap<String, String>,
    minv: f64,
    maxv: f64,
) -> Result<OneDimensionalBoundary, String> {
    let cyclic = option_bool(options, "cyclic")
        .or_else(|| option_bool(options, "periodic"))
        .unwrap_or(false);
    if !cyclic {
        return Ok(OneDimensionalBoundary::Open);
    }
    let start = match option_numeric_expr(options, "period_start")? {
        Some(v) => v,
        None => option_numeric_expr(options, "start")?.unwrap_or(minv),
    };
    let end = match option_numeric_expr(options, "period_end")? {
        Some(v) => v,
        None => option_numeric_expr(options, "end")?.unwrap_or(maxv),
    };
    if end <= start {
        return Err(format!(
            "cyclic smooth requires period_end/end ({end}) > period_start/start ({start})"
        ));
    }
    Ok(OneDimensionalBoundary::Cyclic { start, end })
}

/// Parse the periodic-uniform domain for a one-dimensional cyclic smooth.
///
/// Returns the `(domain_start, period)` pair derived from
/// `period_start` / `start`, `period_end` / `end`, falling back to the
/// data range `[minv, maxv)` when neither bound is provided. The period
/// must be strictly positive.
pub fn parse_periodic_domain_1d(
    options: &BTreeMap<String, String>,
    minv: f64,
    maxv: f64,
) -> Result<(f64, f64), String> {
    let start = match option_numeric_expr(options, "period_start")? {
        Some(v) => v,
        None => option_numeric_expr(options, "start")?.unwrap_or(minv),
    };
    let end = match option_numeric_expr(options, "period_end")? {
        Some(v) => v,
        None => option_numeric_expr(options, "end")?.unwrap_or(maxv),
    };
    if !(start.is_finite() && end.is_finite()) {
        return Err(format!(
            "periodic smooth domain requires finite endpoints, got ({start}, {end})"
        ));
    }
    if end <= start {
        return Err(format!(
            "periodic smooth requires period_end/end ({end}) > period_start/start ({start})"
        ));
    }
    Ok((start, end - start))
}

fn parse_matern_nu(raw: &str) -> Result<MaternNu, String> {
    let trimmed = raw.trim();
    let lowered = trimmed.to_ascii_lowercase();
    match lowered.as_str() {
        "1/2" | "0.5" | "half" => return Ok(MaternNu::Half),
        "3/2" | "1.5" => return Ok(MaternNu::ThreeHalves),
        "5/2" | "2.5" => return Ok(MaternNu::FiveHalves),
        "7/2" | "3.5" => return Ok(MaternNu::SevenHalves),
        "9/2" | "4.5" => return Ok(MaternNu::NineHalves),
        _ => {}
    }

    let value = if let Some((num, den)) = trimmed.split_once('/') {
        let num = num
            .trim()
            .parse::<f64>()
            .map_err(|err| format!("{}: {err}", unsupported_matern_nu_message(raw)))?;
        let den = den
            .trim()
            .parse::<f64>()
            .map_err(|err| format!("{}: {err}", unsupported_matern_nu_message(raw)))?;
        if den == 0.0 || !num.is_finite() || !den.is_finite() {
            return Err(unsupported_matern_nu_message(raw));
        }
        num / den
    } else {
        trimmed
            .parse::<f64>()
            .map_err(|err| format!("{}: {err}", unsupported_matern_nu_message(raw)))?
    };

    const TOL: f64 = 1e-12;
    if (value - 0.5).abs() <= TOL {
        Ok(MaternNu::Half)
    } else if (value - 1.5).abs() <= TOL {
        Ok(MaternNu::ThreeHalves)
    } else if (value - 2.5).abs() <= TOL {
        Ok(MaternNu::FiveHalves)
    } else if (value - 3.5).abs() <= TOL {
        Ok(MaternNu::SevenHalves)
    } else if (value - 4.5).abs() <= TOL {
        Ok(MaternNu::NineHalves)
    } else {
        Err(unsupported_matern_nu_message(raw))
    }
}

fn unsupported_matern_nu_message(raw: &str) -> String {
    TermBuilderError::unsupported_feature(format!(
        "unsupported Matern nu '{raw}'; supported half-integer values are 1/2, 3/2, 5/2, 7/2, and 9/2"
    ))
    .to_string()
}

#[derive(Clone, Debug, serde::Serialize, serde::Deserialize)]
pub enum DuchonPowerPolicy {
    Explicit(f64),
    /// No explicit `power=` given: defer to the cubic structural default, which
    /// the builder resolves dimension-aware as `s = (d − 1)/2` (so `φ(r) = r³`
    /// in every dimension). There is no triple-operator minimum any more.
    CubicStructuralDefault,
}

pub fn parse_duchon_power_policy(
    options: &BTreeMap<String, String>,
) -> Result<DuchonPowerPolicy, String> {
    if let Some(raw_nu) = options.get("nu") {
        return Err(TermBuilderError::incompatible_config(format!(
            "Duchon smooths use power=<number>, not nu='{}'. Use power=1.5, power=2, etc.",
            raw_nu
        ))
        .to_string());
    }
    match options.get("power") {
        Some(raw) => {
            let value = raw.parse::<f64>().map_err(|err| {
                TermBuilderError::invalid_option(format!(
                    "invalid Duchon power '{}'; expected a non-negative number such as power=1.5 or power=2: {}",
                    raw, err
                ))
                .to_string()
            })?;
            if !value.is_finite() || value < 0.0 {
                return Err(TermBuilderError::invalid_option(format!(
                    "invalid Duchon power '{}'; expected a finite non-negative number such as power=1.5 or power=2",
                    raw
                ))
                .to_string());
            }
            Ok(DuchonPowerPolicy::Explicit(value))
        }
        None => Ok(DuchonPowerPolicy::CubicStructuralDefault),
    }
}

pub fn parse_duchon_power(options: &BTreeMap<String, String>) -> Result<f64, String> {
    match parse_duchon_power_policy(options)? {
        DuchonPowerPolicy::Explicit(power) => Ok(power),
        // Context-free placeholder: the bare option parser has no column count,
        // so it cannot compute the dimension-aware cubic power `s = (d − 1)/2`.
        // The dimension-aware resolution happens later in `build_smooth_basis`;
        // this 1.5 is only a stand-in for callers that need a concrete number
        // without data context (e.g. round-trip parser tests).
        DuchonPowerPolicy::CubicStructuralDefault => Ok(1.5),
    }
}

pub fn parse_duchon_order(
    options: &BTreeMap<String, String>,
) -> Result<DuchonNullspaceOrder, String> {
    match options.get("order") {
        // Structural cubic Duchon is affine-by-default: an unspecified order is
        // the `Linear` (constant + linear) null space, matching the magic
        // default. An explicit `order=0` still selects the constant-only space.
        None => Ok(DuchonNullspaceOrder::Linear),
        Some(raw) => match raw.parse::<usize>() {
            Ok(0) => Ok(DuchonNullspaceOrder::Zero),
            Ok(1) => Ok(DuchonNullspaceOrder::Linear),
            Ok(other) => Ok(DuchonNullspaceOrder::Degree(other)),
            Err(_) => Err(TermBuilderError::invalid_option(format!(
                "invalid Duchon order '{}'; expected a non-negative integer such as order=0, order=1, or order=2",
                raw
            ))
            .to_string()),
        },
    }
}

fn parse_matern_identifiability(
    options: &BTreeMap<String, String>,
) -> Result<MaternIdentifiability, TermBuilderError> {
    let Some(raw) = options.get("identifiability").map(String::as_str) else {
        return Ok(MaternIdentifiability::default());
    };
    match raw.trim().to_ascii_lowercase().as_str() {
        "none" => Ok(MaternIdentifiability::None),
        "sum_tozero" | "sum-to-zero" | "center_sum_tozero" | "center-sum-to-zero" | "centered" => {
            Ok(MaternIdentifiability::CenterSumToZero)
        }
        "linear" | "center_linear_orthogonal" | "center-linear-orthogonal" => {
            Ok(MaternIdentifiability::CenterLinearOrthogonal)
        }
        other => Err(TermBuilderError::unsupported_feature(format!(
            "invalid Matérn identifiability '{other}'; expected one of: none, sum_tozero, linear"
        ))),
    }
}

fn parse_spatial_identifiability(
    options: &BTreeMap<String, String>,
) -> Result<SpatialIdentifiability, TermBuilderError> {
    let Some(raw) = options.get("identifiability").map(String::as_str) else {
        return Ok(SpatialIdentifiability::default());
    };
    match raw.trim().to_ascii_lowercase().as_str() {
        "none" => Ok(SpatialIdentifiability::None),
        "orthogonal"
        | "orthogonal_to_parametric"
        | "orthogonal-to-parametric"
        | "parametric_orthogonal" => Ok(SpatialIdentifiability::OrthogonalToParametric),
        "frozen" => Err(TermBuilderError::unsupported_feature(
            "spatial identifiability 'frozen' is internal-only; use none or orthogonal_to_parametric",
        )),
        other => Err(TermBuilderError::unsupported_feature(format!(
            "invalid spatial identifiability '{other}'; expected one of: none, orthogonal_to_parametric"
        ))),
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::inference::formula_dsl::parse_formula;
    use crate::inference::model::{DataSchema, SchemaColumn};
    use ndarray::Array2;
    use std::collections::BTreeMap;

    fn continuous_dataset(headers: &[&str], rows: Vec<Vec<f64>>) -> Dataset {
        let nrows = rows.len();
        let ncols = headers.len();
        let values = Array2::from_shape_vec(
            (nrows, ncols),
            rows.into_iter().flat_map(|row| row.into_iter()).collect(),
        )
        .expect("rectangular test data");
        Dataset {
            headers: headers.iter().map(|name| name.to_string()).collect(),
            values,
            schema: DataSchema {
                columns: headers
                    .iter()
                    .map(|name| SchemaColumn {
                        name: name.to_string(),
                        kind: ColumnKindTag::Continuous,
                        levels: vec![],
                    })
                    .collect(),
            },
            column_kinds: vec![ColumnKindTag::Continuous; ncols],
        }
    }

    fn inferred_tensor_basis_product(ds: &Dataset) -> usize {
        let parsed = parse_formula("y ~ te(theta, h)").expect("parse tensor formula");
        let col_map = ds.column_map();
        let mut notes = Vec::new();
        let terms = build_termspec(
            &parsed.terms,
            ds,
            &col_map,
            &mut notes,
            &ResourcePolicy::default_library(),
        )
        .expect("build tensor termspec");
        let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
            panic!("expected tensor smooth");
        };
        spec.marginalspecs
            .iter()
            .map(|marginal| match marginal.knotspec {
                BSplineKnotSpec::Generate {
                    num_internal_knots, ..
                } => num_internal_knots + marginal.degree + 1,
                BSplineKnotSpec::PeriodicUniform { num_basis, .. } => num_basis,
                BSplineKnotSpec::Automatic {
                    num_internal_knots: Some(num_internal_knots),
                    ..
                } => num_internal_knots + marginal.degree + 1,
                BSplineKnotSpec::Automatic {
                    num_internal_knots: None,
                    ..
                } => panic!("test helper cannot infer automatic knot count"),
                BSplineKnotSpec::Provided(ref knots) => {
                    knots.len().saturating_sub(marginal.degree + 1)
                }
            })
            .product()
    }

    #[test]
    fn parse_cylinder_periodic_options_match_requested_forms() {
        let mut opts = BTreeMap::new();
        opts.insert("periodic".to_string(), "[0]".to_string());
        opts.insert("period".to_string(), "[2*pi, None]".to_string());
        let axes = parse_periodic_axes(&opts, 2).expect("axes");
        let periods = parse_periods(&opts, &axes).expect("periods");
        assert_eq!(axes, vec![true, false]);
        assert!((periods[0].unwrap() - 2.0 * std::f64::consts::PI).abs() < 1e-12);
        assert_eq!(periods[1], None);

        let mut boundary_opts = BTreeMap::new();
        boundary_opts.insert(
            "boundary".to_string(),
            "['periodic', 'natural']".to_string(),
        );
        boundary_opts.insert("period".to_string(), "[2*pi, None]".to_string());
        let boundary_axes = parse_periodic_axes(&boundary_opts, 2).expect("boundary axes");
        let boundary_periods =
            parse_periods(&boundary_opts, &boundary_axes).expect("boundary periods");
        assert_eq!(boundary_axes, vec![true, false]);
        assert!((boundary_periods[0].unwrap() - 2.0 * std::f64::consts::PI).abs() < 1e-12);
        assert_eq!(boundary_periods[1], None);

        let mut unicode_opts = BTreeMap::new();
        unicode_opts.insert("periodic".to_string(), "[0,1]".to_string());
        unicode_opts.insert("period".to_string(), "[2π, τ]".to_string());
        let unicode_axes = parse_periodic_axes(&unicode_opts, 2).expect("unicode axes");
        let unicode_periods = parse_periods(&unicode_opts, &unicode_axes).expect("unicode periods");
        assert_eq!(unicode_axes, vec![true, true]);
        assert!((unicode_periods[0].unwrap() - 2.0 * std::f64::consts::PI).abs() < 1e-12);
        assert!((unicode_periods[1].unwrap() - std::f64::consts::TAU).abs() < 1e-12);
    }

    #[test]
    fn parse_single_axis_periodic_zero_as_axis_not_false() {
        let mut opts = BTreeMap::new();
        opts.insert("periodic".to_string(), "[0]".to_string());
        opts.insert("period".to_string(), "2*pi".to_string());
        opts.insert("origin".to_string(), "0".to_string());
        let axes = parse_periodic_axes(&opts, 1).expect("axes");
        let periods = parse_periods(&opts, &axes).expect("periods");
        let origins = parse_period_origins(&opts, &axes).expect("origins");
        assert_eq!(axes, vec![true]);
        assert!((periods[0].unwrap() - 2.0 * std::f64::consts::PI).abs() < 1e-12);
        assert_eq!(origins[0], Some(0.0));
    }

    #[test]
    fn one_dimensional_bspline_accepts_boundary_periodic() {
        let ds = continuous_dataset(
            &["y", "theta"],
            (0..16)
                .map(|i| {
                    let theta = std::f64::consts::TAU * i as f64 / 16.0;
                    vec![theta.sin(), theta]
                })
                .collect(),
        );
        let parsed = parse_formula("y ~ s(theta, boundary=periodic, period=2*pi, origin=0, k=8)")
            .expect("parse");
        let col_map = ds.column_map();
        let mut notes = Vec::new();
        let terms = build_termspec(
            &parsed.terms,
            &ds,
            &col_map,
            &mut notes,
            &crate::resource::ResourcePolicy::default_library(),
        )
        .expect("periodic boundary should build");
        let SmoothBasisSpec::BSpline1D { spec, .. } = &terms.smooth_terms[0].basis else {
            panic!("expected 1D B-spline");
        };
        assert!(matches!(
            &spec.knotspec,
            BSplineKnotSpec::PeriodicUniform {
                data_range,
                num_basis: 8
            } if *data_range == (0.0, std::f64::consts::TAU)
        ));
    }

    #[test]
    fn matern_formula_defaults_to_no_double_penalty_but_honors_explicit() {
        // gam#787: the Matérn kernel is strictly positive-definite (no structural
        // polynomial nullspace), so its double-penalty nullspace ridge is spurious
        // and was flipping the learned-penalty count across kappa rebuilds. The
        // formula default must therefore be `double_penalty=false` for Matérn, and
        // an explicit user `double_penalty=true` must still win.
        let ds = continuous_dataset(
            &["y", "x1", "x2"],
            (0..40)
                .map(|i| {
                    let t = i as f64 / 40.0;
                    vec![t.sin(), t, (2.0 * t).cos()]
                })
                .collect(),
        );
        let col_map = ds.column_map();
        let policy = crate::resource::ResourcePolicy::default_library();

        let parsed_default =
            parse_formula("y ~ matern(x1, x2, centers=6)").expect("parse default matern");
        let mut notes = Vec::new();
        let terms_default =
            build_termspec(&parsed_default.terms, &ds, &col_map, &mut notes, &policy)
                .expect("default matern should build");
        let SmoothBasisSpec::Matern { spec, .. } = &terms_default.smooth_terms[0].basis else {
            panic!("expected Matérn basis");
        };
        assert!(
            !spec.double_penalty,
            "Matérn (strictly PD, no structural nullspace) must default to double_penalty=false"
        );

        let parsed_explicit = parse_formula("y ~ matern(x1, x2, centers=6, double_penalty=true)")
            .expect("parse explicit matern");
        let mut notes_explicit = Vec::new();
        let terms_explicit = build_termspec(
            &parsed_explicit.terms,
            &ds,
            &col_map,
            &mut notes_explicit,
            &policy,
        )
        .expect("explicit double_penalty matern should build");
        let SmoothBasisSpec::Matern { spec, .. } = &terms_explicit.smooth_terms[0].basis else {
            panic!("expected Matérn basis");
        };
        assert!(
            spec.double_penalty,
            "an explicit user double_penalty=true must still be honored for Matérn"
        );
    }

    #[test]
    fn formula_shape_constraint_round_trips_and_rejects_bogus() {
        let ds = continuous_dataset(
            &["y", "x"],
            (0..32)
                .map(|i| {
                    let x = i as f64 / 31.0;
                    vec![x * x, x]
                })
                .collect(),
        );
        let col_map = ds.column_map();

        let parsed =
            parse_formula("y ~ s(x, shape=monotone_increasing)").expect("parse monotone smooth");
        let mut notes = Vec::new();
        let terms = build_termspec(
            &parsed.terms,
            &ds,
            &col_map,
            &mut notes,
            &crate::resource::ResourcePolicy::default_library(),
        )
        .expect("monotone smooth should build");
        assert_eq!(
            terms.smooth_terms[0].shape,
            ShapeConstraint::MonotoneIncreasing
        );

        let parsed_bad = parse_formula("y ~ s(x, shape=bogus)").expect("parse bogus shape");
        let mut notes_bad = Vec::new();
        let err = build_termspec(
            &parsed_bad.terms,
            &ds,
            &col_map,
            &mut notes_bad,
            &crate::resource::ResourcePolicy::default_library(),
        )
        .expect_err("bogus shape must error");
        assert!(
            format!("{err:?}").contains("unknown shape constraint"),
            "got: {err:?}"
        );
    }

    #[test]
    fn default_sphere_smooth_uses_spherical_farthest_point_centers() {
        let ds = continuous_dataset(
            &["y", "lat", "lon"],
            (0..24)
                .map(|i| {
                    let t = i as f64 / 24.0;
                    let lat = -60.0 + 120.0 * t;
                    let lon = -180.0 + 360.0 * ((7 * i) % 24) as f64 / 24.0;
                    vec![lat.to_radians().sin(), lat, lon]
                })
                .collect(),
        );
        let parsed = parse_formula("y ~ sphere(lat, lon)").expect("parse");
        let col_map = ds.column_map();
        let mut notes = Vec::new();
        let terms = build_termspec(
            &parsed.terms,
            &ds,
            &col_map,
            &mut notes,
            &crate::resource::ResourcePolicy::default_library(),
        )
        .expect("build sphere termspec");
        let SmoothBasisSpec::Sphere { spec, .. } = &terms.smooth_terms[0].basis else {
            panic!("expected sphere term");
        };
        assert!(matches!(
            spec.center_strategy,
            CenterStrategy::FarthestPoint { .. }
        ));
    }

    #[test]
    fn one_dimensional_duchon_defaults_to_scale_free_length_scale() {
        let ds = continuous_dataset(
            &["y", "x"],
            (0..32)
                .map(|i| {
                    let x = i as f64 / 31.0;
                    vec![(std::f64::consts::TAU * x).sin(), x]
                })
                .collect(),
        );
        let parsed = parse_formula("y ~ duchon(x)").expect("parse");
        let col_map = ds.column_map();
        let mut notes = Vec::new();
        let terms = build_termspec(
            &parsed.terms,
            &ds,
            &col_map,
            &mut notes,
            &crate::resource::ResourcePolicy::default_library(),
        )
        .expect("build default duchon termspec");
        let SmoothBasisSpec::Duchon { spec, .. } = &terms.smooth_terms[0].basis else {
            panic!("expected Duchon term");
        };
        assert_eq!(spec.length_scale, None);
    }

    #[test]
    fn one_dimensional_duchon_length_scale_opts_into_hybrid_mode() {
        let ds = continuous_dataset(
            &["y", "x"],
            (0..32)
                .map(|i| {
                    let x = i as f64 / 31.0;
                    vec![(std::f64::consts::TAU * x).sin(), x]
                })
                .collect(),
        );
        let parsed = parse_formula("y ~ duchon(x, length_scale=0.25)").expect("parse");
        let col_map = ds.column_map();
        let mut notes = Vec::new();
        let terms = build_termspec(
            &parsed.terms,
            &ds,
            &col_map,
            &mut notes,
            &crate::resource::ResourcePolicy::default_library(),
        )
        .expect("build hybrid duchon termspec");
        let SmoothBasisSpec::Duchon { spec, .. } = &terms.smooth_terms[0].basis else {
            panic!("expected Duchon term");
        };
        assert_eq!(spec.length_scale, Some(0.25));
    }

    #[test]
    fn parse_matern_nu_accepts_equivalent_half_integer_forms() {
        let cases = [
            ("1/2", MaternNu::Half),
            (" 1 / 2 ", MaternNu::Half),
            (".5", MaternNu::Half),
            ("0.50", MaternNu::Half),
            ("half", MaternNu::Half),
            ("3 / 2", MaternNu::ThreeHalves),
            ("1.50", MaternNu::ThreeHalves),
            ("5 / 2", MaternNu::FiveHalves),
            ("2.500000000000", MaternNu::FiveHalves),
            ("7 / 2", MaternNu::SevenHalves),
            ("3.50", MaternNu::SevenHalves),
            ("9 / 2", MaternNu::NineHalves),
            ("4.50", MaternNu::NineHalves),
        ];
        for (raw, expected) in cases {
            let parsed = parse_matern_nu(raw).expect(raw);
            assert!(
                matches!(
                    (parsed, expected),
                    (MaternNu::Half, MaternNu::Half)
                        | (MaternNu::ThreeHalves, MaternNu::ThreeHalves)
                        | (MaternNu::FiveHalves, MaternNu::FiveHalves)
                        | (MaternNu::SevenHalves, MaternNu::SevenHalves)
                        | (MaternNu::NineHalves, MaternNu::NineHalves)
                ),
                "parsed {raw:?} as {parsed:?}, expected {expected:?}"
            );
        }
    }

    #[test]
    fn parse_matern_nu_rejects_unsupported_or_invalid_values() {
        for raw in ["1", "2", "11/2", "1/0", "nan", "fast"] {
            let err = parse_matern_nu(raw).expect_err(raw);
            assert!(
                err.contains("supported half-integer values"),
                "unexpected error for {raw:?}: {err}"
            );
        }
    }

    #[test]
    fn parse_ps_k_promotes_underexpressive_cubic_basis() {
        let mut opts = BTreeMap::new();
        opts.insert("k".to_string(), "4".to_string());
        let (internal, inferred) = parse_ps_internal_knots(&opts, 3, 20).expect("k=4");
        assert_eq!(internal, 2);
        assert!(!inferred);

        opts.insert("k".to_string(), "6".to_string());
        let (internal, inferred) = parse_ps_internal_knots(&opts, 3, 20).expect("k=6");
        assert_eq!(internal, 2);
        assert!(!inferred);

        opts.insert("k".to_string(), "10".to_string());
        let (internal, inferred) = parse_ps_internal_knots(&opts, 3, 20).expect("k=10");
        assert_eq!(internal, 6);
        assert!(!inferred);
    }

    #[test]
    fn parse_tensor_periods_and_origins_aliases() {
        let mut opts = BTreeMap::new();
        opts.insert(
            "boundary".to_string(),
            "['periodic', 'periodic']".to_string(),
        );
        opts.insert("periods".to_string(), "[7, 24]".to_string());
        opts.insert("origins".to_string(), "[0, -12]".to_string());
        let axes = parse_periodic_axes(&opts, 2).expect("axes");
        let periods = parse_periods(&opts, &axes).expect("periods");
        let origins = parse_period_origins(&opts, &axes).expect("origins");
        assert_eq!(axes, vec![true, true]);
        assert_eq!(periods, vec![Some(7.0), Some(24.0)]);
        assert_eq!(origins, vec![Some(0.0), Some(-12.0)]);
    }

    #[test]
    fn tensor_smooth_honors_per_margin_k_list() {
        let ds = continuous_dataset(
            &["y", "theta", "h"],
            (0..20)
                .map(|i| {
                    let theta = std::f64::consts::TAU * i as f64 / 20.0;
                    let h = -1.0 + 2.0 * (i % 5) as f64 / 4.0;
                    vec![theta.cos() + h, theta, h]
                })
                .collect(),
        );
        let parsed = parse_formula(
            "y ~ te(theta, h, periodic=[0], period=[2*pi, None], origin=[0, None], k=[9,5])",
        )
        .expect("parse tensor formula");
        let col_map = ds.column_map();
        let mut notes = Vec::new();
        let terms = build_termspec(
            &parsed.terms,
            &ds,
            &col_map,
            &mut notes,
            &crate::resource::ResourcePolicy::default_library(),
        )
        .expect("build tensor terms");
        let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
            panic!("expected tensor B-spline");
        };
        let dims = spec
            .marginalspecs
            .iter()
            .map(|m| match m.knotspec {
                BSplineKnotSpec::PeriodicUniform { num_basis, .. } => num_basis,
                BSplineKnotSpec::Generate {
                    num_internal_knots, ..
                } => num_internal_knots + m.degree + 1,
                _ => panic!("unexpected tensor marginal knotspec"),
            })
            .collect::<Vec<_>>();
        assert_eq!(dims, vec![9, 5]);
    }

    #[test]
    fn explicit_basis_sizes_are_not_small_n_clamped() {
        let ds = continuous_dataset(
            &["y", "x1", "x2", "x3", "x4", "x5"],
            (0..12)
                .map(|i| {
                    let x = i as f64 / 11.0;
                    vec![x.sin(), x, x * x, x + 0.1, 1.0 - x, (2.0 * x).sin()]
                })
                .collect(),
        );
        let parsed = parse_formula("y ~ s(x1, k=10) + s(x2) + s(x3) + s(x4) + s(x5)")
            .expect("parse multi-smooth formula");
        let col_map = ds.column_map();
        let mut notes = Vec::new();
        let terms = build_termspec(
            &parsed.terms,
            &ds,
            &col_map,
            &mut notes,
            &crate::resource::ResourcePolicy::default_library(),
        )
        .expect("build multi-smooth terms");
        let SmoothBasisSpec::BSpline1D { spec, .. } = &terms.smooth_terms[0].basis else {
            panic!("expected first smooth to be B-spline");
        };
        assert!(matches!(
            &spec.knotspec,
            BSplineKnotSpec::Generate {
                num_internal_knots: 6,
                ..
            }
        ));
    }

    #[test]
    fn explicit_duchon_centers_are_not_small_n_bumped() {
        let ds = continuous_dataset(
            &["y", "x1", "x2", "x3", "x4", "x5"],
            (0..12)
                .map(|i| {
                    let x = i as f64 / 11.0;
                    vec![x.sin(), x, x * x, x + 0.1, 1.0 - x, (2.0 * x).sin()]
                })
                .collect(),
        );
        // Pure 1D Duchon at default options resolves the nullspace to Linear
        // (2s < d forces escalation), giving 2 polynomial nullspace columns;
        // the well-posedness gate requires num_centers > polynomial_cols, so
        // 3 is the smallest valid count. It is still well below the small-N
        // bump target of polynomial_cols + 4 = 6, so this exercises the
        // "explicit value is honored" path the test name advertises.
        let parsed = parse_formula("y ~ duchon(x1, centers=3) + s(x2) + s(x3) + s(x4) + s(x5)")
            .expect("parse multi-smooth formula");
        let col_map = ds.column_map();
        let mut notes = Vec::new();
        let terms = build_termspec(
            &parsed.terms,
            &ds,
            &col_map,
            &mut notes,
            &crate::resource::ResourcePolicy::default_library(),
        )
        .expect("build multi-smooth terms");
        let SmoothBasisSpec::Duchon { spec, .. } = &terms.smooth_terms[0].basis else {
            panic!("expected first smooth to be Duchon");
        };
        assert!(matches!(
            spec.center_strategy,
            CenterStrategy::FarthestPoint { num_centers: 3 }
        ));
    }

    #[test]
    fn inferred_tensor_basis_cap_uses_coordinate_support_not_duplicate_rows() {
        let mut unique_rows = Vec::new();
        for i in 0..50 {
            let theta = i as f64 / 50.0;
            for j in 0..16 {
                let h = -1.0 + 2.0 * (j as f64) / 15.0;
                let y = theta.cos() + h;
                unique_rows.push(vec![y, theta, h]);
            }
        }
        let mut repeated_rows = Vec::new();
        for _ in 0..12 {
            repeated_rows.extend(unique_rows.iter().cloned());
        }

        let unique = continuous_dataset(&["y", "theta", "h"], unique_rows);
        let repeated = continuous_dataset(&["y", "theta", "h"], repeated_rows);

        let unique_basis = inferred_tensor_basis_product(&unique);
        let repeated_basis = inferred_tensor_basis_product(&repeated);

        assert_eq!(
            unique_basis, repeated_basis,
            "duplicating existing tensor coordinates must not inflate inferred basis width"
        );
    }

    #[test]
    fn parse_bspline_boundary_conditions_and_side_selector() {
        // Non-zero anchors are rejected at parse time; the diagnostic must
        // name the side and value, which doubles as a check that the
        // `side=left` filter routes the global `anchor=` value to the
        // left endpoint (not the right).
        let mut opts = BTreeMap::new();
        opts.insert("boundary_conditions".to_string(), "anchored".to_string());
        opts.insert("side".to_string(), "left".to_string());
        opts.insert("anchor".to_string(), "2.5".to_string());
        let err = parse_bspline_boundary_conditions(&opts)
            .expect_err("non-zero left anchor must be rejected")
            .to_string();
        assert!(
            err.contains("left") && err.contains("2.5"),
            "rejection should name the affected side and value: {err}"
        );

        // Side-specific aliases (`start_bc`/`end_bc`) plus the side-specific
        // anchor key (`right_anchor`) must funnel the value onto the right
        // endpoint — verified through the rejection diagnostic.
        let mut opts = BTreeMap::new();
        opts.insert("start_bc".to_string(), "clamped".to_string());
        opts.insert("end_bc".to_string(), "zero".to_string());
        opts.insert("right_anchor".to_string(), "-1.0".to_string());
        let err = parse_bspline_boundary_conditions(&opts)
            .expect_err("non-zero right anchor must be rejected")
            .to_string();
        assert!(
            err.contains("right") && err.contains("-1"),
            "rejection should name the affected side and value: {err}"
        );

        // With anchors at zero the basis builder accepts the configuration,
        // so the same alias plumbing yields a clean `Anchored { value: 0.0 }`
        // on the right and `Clamped` on the left.
        let mut opts = BTreeMap::new();
        opts.insert("start_bc".to_string(), "clamped".to_string());
        opts.insert("end_bc".to_string(), "zero".to_string());
        let parsed = parse_bspline_boundary_conditions(&opts).expect("boundary conditions");
        assert!(matches!(
            parsed.left,
            BSplineEndpointBoundaryCondition::Clamped
        ));
        assert!(matches!(
            parsed.right,
            BSplineEndpointBoundaryCondition::Anchored { value } if value.abs() < 1e-12
        ));
    }
}
