#[cfg(test)]
mod design_assembly_constraint_tests {
use super::*;
use gam_terms::basis::{
BSplineBasisSpec, BSplineBoundaryConditions, BSplineEndpointBoundaryCondition,
BSplineIdentifiability, BSplineKnotSpec, BasisOptions, Dense, DuchonBasisSpec,
DuchonNullspaceOrder, DuchonOperatorPenaltySpec, KnotSource, MaternBasisSpec,
MaternNu, OneDimensionalBoundary, SpatialIdentifiability, SphericalSplineBasisSpec,
ThinPlateBasisSpec, build_bspline_basis_1d,
};
use gam_linalg::faer_ndarray::{FaerEigh, FaerSvd};
use gam_model_api::OuterEvalOrder;
use ndarray::{Axis, array};
use rand::RngExt as _;
use rand::SeedableRng as _;
use rand::rngs::StdRng;
fn remap_test_bspline(feature_col: usize) -> SmoothBasisSpec {
SmoothBasisSpec::BSpline1D {
feature_col,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 4,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary_conditions: BSplineBoundaryConditions::default(),
boundary: OneDimensionalBoundary::Open,
},
}
}
#[test]
fn superseded_fit_options_skip_only_rho_posterior_inference() {
let options = FitOptions {
compute_inference: true,
max_iter: 17,
..FitOptions::default()
};
let superseded = superseded_fit_options(&options);
assert!(superseded.compute_inference);
assert!(superseded.skip_rho_posterior_inference);
assert_eq!(superseded.max_iter, 17);
assert!(!options.skip_rho_posterior_inference);
}
fn structural_shape_hex(spec: &TermCollectionSpec) -> String {
let mut h = gam_runtime::warm_start::Fingerprinter::new();
spec.write_structural_shape_hash(&mut h);
h.finish_hex()
}
fn smooth_only_collection(basis: SmoothBasisSpec) -> TermCollectionSpec {
TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "s".to_string(),
basis,
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
}
}
#[test]
fn structural_shape_hash_separates_topologies_but_repeats_for_same_topology() {
let bspline = smooth_only_collection(remap_test_bspline(0));
let sphere = smooth_only_collection(SmoothBasisSpec::Sphere {
feature_cols: vec![0, 1, 2],
spec: SphericalSplineBasisSpec::default(),
});
assert_ne!(
structural_shape_hex(&bspline),
structural_shape_hex(&sphere),
"bspline and sphere topologies must key distinctly"
);
let bspline_axis1 = smooth_only_collection(remap_test_bspline(1));
assert_ne!(
structural_shape_hex(&bspline),
structural_shape_hex(&bspline_axis1),
"same basis kind on a different feature column must key distinctly"
);
let bspline_again = smooth_only_collection(remap_test_bspline(0));
assert_eq!(
structural_shape_hex(&bspline),
structural_shape_hex(&bspline_again),
"identical topology must reuse the same warm-start key"
);
}
#[test]
fn structural_kind_and_feature_cols_track_basis_identity() {
let bspline = remap_test_bspline(2);
let sphere = SmoothBasisSpec::Sphere {
feature_cols: vec![0, 1, 2],
spec: SphericalSplineBasisSpec::default(),
};
assert_ne!(bspline.structural_kind(), sphere.structural_kind());
assert_eq!(bspline.structural_kind(), "bspline_1d");
assert_eq!(sphere.structural_kind(), "sphere");
assert_eq!(bspline.structural_feature_cols(), vec![2]);
assert_eq!(sphere.structural_feature_cols(), vec![0, 1, 2]);
}
#[test]
fn remap_feature_columns_rewrites_every_index_bearing_field() {
let spec = TermCollectionSpec {
linear_terms: vec![LinearTermSpec {
name: "lin".to_string(),
feature_col: 1,
feature_cols: vec![1, 12],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
}],
random_effect_terms: vec![RandomEffectTermSpec {
name: "re".to_string(),
feature_col: 2,
drop_first_level: false,
penalized: true,
frozen_levels: Some(vec![0, 1]),
}],
smooth_terms: vec![
SmoothTermSpec {
name: "bspline".to_string(),
basis: remap_test_bspline(3),
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
SmoothTermSpec {
name: "by_variable".to_string(),
basis: SmoothBasisSpec::ByVariable {
inner: Box::new(remap_test_bspline(4)),
by_col: 5,
kind: BySmoothKind::Numeric,
by: ByVariableSpec::Numeric,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
SmoothTermSpec {
name: "by_smooth".to_string(),
basis: SmoothBasisSpec::BySmooth {
smooth: Box::new(remap_test_bspline(6)),
by_kind: ByVarKind::Numeric { feature_col: 7 },
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
SmoothTermSpec {
name: "factor_smooth".to_string(),
basis: SmoothBasisSpec::FactorSmooth {
spec: FactorSmoothSpec {
continuous_cols: vec![8],
group_col: 9,
marginal: match remap_test_bspline(0) {
SmoothBasisSpec::BSpline1D { spec, .. } => spec,
_ => unreachable!(),
},
flavour: FactorSmoothFlavour::Sz,
group_frozen_levels: Some(vec![0, 1]),
frozen_global_orthogonality: None,
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
SmoothTermSpec {
name: "thin_plate".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![10, 11],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: 1.0,
double_penalty: true,
identifiability: SpatialIdentifiability::default(),
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
],
};
let remapped: TermCollectionSpec = spec
.remap_feature_columns(|i| Ok::<usize, String>(i + 100))
.unwrap_or_else(|e| panic!("{} failed: {:?}", "remap must succeed", e));
assert_eq!(remapped.linear_terms[0].feature_col, 101);
assert_eq!(remapped.linear_terms[0].feature_cols, vec![101, 112]);
assert_eq!(remapped.random_effect_terms[0].feature_col, 102);
let collected = collect_feature_columns(&remapped);
let mut expected: Vec<usize> = (3..=11).map(|i| i + 100).collect();
expected.sort_unstable();
assert_eq!(
collected, expected,
"every smooth-basis column index must be remapped exactly once"
);
let err = spec.remap_feature_columns(|_| Err::<usize, String>("boom".to_string()));
assert_eq!(err.unwrap_err(), "boom");
}
fn collect_feature_columns(spec: &TermCollectionSpec) -> Vec<usize> {
fn walk(basis: &SmoothBasisSpec, out: &mut Vec<usize>) {
match basis {
SmoothBasisSpec::ByVariable { inner, by_col, .. }
| SmoothBasisSpec::FactorSumToZero { inner, by_col, .. } => {
out.push(*by_col);
walk(inner, out);
}
SmoothBasisSpec::BSpline1D { feature_col, .. } => out.push(*feature_col),
SmoothBasisSpec::BySmooth { smooth, by_kind } => {
match by_kind {
ByVarKind::Numeric { feature_col } | ByVarKind::Factor { feature_col, .. } => {
out.push(*feature_col)
}
}
walk(smooth, out);
}
SmoothBasisSpec::FactorSmooth { spec } => {
out.extend(spec.continuous_cols.iter().copied());
out.push(spec.group_col);
}
SmoothBasisSpec::ThinPlate { feature_cols, .. }
| SmoothBasisSpec::Sphere { feature_cols, .. }
| SmoothBasisSpec::ConstantCurvature { feature_cols, .. }
| SmoothBasisSpec::Matern { feature_cols, .. }
| SmoothBasisSpec::MeasureJet { feature_cols, .. }
| SmoothBasisSpec::Duchon { feature_cols, .. }
| SmoothBasisSpec::Pca { feature_cols, .. }
| SmoothBasisSpec::TensorBSpline { feature_cols, .. } => {
out.extend(feature_cols.iter().copied())
}
}
}
let mut out = Vec::new();
for st in &spec.smooth_terms {
walk(&st.basis, &mut out);
}
out.sort_unstable();
out
}
#[test]
fn bspline_nonzero_anchor_rejected_homogeneous_pin_baked_structurally() {
let x = Array1::linspace(0.0, 1.0, 25);
let data = x.clone().insert_axis(Axis(1));
let mk = |left, right| SmoothTermSpec {
name: "s(x, anchored_left, clamped_right)".to_string(),
basis: SmoothBasisSpec::BSpline1D {
feature_col: 0,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 4,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary_conditions: BSplineBoundaryConditions { left, right },
boundary: OneDimensionalBoundary::Open,
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let err = build_smooth_design(
data.view(),
&[mk(
BSplineEndpointBoundaryCondition::Anchored { value: 2.0 },
BSplineEndpointBoundaryCondition::Clamped,
)],
)
.expect_err("non-zero anchor must be rejected post-#1238");
match err {
gam_terms::basis::BasisError::InvalidInput(msg) => assert!(
msg.contains("anchored B-spline boundary value must be zero"),
"unexpected rejection message: {msg}"
),
other => panic!("expected InvalidInput for non-zero anchor, got {other:?}"),
}
let free = build_smooth_design(
data.view(),
&[mk(
BSplineEndpointBoundaryCondition::Free,
BSplineEndpointBoundaryCondition::Free,
)],
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "free basis builds", e));
let raw_cols = free.total_smooth_cols();
assert_eq!(raw_cols, 8);
let design = build_smooth_design(
data.view(),
&[mk(
BSplineEndpointBoundaryCondition::Anchored { value: 0.0 },
BSplineEndpointBoundaryCondition::Clamped,
)],
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "homogeneous anchored/clamped basis builds", e));
assert!(
design.linear_constraints.is_none(),
"boundary conditions are baked structurally, not emitted as linear constraints"
);
assert_eq!(
design.total_smooth_cols(),
raw_cols - 3,
"left value+slope and right slope pins drop three raw columns"
);
let BasisMetadata::BSpline1D {
knots,
identifiability_transform,
..
} = &design.terms[0].metadata
else {
panic!("expected B-spline metadata");
};
let z = identifiability_transform
.as_ref()
.unwrap_or_else(|| panic!("{} failed", "baked boundary transform recorded for replay"));
assert_eq!(z.nrows(), raw_cols);
assert_eq!(z.ncols(), design.total_smooth_cols());
let (left_value, _) = gam_terms::basis::create_basis::<Dense>(
array![0.0].view(),
KnotSource::Provided(knots.view()),
3,
BasisOptions::value(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "left endpoint value basis", e));
let (left_slope, _) = gam_terms::basis::create_basis::<Dense>(
array![0.0].view(),
KnotSource::Provided(knots.view()),
3,
BasisOptions::first_derivative(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "left endpoint derivative basis", e));
let (right_slope, _) = gam_terms::basis::create_basis::<Dense>(
array![1.0].view(),
KnotSource::Provided(knots.view()),
3,
BasisOptions::first_derivative(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "right endpoint derivative basis", e));
let baked_left_value = left_value.row(0).dot(z);
let baked_left_slope = left_slope.row(0).dot(z);
let baked_right_slope = right_slope.row(0).dot(z);
let maxabs = |v: &Array1<f64>| v.iter().fold(0.0_f64, |a, &x| a.max(x.abs()));
assert!(
maxabs(&baked_left_value) < 1e-10,
"left endpoint value must be structurally pinned to zero, got {:e}",
maxabs(&baked_left_value)
);
assert!(
maxabs(&baked_left_slope) < 1e-10,
"left endpoint slope (Hermite C1) must be structurally pinned to zero, got {:e}",
maxabs(&baked_left_slope)
);
assert!(
maxabs(&baked_right_slope) < 1e-10,
"right endpoint slope (Clamped) must be structurally pinned to zero, got {:e}",
maxabs(&baked_right_slope)
);
}
#[test]
fn bspline_boundary_conditions_follow_frozen_identifiability_transform() {
let x = Array1::linspace(0.0, 1.0, 20);
let data = x.insert_axis(Axis(1));
let raw_cols = 8;
let mut z = Array2::<f64>::zeros((raw_cols, raw_cols - 1));
for j in 0..(raw_cols - 1) {
z[[j, j]] = 1.0;
z[[raw_cols - 1, j]] = -1.0;
}
let spec = SmoothTermSpec {
name: "half-open anchored smooth".to_string(),
basis: SmoothBasisSpec::BSpline1D {
feature_col: 0,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 4,
},
double_penalty: false,
identifiability: BSplineIdentifiability::FrozenTransform {
transform: z.clone(),
},
boundary_conditions: BSplineBoundaryConditions {
left: BSplineEndpointBoundaryCondition::Anchored { value: 0.0 },
right: BSplineEndpointBoundaryCondition::Free,
},
boundary: OneDimensionalBoundary::Open,
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let design = build_smooth_design(data.view(), &[spec]).unwrap_or_else(|e| panic!("{} failed: {:?}", "boundary design", e));
assert!(
design.linear_constraints.is_none(),
"frozen-transform boundary smooth emits no linear constraints"
);
assert_eq!(design.total_smooth_cols(), raw_cols - 1);
let BasisMetadata::BSpline1D {
knots,
identifiability_transform,
..
} = &design.terms[0].metadata
else {
panic!("expected B-spline metadata");
};
let stored = identifiability_transform
.as_ref()
.unwrap_or_else(|| panic!("{} failed", "frozen transform recorded in metadata"));
assert_eq!(stored.dim(), z.dim());
for ((i, j), &v) in z.indexed_iter() {
assert!(
(stored[[i, j]] - v).abs() < 1e-12,
"stored transform must equal the frozen transform verbatim at [{i},{j}]"
);
}
let (raw_basis, _) = gam_terms::basis::create_basis::<Dense>(
data.column(0),
KnotSource::Provided(knots.view()),
3,
BasisOptions::value(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "raw basis", e));
let expected_design = raw_basis.dot(&z);
let shipped = design.term_designs[0].to_dense();
assert_eq!(shipped.dim(), expected_design.dim());
for ((i, j), &v) in expected_design.indexed_iter() {
assert!(
(shipped[[i, j]] - v).abs() < 1e-9,
"shipped design must equal raw basis · frozen transform at [{i},{j}]"
);
}
}
fn assert_spatial_derivative_width(
label: &str,
dense: &Array2<f64>,
implicit: Option<&gam_terms::basis::ImplicitDesignPsiDerivative>,
expected: usize,
) {
if let Some(op) = implicit {
assert_eq!(
op.p_out(),
expected,
"{label} implicit derivative width should match term coefficient width"
);
} else {
assert_eq!(
dense.ncols(),
expected,
"{label} dense derivative width should match term coefficient width"
);
}
}
fn numerical_rank(x: &Array2<f64>) -> usize {
let (_, s, _) = x
.svd(false, false)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "SVD should succeed in rank test", e));
let sigma_max = s.iter().copied().fold(0.0_f64, f64::max);
let tol = (x.nrows().max(x.ncols()).max(1) as f64) * f64::EPSILON * sigma_max.max(1.0);
s.iter().filter(|&&sv| sv > tol).count()
}
fn residual_norm_to_column_space(x: &Array2<f64>, y: &Array1<f64>) -> f64 {
let (u_opt, _, _) = x
.svd(true, false)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "SVD should succeed in projection residual test", e));
let u = u_opt.unwrap_or_else(|| panic!("{} failed", "left singular vectors should be present"));
let rank = numerical_rank(x);
let mut proj = Array1::<f64>::zeros(y.len());
for j in 0..rank.min(u.ncols()) {
let uj = u.column(j);
let coeff = uj.dot(y);
proj.scaled_add(coeff, &uj);
}
let resid = y - &proj;
resid.dot(&resid).sqrt()
}
fn spatial_log_kappa_bounds_from_options(
dims_per_term: &[usize],
options: &SpatialLengthScaleOptimizationOptions,
lower: bool,
) -> SpatialLogKappaCoords {
let total: usize = dims_per_term.iter().sum();
let value = if lower {
-options.max_length_scale.ln()
} else {
-options.min_length_scale.ln()
};
SpatialLogKappaCoords::new_with_dims(
Array1::<f64>::from_elem(total, value),
dims_per_term.to_vec(),
)
}
pub(super) fn two_block_exact_joint_hyper_setup(
meanspec: &TermCollectionSpec,
noisespec: &TermCollectionSpec,
kappa_options: &SpatialLengthScaleOptimizationOptions,
) -> ExactJointHyperSetup {
let mean_terms = spatial_length_scale_term_indices(meanspec);
let noise_terms = spatial_length_scale_term_indices(noisespec);
let mean_dims_per_term = spatial_dims_per_term(meanspec, &mean_terms);
let noise_dims_per_term = spatial_dims_per_term(noisespec, &noise_terms);
let mean_use_aniso = has_aniso_terms(meanspec, &mean_terms);
let noise_use_aniso = has_aniso_terms(noisespec, &noise_terms);
let mean_log_kappa = if mean_use_aniso {
SpatialLogKappaCoords::from_length_scales_aniso(meanspec, &mean_terms, kappa_options)
} else {
SpatialLogKappaCoords::from_length_scales(meanspec, &mean_terms, kappa_options)
};
let noise_log_kappa = if noise_use_aniso {
SpatialLogKappaCoords::from_length_scales_aniso(noisespec, &noise_terms, kappa_options)
} else {
SpatialLogKappaCoords::from_length_scales(noisespec, &noise_terms, kappa_options)
};
let dims_per_term = mean_log_kappa
.dims_per_term()
.iter()
.copied()
.chain(noise_log_kappa.dims_per_term().iter().copied())
.collect::<Vec<_>>();
assert_eq!(
dims_per_term,
mean_dims_per_term
.iter()
.copied()
.chain(noise_dims_per_term.iter().copied())
.collect::<Vec<_>>()
);
let log_kappa0 = SpatialLogKappaCoords::new_with_dims(
Array1::from_iter(
mean_log_kappa
.as_array()
.iter()
.chain(noise_log_kappa.as_array().iter())
.copied(),
),
dims_per_term.clone(),
);
ExactJointHyperSetup::new(
Array1::zeros(0),
Array1::zeros(0),
Array1::zeros(0),
log_kappa0,
spatial_log_kappa_bounds_from_options(&dims_per_term, kappa_options, true),
spatial_log_kappa_bounds_from_options(&dims_per_term, kappa_options, false),
)
}
fn max_abs_diff_matrix(a: &Array2<f64>, b: &Array2<f64>) -> f64 {
assert_eq!(a.dim(), b.dim());
a.iter()
.zip(b.iter())
.map(|(&x, &y)| (x - y).abs())
.fold(0.0_f64, f64::max)
}
fn assert_frozen_replay_matches_fit(
data: ArrayView2<'_, f64>,
spec: &TermCollectionSpec,
label: &str,
) {
let fit_design = build_term_collection_design(data, spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "fit-time design", e));
let frozen =
freeze_term_collection_from_design(spec, &fit_design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze term collection", e));
let replay_design = build_term_collection_design(data, &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "replay design", e));
let max_abs = max_abs_diff_matrix(
&fit_design.design.to_dense(),
&replay_design.design.to_dense(),
);
assert!(
max_abs <= 1e-10,
"{label} frozen replay changed realized design: max_abs={max_abs}"
);
}
fn dense_kronecker_pseudo_logdet_reference(
marginal_penalties: &[Array2<f64>],
lambdas: &[f64],
ridge: f64,
) -> (f64, Array1<f64>, Array2<f64>) {
let p_total: usize = marginal_penalties
.iter()
.map(|penalty| penalty.nrows())
.product();
let mut s_dense = Array2::<f64>::zeros((p_total, p_total));
for (axis, penalty) in marginal_penalties.iter().enumerate() {
let mut kron_term = Array2::<f64>::eye(1);
for (other_axis, other_penalty) in marginal_penalties.iter().enumerate() {
let factor = if axis == other_axis {
penalty.clone()
} else {
Array2::<f64>::eye(other_penalty.nrows())
};
kron_term = gam_terms::construction::kronecker_product(&kron_term, &factor);
}
s_dense.scaled_add(lambdas[axis], &kron_term);
}
if ridge > 0.0 {
for idx in 0..p_total {
s_dense[[idx, idx]] += ridge;
}
}
let (evals_dense, evecs_dense): (Array1<f64>, Array2<f64>) = s_dense
.eigh(faer::Side::Lower)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "dense Kronecker eigh", e));
let tol = 1e-12;
let positive_indices: Vec<usize> = evals_dense
.iter()
.enumerate()
.filter_map(|(idx, &value)| (value > tol).then_some(idx))
.collect();
let logdet = positive_indices
.iter()
.map(|&idx| evals_dense[idx].ln())
.sum();
let mut grad = Array1::<f64>::zeros(lambdas.len());
let mut hess = Array2::<f64>::zeros((lambdas.len(), lambdas.len()));
for (axis, penalty) in marginal_penalties.iter().enumerate() {
let mut kron_term = Array2::<f64>::eye(1);
for (other_axis, other_penalty) in marginal_penalties.iter().enumerate() {
let factor = if axis == other_axis {
penalty.clone()
} else {
Array2::<f64>::eye(other_penalty.nrows())
};
kron_term = gam_terms::construction::kronecker_product(&kron_term, &factor);
}
for &eig_idx in &positive_indices {
let eigval = evals_dense[eig_idx];
let eigvec = evecs_dense.column(eig_idx).to_owned();
let projected = kron_term.dot(&eigvec);
let ck = lambdas[axis] * eigvec.dot(&projected);
grad[axis] += ck / eigval;
hess[[axis, axis]] += ck / eigval - (ck * ck) / (eigval * eigval);
for other_axis in (axis + 1)..lambdas.len() {
let mut other_kron = Array2::<f64>::eye(1);
for (inner_axis, inner_penalty) in marginal_penalties.iter().enumerate() {
let factor = if other_axis == inner_axis {
inner_penalty.clone()
} else {
Array2::<f64>::eye(inner_penalty.nrows())
};
other_kron = gam_terms::construction::kronecker_product(&other_kron, &factor);
}
let other_projected = other_kron.dot(&eigvec);
let cl = lambdas[other_axis] * eigvec.dot(&other_projected);
let off = -(ck * cl) / (eigval * eigval);
hess[[axis, other_axis]] += off;
hess[[other_axis, axis]] += off;
}
}
}
(logdet, grad, hess)
}
fn max_abs_diff_vector(a: &Array1<f64>, b: &Array1<f64>) -> f64 {
assert_eq!(a.len(), b.len());
a.iter()
.zip(b.iter())
.map(|(&x, &y)| (x - y).abs())
.fold(0.0_f64, f64::max)
}
#[test]
fn kronecker_penalty_system_logdet_matches_dense_reference() {
let q1 = 3usize;
let q2 = 4usize;
let s1 = array![[1.0, -1.0, 0.0], [-1.0, 2.0, -1.0], [0.0, -1.0, 1.0]];
let s2 = array![
[1.0, -1.0, 0.0, 0.0],
[-1.0, 2.0, -1.0, 0.0],
[0.0, -1.0, 2.0, -1.0],
[0.0, 0.0, -1.0, 1.0]
];
let marginal_penalties = vec![s1, s2];
let lambdas = vec![2.5, 1.3];
let ridge = 0.0;
let system = KroneckerPenaltySystem::new(marginal_penalties.clone(), vec![q1, q2], false)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "KroneckerPenaltySystem", e));
let (logdet, grad, hess) = system.logdet_and_derivatives(&lambdas, ridge);
let (dense_logdet, dense_grad, dense_hess) =
dense_kronecker_pseudo_logdet_reference(&marginal_penalties, &lambdas, ridge);
assert!(
(logdet - dense_logdet).abs() < 1e-8,
"KroneckerPenaltySystem logdet mismatch: factored={} dense={}",
logdet,
dense_logdet
);
let grad_diff = max_abs_diff_vector(&grad, &dense_grad);
assert!(
grad_diff < 1e-8,
"KroneckerPenaltySystem gradient mismatch: max diff={grad_diff}"
);
let hess_diff = max_abs_diff_matrix(&hess, &dense_hess);
assert!(
hess_diff < 1e-8,
"KroneckerPenaltySystem Hessian mismatch: max diff={hess_diff}"
);
}
pub(super) fn assert_term_collection_designs_match(
left: &TermCollectionDesign,
right: &TermCollectionDesign,
label: &str,
) {
let left_design = left.design.to_dense();
let right_design = right.design.to_dense();
let design_diff = max_abs_diff_matrix(&left_design, &right_design);
assert!(
design_diff <= 1e-10,
"{label} design mismatch max_abs={design_diff}"
);
assert_eq!(
left.penalties.len(),
right.penalties.len(),
"{label} penalty count mismatch"
);
for (idx, (lp, rp)) in left
.penalties
.iter()
.zip(right.penalties.iter())
.enumerate()
{
assert_eq!(
lp.col_range, rp.col_range,
"{label} penalty {idx} col_range mismatch"
);
let penalty_diff = max_abs_diff_matrix(&lp.local, &rp.local);
assert!(
penalty_diff <= 1e-10,
"{label} penalty {idx} mismatch max_abs={penalty_diff}"
);
}
assert_eq!(
left.nullspace_dims, right.nullspace_dims,
"{label} nullspace dims mismatch"
);
assert_eq!(
left.penaltyinfo.len(),
right.penaltyinfo.len(),
"{label} penaltyinfo length mismatch"
);
for (idx, (linfo, rinfo)) in left
.penaltyinfo
.iter()
.zip(right.penaltyinfo.iter())
.enumerate()
{
assert_eq!(
linfo.termname, rinfo.termname,
"{label} penaltyinfo termname mismatch at {idx}"
);
assert_eq!(
linfo.penalty.source, rinfo.penalty.source,
"{label} penalty source mismatch at {idx}"
);
assert_eq!(
linfo.penalty.active, rinfo.penalty.active,
"{label} penalty active mismatch at {idx}"
);
assert_eq!(
linfo.penalty.effective_rank, rinfo.penalty.effective_rank,
"{label} penalty rank mismatch at {idx}"
);
assert_eq!(
linfo.penalty.nullspace_dim_hint, rinfo.penalty.nullspace_dim_hint,
"{label} penalty nullspace hint mismatch at {idx}"
);
assert!(
(linfo.penalty.normalization_scale - rinfo.penalty.normalization_scale).abs() <= 1e-10,
"{label} penalty normalization mismatch at {idx}"
);
}
match (
left.coefficient_lower_bounds.as_ref(),
right.coefficient_lower_bounds.as_ref(),
) {
(Some(lb_left), Some(lb_right)) => {
let diff = max_abs_diff_vector(lb_left, lb_right);
assert!(diff <= 1e-10, "{label} lower-bound mismatch max_abs={diff}");
}
(None, None) => {}
_ => panic!("{label} lower-bound presence mismatch"),
}
match (
left.linear_constraints.as_ref(),
right.linear_constraints.as_ref(),
) {
(Some(c_left), Some(c_right)) => {
let a_diff = max_abs_diff_matrix(&c_left.a, &c_right.a);
let b_diff = max_abs_diff_vector(&c_left.b, &c_right.b);
assert!(
a_diff <= 1e-10,
"{label} linear-constraint A mismatch max_abs={a_diff}"
);
assert!(
b_diff <= 1e-10,
"{label} linear-constraint b mismatch max_abs={b_diff}"
);
}
(None, None) => {}
_ => panic!("{label} linear-constraint presence mismatch"),
}
}
#[test]
fn smooth_design_assembles_terms_and_penalties() {
let data = array![
[0.0, 0.0, 0.2],
[0.2, 0.1, 0.4],
[0.4, 0.2, 0.6],
[0.6, 0.4, 0.7],
[0.8, 0.7, 0.9],
[1.0, 1.0, 1.1]
];
let terms = vec![
SmoothTermSpec {
name: "s_x0".to_string(),
basis: SmoothBasisSpec::BSpline1D {
feature_col: 0,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 4,
},
double_penalty: true,
identifiability: BSplineIdentifiability::default(),
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
SmoothTermSpec {
name: "tps_x1x2".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![1, 2],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: 1.0,
double_penalty: true,
identifiability: SpatialIdentifiability::default(),
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
];
let sd = build_smooth_design(data.view(), &terms).unwrap();
assert_eq!(sd.nrows(), data.nrows());
assert_eq!(sd.terms.len(), 2);
assert_eq!(sd.penalties.len(), 4);
assert_eq!(sd.nullspace_dims.len(), 4);
for bp in &sd.penalties {
assert_eq!(bp.local.nrows(), bp.block_size());
assert_eq!(bp.local.ncols(), bp.block_size());
assert!(bp.col_range.end <= sd.total_smooth_cols());
}
}
#[test]
fn shape_mapping_monotone_increasing_is_non_decreasing() {
let theta = array![-1.0, 0.5, -0.2, 0.3];
let beta =
SmoothDesign::map_term_coefficients(&theta, ShapeConstraint::MonotoneIncreasing).unwrap();
for i in 1..beta.len() {
assert!(beta[i] >= beta[i - 1]);
}
}
#[test]
fn build_smooth_design_rejectsmultiaxis_spatial_shape_constraints() {
let data = array![[0.0, 0.0], [0.5, 0.2], [1.0, 0.4], [1.5, 0.6],];
let terms = vec![SmoothTermSpec {
name: "tps_shape".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 3 },
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::default(),
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::MonotoneIncreasing,
joint_null_rotation: None,
}];
let err = build_smooth_design(data.view(), &terms).expect_err("shape should be rejected");
match err {
BasisError::InvalidInput(msg) => {
assert!(msg.contains("requires exactly 1 feature axis"));
}
other => panic!("unexpected error: {other:?}"),
}
}
#[test]
fn build_smooth_design_accepts_monotone_thin_plate_1dwith_linear_constraints() {
let data = array![[0.0], [0.15], [0.35], [0.5], [0.65], [0.85], [1.0]];
let terms = vec![SmoothTermSpec {
name: "mono_tps".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::default(),
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::MonotoneIncreasing,
joint_null_rotation: None,
}];
let sd = build_smooth_design(data.view(), &terms).unwrap_or_else(|e| panic!("{} failed: {:?}", "shape-constrained thin-plate", e));
assert!(sd.coefficient_lower_bounds.is_none());
let lin = sd
.linear_constraints
.as_ref()
.unwrap_or_else(|| panic!("{} failed", "linear constraints should be generated"));
assert!(lin.a.nrows() > 0);
assert_eq!(lin.a.ncols(), sd.total_smooth_cols());
assert_eq!(lin.b.len(), lin.a.nrows());
}
#[test]
fn build_smooth_design_auto_promotes_thin_plate_below_canonical_polynomial_dimension() {
let data = array![
[0.0, 0.0, 0.0],
[0.2, 0.1, 0.3],
[0.4, 0.3, 0.5],
[0.7, 0.6, 0.8],
];
let terms = vec![SmoothTermSpec {
name: "thinplate(pc1, pc2, pc3)".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1, 2],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 3 },
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::default(),
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}];
let sd = build_smooth_design(data.view(), &terms)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "auto-promotion to Duchon should succeed at infeasible canonical (d, k)", e));
let metadata = &sd.terms.first().expect("at least one smooth term").metadata;
assert!(
matches!(metadata, BasisMetadata::Duchon { .. }),
"expected Duchon metadata after auto-promotion, got {metadata:?}"
);
}
#[test]
fn freeze_term_collection_handles_thin_plate_auto_promotion_to_duchon() {
let mut rng = StdRng::seed_from_u64(20260504);
let n = 200usize;
let mut data = Array2::<f64>::zeros((n, 5));
for i in 0..n {
for j in 0..5 {
data[[i, j]] = rng.random_range(-1.0..1.0);
}
}
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "thinplate(pc1, pc2, pc3, pc4, pc5)".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1, 2, 3, 4],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 10 },
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::default(),
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let fit_design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "fit-time design", e));
let metadata = &fit_design
.smooth
.terms
.first()
.unwrap_or_else(|| panic!("{} failed", "at least one smooth term"))
.metadata;
assert!(
matches!(metadata, BasisMetadata::Duchon { .. }),
"expected auto-promotion to Duchon, got {metadata:?}"
);
let frozen = freeze_term_collection_from_design(&spec, &fit_design).unwrap_or_else(|e| panic!("{} failed: {:?}",
"freeze must succeed across the auto-promoted (ThinPlate spec, Duchon metadata) pair", e
));
assert!(
matches!(frozen.smooth_terms[0].basis, SmoothBasisSpec::Duchon { .. }),
"frozen spec should reflect the auto-promotion as a Duchon variant"
);
let replay_design = build_term_collection_design(data.view(), &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "replay design", e));
let max_abs = fit_design
.design
.to_dense()
.iter()
.zip(replay_design.design.to_dense().iter())
.map(|(&a, &b)| (a - b).abs())
.fold(0.0_f64, f64::max);
assert!(
max_abs <= 1e-10,
"auto-promoted frozen replay changed realized design: max_abs={max_abs}"
);
}
#[test]
fn build_smooth_design_accepts_monotone_matern_1dwith_linear_constraints() {
let data = array![[0.0], [0.2], [0.4], [0.6], [0.8], [1.0]];
let terms = vec![SmoothTermSpec {
name: "mono_matern".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: 0.7,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: false,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::MonotoneIncreasing,
joint_null_rotation: None,
}];
let sd = build_smooth_design(data.view(), &terms).unwrap_or_else(|e| panic!("{} failed: {:?}", "shape-constrained Matérn", e));
assert!(sd.coefficient_lower_bounds.is_none());
let lin = sd
.linear_constraints
.as_ref()
.unwrap_or_else(|| panic!("{} failed", "linear constraints should be generated"));
assert!(lin.a.nrows() > 0);
assert_eq!(lin.a.ncols(), sd.total_smooth_cols());
assert_eq!(lin.b.len(), lin.a.nrows());
}
#[test]
fn build_smooth_design_accepts_monotone_duchon_1dwith_linear_constraints() {
let data = array![[0.0], [0.2], [0.4], [0.6], [0.8], [1.0]];
let terms = vec![SmoothTermSpec {
name: "mono_duchon".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: Some(0.9),
power: 5.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::OrthogonalToParametric,
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::MonotoneIncreasing,
joint_null_rotation: None,
}];
let sd = build_smooth_design(data.view(), &terms).unwrap_or_else(|e| panic!("{} failed: {:?}", "shape-constrained Duchon", e));
assert!(sd.coefficient_lower_bounds.is_none());
let lin = sd
.linear_constraints
.as_ref()
.unwrap_or_else(|| panic!("{} failed", "linear constraints should be generated"));
assert!(lin.a.nrows() > 0);
assert_eq!(lin.a.ncols(), sd.total_smooth_cols());
assert_eq!(lin.b.len(), lin.a.nrows());
}
#[test]
fn build_smooth_design_accepts_monotone_bsplinewith_bounds() {
let data = array![[0.0], [0.25], [0.5], [0.75], [1.0]];
let terms = vec![SmoothTermSpec {
name: "mono_bs".to_string(),
basis: SmoothBasisSpec::BSpline1D {
feature_col: 0,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 3,
},
double_penalty: false,
identifiability: BSplineIdentifiability::default(),
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
},
},
shape: ShapeConstraint::MonotoneIncreasing,
joint_null_rotation: None,
}];
let sd = build_smooth_design(data.view(), &terms).unwrap_or_else(|e| panic!("{} failed: {:?}", "shape-constrained bspline", e));
let lb = sd
.coefficient_lower_bounds
.as_ref()
.unwrap_or_else(|| panic!("{} failed", "lower bounds should be generated"));
assert_eq!(lb.len(), sd.total_smooth_cols());
assert!(lb[0].is_infinite() && lb[0].is_sign_negative());
for j in 1..lb.len() {
assert_eq!(lb[j], 0.0);
}
}
#[test]
fn term_collection_design_combines_linear_and_smooth() {
let data = array![
[0.0, 0.0, 0.2],
[0.2, 0.1, 0.4],
[0.4, 0.2, 0.6],
[0.6, 0.4, 0.7],
[0.8, 0.7, 0.9],
[1.0, 1.0, 1.1]
];
let spec = TermCollectionSpec {
linear_terms: vec![LinearTermSpec {
name: "lin_x0".to_string(),
feature_col: 0,
feature_cols: vec![0],
categorical_levels: vec![],
double_penalty: true,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
}],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "tps_x1x2".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![1, 2],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: 1.0,
double_penalty: true,
identifiability: SpatialIdentifiability::default(),
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap();
let design_dense = design.design.to_dense();
assert_eq!(design.design.nrows(), data.nrows());
assert_eq!(design.intercept_range, 0..1);
assert!(
design_dense
.column(design.intercept_range.start)
.iter()
.all(|&v: &f64| (v - 1.0).abs() < 1e-12)
);
assert!(design.design.ncols() >= 2);
assert_eq!(design.linear_ranges.len(), 1);
assert_eq!(design.random_effect_ranges.len(), 0);
assert_eq!(design.penalties.len(), 3); assert_eq!(design.nullspace_dims.len(), 3);
}
#[test]
fn spatial_smooth_columns_do_not_duplicate_global_intercept() {
let data = array![
[0.0, 0.0],
[0.2, 0.1],
[0.4, 0.3],
[0.6, 0.6],
[0.8, 0.7],
[1.0, 1.0],
];
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "tps_xy".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::default(),
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap();
let design_dense = design.design.to_dense();
let smooth_start = 1usize;
let smooth_end = smooth_start + design.smooth.total_smooth_cols();
for col in smooth_start..smooth_end {
let is_all_ones = design_dense
.column(col)
.iter()
.all(|&v: &f64| (v - 1.0).abs() < 1e-12);
assert!(
!is_all_ones,
"smooth column {col} unexpectedly duplicated intercept"
);
}
}
#[test]
fn spatial_smooth_drops_matching_linear_trend_columns() {
let data = array![
[0.0, 0.1],
[0.2, 0.0],
[0.3, 0.4],
[0.5, 0.2],
[0.7, 0.9],
[1.0, 0.8],
[1.2, 1.1],
[1.4, 1.3],
];
let spec = TermCollectionSpec {
linear_terms: vec![LinearTermSpec {
name: "lin_x0".to_string(),
feature_col: 0,
feature_cols: vec![0],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
}],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "tps_xy".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::default(),
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap();
assert_eq!(design.smooth.total_smooth_cols(), 2);
let dense = design.design.to_dense_cow();
let lin_col = design.linear_ranges[0].1.start;
let linvalues = dense.column(lin_col).to_owned();
let smooth_start = 1 + spec.linear_terms.len();
let smooth_end = smooth_start + design.smooth.total_smooth_cols();
for col in smooth_start..smooth_end {
let same_as_linear = dense
.column(col)
.iter()
.zip(linvalues.iter())
.all(|(&a, &b)| (a - b).abs() < 1e-12);
assert!(
!same_as_linear,
"smooth column {col} unexpectedly duplicated linear term column"
);
}
}
#[test]
fn spatial_option5_is_orthogonal_to_parametric_block() {
let data = array![
[0.0, 0.1],
[0.2, 0.0],
[0.3, 0.4],
[0.5, 0.2],
[0.7, 0.9],
[1.0, 0.8],
];
let spec = TermCollectionSpec {
linear_terms: vec![LinearTermSpec {
name: "lin_x0".to_string(),
feature_col: 0,
feature_cols: vec![0],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
}],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "tps_xy".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::OrthogonalToParametric,
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap();
let design_dense = design.design.to_dense();
let n = data.nrows();
let mut c = Array2::<f64>::zeros((n, 2));
c.column_mut(0).fill(1.0);
c.column_mut(1).assign(&data.column(0));
let smooth_start = 1 + spec.linear_terms.len();
let b = design_dense
.slice(s![
..,
smooth_start..(smooth_start + design.smooth.total_smooth_cols())
])
.to_owned();
let cross = b.t().dot(&c);
let num = cross.iter().map(|v| v * v).sum::<f64>().sqrt();
let b_norm = b.iter().map(|v| v * v).sum::<f64>().sqrt();
let c_norm = c.iter().map(|v| v * v).sum::<f64>().sqrt();
let rel = num / (b_norm * c_norm).max(1e-300);
assert!(
rel <= 1e-10,
"smooth residual against model-owned parametric block too large: {rel}"
);
}
#[test]
fn thin_plate_default_identifiability_centers_against_intercept_only_without_linear_terms() {
let data = array![
[-1.9, -1.2],
[-1.3, -0.7],
[-0.8, -0.4],
[-0.2, 0.1],
[0.0, 0.3],
[0.4, 0.5],
[0.9, 0.8],
[1.4, 1.1],
[1.9, 1.5],
[2.3, 1.8],
];
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: (0..2)
.map(|feature| SmoothTermSpec {
name: format!("tps_x{feature}"),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![feature],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::EqualMass { num_centers: 4 },
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::OrthogonalToParametric,
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
})
.collect(),
};
let design = build_term_collection_design(data.view(), &spec).unwrap();
let design_dense = design.design.to_dense();
let smooth_start = 1 + spec.linear_terms.len();
let intercept = Array2::<f64>::ones((data.nrows(), 1));
for (term_idx, term) in design.smooth.terms.iter().enumerate() {
let block = design_dense
.slice(s![
..,
(smooth_start + term.coeff_range.start)..(smooth_start + term.coeff_range.end)
])
.to_owned();
let cross = block.t().dot(&intercept);
let num = cross.iter().map(|v| v * v).sum::<f64>().sqrt();
let block_norm = block.iter().map(|v| v * v).sum::<f64>().sqrt();
let intercept_norm = intercept.iter().map(|v| v * v).sum::<f64>().sqrt();
let rel = num / (block_norm * intercept_norm).max(1e-300);
assert!(
rel <= 1e-10,
"ThinPlate term {term_idx} should be centered against the intercept (no linear terms in formula); got rel={rel:.3e}"
);
}
}
#[test]
fn spatial_option5_does_not_overconstrain_on_nonoverlapping_linear_terms() {
let n = 40usize;
let p = 16usize;
let mut data = Array2::<f64>::zeros((n, p));
for i in 0..n {
for j in 0..p {
data[[i, j]] = (i as f64) * 0.03 + (j as f64) * 0.11 + ((i * (j + 1)) as f64) * 1e-3;
}
}
let spec = TermCollectionSpec {
linear_terms: (5..16)
.map(|j| LinearTermSpec {
name: format!("pc{j}"),
feature_col: j,
feature_cols: vec![j],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
})
.collect(),
random_effect_terms: vec![],
smooth_terms: vec![
SmoothTermSpec {
name: "tps_pc1".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![1],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 12 },
length_scale: 1.0,
double_penalty: true,
identifiability: SpatialIdentifiability::OrthogonalToParametric,
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
SmoothTermSpec {
name: "tps_pc2".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![2],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 12 },
length_scale: 1.0,
double_penalty: true,
identifiability: SpatialIdentifiability::OrthogonalToParametric,
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
],
};
let out = build_term_collection_design(data.view(), &spec);
assert!(
out.is_ok(),
"term-local Option 5 should not over-constrain non-overlapping smooth/linear terms: {:?}",
out.err()
);
}
#[test]
fn overlapping_linear_term_residualizes_bspline_smooth() {
let data = array![
[0.0],
[0.1],
[0.2],
[0.3],
[0.4],
[0.5],
[0.6],
[0.7],
[0.8],
[0.9],
[1.0],
];
let spec = TermCollectionSpec {
linear_terms: vec![LinearTermSpec {
name: "x".to_string(),
feature_col: 0,
feature_cols: vec![0],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
}],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "s_x".to_string(),
basis: SmoothBasisSpec::BSpline1D {
feature_col: 0,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 4,
},
double_penalty: false,
identifiability: BSplineIdentifiability::default(),
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "bspline design", e));
let mut c = Array2::<f64>::zeros((data.nrows(), 2));
c.column_mut(0).fill(1.0);
c.column_mut(1).assign(&data.column(0));
let rel = orthogonality_relative_residual_for_design(&design.smooth.term_designs[0], c.view())
.unwrap_or_else(|e| panic!("{} failed: {:?}", "orthogonality residual", e));
assert!(
rel <= 1e-10,
"B-spline smooth should be orthogonal to [1, x] when linear(x) is present; rel={rel}"
);
}
#[test]
fn standalone_tps_keeps_centered_linear_nullspace() {
let data = array![[-1.5], [-0.7], [0.2], [0.8], [1.6]];
let centers = array![[-1.5], [0.2], [1.6]];
let smooth = SmoothTermSpec {
name: "s_x".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::UserProvided(centers),
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::OrthogonalToParametric,
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![smooth],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "tps design", e));
assert_eq!(design.smooth.term_designs[0].ncols(), 2);
assert_eq!(design.smooth.nullspace_dims, vec![1]);
let intercept = Array2::<f64>::ones((data.nrows(), 1));
let rel = orthogonality_relative_residual_for_design(
&design.smooth.term_designs[0],
intercept.view(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "intercept residual", e));
assert!(
rel <= 1e-10,
"standalone TPS should be centered against the intercept while retaining its linear nullspace; rel={rel}"
);
}
#[test]
fn spatial_parametric_ownership_projects_only_explicit_linear_axes() {
let term = SmoothTermSpec {
name: "s_xy".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::EqualMass { num_centers: 4 },
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::OrthogonalToParametric,
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let linear_terms = vec![LinearTermSpec {
name: "x0".to_string(),
feature_col: 0,
feature_cols: vec![0],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
}];
assert_eq!(
smooth_intrinsic_parametric_feature_cols(&linear_terms, &term),
vec![0],
"a linear term on x0 should not claim the smooth's x1 nullspace"
);
}
#[test]
fn hierarchical_smooth_ownership_is_order_independent_for_bspline_and_duchon() {
let data = array![
[0.00, 0.00],
[0.10, 0.15],
[0.18, 0.30],
[0.27, 0.10],
[0.35, 0.55],
[0.46, 0.25],
[0.54, 0.70],
[0.63, 0.40],
[0.72, 0.85],
[0.81, 0.60],
[0.90, 0.95],
[1.00, 0.75],
];
let bspline_term = SmoothTermSpec {
name: "s_x".to_string(),
basis: SmoothBasisSpec::BSpline1D {
feature_col: 0,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 5,
},
double_penalty: false,
identifiability: BSplineIdentifiability::default(),
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let duchon_term = SmoothTermSpec {
name: "duchon_xy".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0, 1],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
length_scale: Some(1.0),
power: 5.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let spec_a = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![duchon_term.clone(), bspline_term.clone()],
};
let spec_b = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![bspline_term, duchon_term],
};
let design_a = build_term_collection_design(data.view(), &spec_a).unwrap_or_else(|e| panic!("{} failed: {:?}", "design a", e));
let design_b = build_term_collection_design(data.view(), &spec_b).unwrap_or_else(|e| panic!("{} failed: {:?}", "design b", e));
for design in [&design_a, &design_b] {
let owner_idx = design
.smooth
.terms
.iter()
.position(|term| term.name == "s_x")
.unwrap_or_else(|| panic!("{} failed", "owner term"));
let target_idx = design
.smooth
.terms
.iter()
.position(|term| term.name == "duchon_xy")
.unwrap_or_else(|| panic!("{} failed", "target term"));
let owner_dense = design.smooth.term_designs[owner_idx].to_dense();
let rel = orthogonality_relative_residual_for_design(
&design.smooth.term_designs[target_idx],
owner_dense.view(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "orthogonality residual", e));
assert!(
rel <= 1e-10,
"multivariate Duchon term should be residualized against owned 1D spline space; rel={rel}"
);
}
let duchon_a_idx = design_a
.smooth
.terms
.iter()
.position(|term| term.name == "duchon_xy")
.unwrap_or_else(|| panic!("{} failed", "duchon in design a"));
let duchon_b_idx = design_b
.smooth
.terms
.iter()
.position(|term| term.name == "duchon_xy")
.unwrap_or_else(|| panic!("{} failed", "duchon in design b"));
let duchon_a = design_a.smooth.term_designs[duchon_a_idx].to_dense();
let duchon_b = design_b.smooth.term_designs[duchon_b_idx].to_dense();
assert_eq!(duchon_a.dim(), duchon_b.dim());
let max_abs = duchon_a
.iter()
.zip(duchon_b.iter())
.map(|(lhs, rhs)| (lhs - rhs).abs())
.fold(0.0_f64, f64::max);
assert!(
max_abs <= 1e-10,
"hierarchical ownership should not depend on user term order; max_abs={max_abs}"
);
}
#[test]
fn freeze_roundtrip_preserves_hierarchical_smooth_transforms() {
let data = array![
[0.00, 0.00],
[0.10, 0.15],
[0.18, 0.30],
[0.27, 0.10],
[0.35, 0.55],
[0.46, 0.25],
[0.54, 0.70],
[0.63, 0.40],
[0.72, 0.85],
[0.81, 0.60],
[0.90, 0.95],
[1.00, 0.75],
];
let spec = TermCollectionSpec {
linear_terms: vec![LinearTermSpec {
name: "x".to_string(),
feature_col: 0,
feature_cols: vec![0],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
}],
random_effect_terms: vec![],
smooth_terms: vec![
SmoothTermSpec {
name: "duchon_xy".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0, 1],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
length_scale: Some(1.0),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
SmoothTermSpec {
name: "s_x".to_string(),
basis: SmoothBasisSpec::BSpline1D {
feature_col: 0,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 5,
},
double_penalty: false,
identifiability: BSplineIdentifiability::default(),
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "fit-time design", e));
let frozen =
freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze hierarchical design", e));
let replay = build_term_collection_design(data.view(), &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "replay design", e));
let dense_fit = design.design.to_dense();
let dense_replay = replay.design.to_dense();
assert_eq!(dense_fit.dim(), dense_replay.dim());
let max_abs = dense_fit
.iter()
.zip(dense_replay.iter())
.map(|(lhs, rhs)| (lhs - rhs).abs())
.fold(0.0_f64, f64::max);
assert!(
max_abs <= 1e-10,
"frozen hierarchical transforms should replay exactly on the training data; max_abs={max_abs}"
);
}
#[test]
fn spatial_option5_preserves_lazy_thin_plate_terms_at_large_scale() {
let n = 200_000usize;
let k = 256usize;
let mut data = Array2::<f64>::zeros((n, 1));
let mut centers = Array2::<f64>::zeros((k, 1));
for i in 0..n {
data[[i, 0]] = i as f64 / (n - 1) as f64;
}
for j in 0..k {
centers[[j, 0]] = j as f64 / (k - 1) as f64;
}
let spec = TermCollectionSpec {
linear_terms: vec![LinearTermSpec {
name: "x".to_string(),
feature_col: 0,
feature_cols: vec![0],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
}],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "tps_x".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::UserProvided(centers),
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::OrthogonalToParametric,
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "large option-5 design", e));
assert!(matches!(
&design.smooth.term_designs[0],
DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::Lazy(_))
));
let mut c = Array2::<f64>::zeros((n, 2));
c.column_mut(0).fill(1.0);
c.column_mut(1).assign(&data.column(0));
let rel = orthogonality_relative_residual_for_design(&design.smooth.term_designs[0], c.view())
.unwrap_or_else(|e| panic!("{} failed: {:?}", "orthogonality residual", e));
assert!(rel <= 1e-8, "lazy option-5 residual too large: {rel}");
}
#[test]
fn spatial_frozen_transform_rebuild_is_exact_on_trainingrows() {
let data = array![
[0.0, 0.1],
[0.2, 0.0],
[0.3, 0.4],
[0.5, 0.2],
[0.7, 0.9],
[1.0, 0.8],
];
let fitspec = TermCollectionSpec {
linear_terms: vec![LinearTermSpec {
name: "lin_x0".to_string(),
feature_col: 0,
feature_cols: vec![0],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
}],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "tps_xy".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: 1.0,
double_penalty: false,
identifiability: SpatialIdentifiability::OrthogonalToParametric,
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let fit_design = build_term_collection_design(data.view(), &fitspec).unwrap();
let term_meta = &fit_design.smooth.terms[0].metadata;
let (centers, length_scale, z) = match term_meta {
BasisMetadata::ThinPlate {
centers,
length_scale,
identifiability_transform,
..
} => (
centers.clone(),
*length_scale,
identifiability_transform
.clone()
.expect("fit-time Option 5 should store transform"),
),
other => panic!("unexpected metadata variant: {other:?}"),
};
let frozenspec = TermCollectionSpec {
linear_terms: fitspec.linear_terms.clone(),
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "tps_xy".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::UserProvided(centers),
length_scale,
double_penalty: false,
identifiability: SpatialIdentifiability::FrozenTransform { transform: z },
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let frozen_design = build_term_collection_design(data.view(), &frozenspec).unwrap();
assert_eq!(
fit_design.smooth.term_designs.len(),
frozen_design.smooth.term_designs.len(),
"frozen transform rebuild term count mismatch"
);
let max_abs = fit_design
.smooth
.term_designs
.iter()
.zip(frozen_design.smooth.term_designs.iter())
.flat_map(|(a, b)| {
let a_dense = a.to_dense();
let b_dense = b.to_dense();
assert_eq!(a_dense.dim(), b_dense.dim());
a_dense
.iter()
.zip(b_dense.iter())
.map(|(&x, &y)| (x - y).abs())
.collect::<Vec<_>>()
})
.fold(0.0_f64, f64::max);
assert!(
max_abs <= 1e-12,
"frozen transform rebuild mismatch max_abs={max_abs}"
);
}
#[test]
fn frozen_spatial_replay_preserves_standardized_length_scale_compensation() {
let n = 16usize;
let mut data = Array2::<f64>::zeros((n, 2));
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = 0.07 * i as f64 + 0.02 * (3.0 * t).sin();
data[[i, 1]] = 4.0 * t + 0.35 * (5.0 * t).cos();
}
let tps_spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "tps_xy".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
length_scale: 1.3,
double_penalty: true,
identifiability: SpatialIdentifiability::OrthogonalToParametric,
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
assert_frozen_replay_matches_fit(data.view(), &tps_spec, "thin-plate");
let matern_spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "matern_xy".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
length_scale: 1.1,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
assert_frozen_replay_matches_fit(data.view(), &matern_spec, "matern");
let duchon_spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "duchon_xy".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0, 1],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
length_scale: Some(1.4),
power: 5.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::OrthogonalToParametric,
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
assert_frozen_replay_matches_fit(data.view(), &duchon_spec, "duchon");
}
#[test]
fn term_collection_design_adds_random_effect_dummy_blockwithridge() {
let data = array![
[0.1, 0.0],
[0.2, 1.0],
[0.3, 0.0],
[0.4, 2.0],
[0.5, 1.0],
[0.6, 2.0],
];
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![RandomEffectTermSpec {
name: "id".to_string(),
feature_col: 1,
drop_first_level: false,
penalized: true,
frozen_levels: None,
}],
smooth_terms: vec![],
};
let design = build_term_collection_design(data.view(), &spec).unwrap();
assert_eq!(design.intercept_range, 0..1);
assert_eq!(design.design.ncols(), 4);
assert_eq!(design.random_effect_ranges.len(), 1);
assert_eq!(design.penalties.len(), 1);
assert_eq!(design.nullspace_dims, vec![0]);
let (_, range) = &design.random_effect_ranges[0];
let dense = design.design.to_dense_cow();
for i in 0..dense.nrows() {
let row_sum: f64 = dense.slice(s![i, range.clone()]).sum();
assert!((row_sum - 1.0).abs() < 1e-12);
}
}
#[test]
fn matern_smooth_buildswith_double_penalty_in_high_dim() {
let n = 12usize;
let d = 10usize;
let mut data = Array2::<f64>::zeros((n, d));
for i in 0..n {
for j in 0..d {
data[[i, j]] = (i as f64) * 0.1 + (j as f64) * 0.03;
}
}
let terms = vec![SmoothTermSpec {
name: "matern_x".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: (0..d).collect(),
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 5 },
length_scale: 0.75,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}];
let sd = build_smooth_design(data.view(), &terms).unwrap();
assert_eq!(sd.nrows(), n);
assert_eq!(sd.terms.len(), 1);
assert_eq!(sd.penalties.len(), 3);
assert_eq!(sd.nullspace_dims.len(), 3);
}
#[test]
fn duchon_linear_nullspace_builds_and_reports_nullspace_dim() {
let n = 20usize;
let d = 10usize;
let mut data = Array2::<f64>::zeros((n, d));
for i in 0..n {
for j in 0..d {
let mut key = (i as u64).wrapping_mul(0x9E3779B97F4A7C15);
key ^= (j as u64).wrapping_mul(0xBF58476D1CE4E5B9);
key = (key ^ (key >> 30)).wrapping_mul(0xBF58476D1CE4E5B9);
key = (key ^ (key >> 27)).wrapping_mul(0x94D049BB133111EB);
let v = ((key ^ (key >> 31)) as f64) / (u64::MAX as f64);
data[[i, j]] = v;
}
}
let terms = vec![SmoothTermSpec {
name: "duchon_x".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: (0..d).collect(),
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 12 },
length_scale: Some(0.9),
power: 5.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}];
let sd = build_smooth_design(data.view(), &terms).unwrap();
assert_eq!(sd.nrows(), n);
assert_eq!(sd.terms.len(), 1);
assert_eq!(sd.penalties.len(), 4);
assert_eq!(sd.nullspace_dims.len(), 4);
}
#[test]
fn joint_duchon_orderzero_raw_smooth_build_preserves_unconstrained_basis() {
let n = 12usize;
let d = 4usize;
let mut data = Array2::<f64>::zeros((n, d));
for i in 0..n {
for j in 0..d {
data[[i, j]] = (i as f64) * 0.13 + (j as f64) * 0.17;
}
}
let terms = vec![SmoothTermSpec {
name: "duchon_joint".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: (0..d).collect(),
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: Some(1.0),
power: 3.0,
nullspace_order: DuchonNullspaceOrder::Zero,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}];
let sd = build_smooth_design(data.view(), &terms).unwrap_or_else(|e| panic!("{} failed: {:?}", "joint duchon build", e));
assert_eq!(sd.total_smooth_cols(), 4);
match &sd.terms[0].metadata {
BasisMetadata::Duchon {
identifiability_transform,
..
} => {
assert!(
identifiability_transform.is_some(),
"raw smooth build should freeze Duchon orthogonality once the basis is built"
);
}
other => panic!("expected Duchon metadata, got {other:?}"),
}
}
#[test]
fn term_collection_joint_duchon_carries_frozen_transform_into_metadata() {
let n = 12usize;
let d = 4usize;
let mut data = Array2::<f64>::zeros((n, d));
for i in 0..n {
for j in 0..d {
data[[i, j]] = (i as f64) * 0.13 + (j as f64) * 0.17;
}
}
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "duchon_joint".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: (0..d).collect(),
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: Some(1.0),
power: 3.0,
nullspace_order: DuchonNullspaceOrder::Zero,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "term collection design", e));
let term = &design.smooth.terms[0];
assert_eq!(term.coeff_range.len(), 3);
match &term.metadata {
BasisMetadata::Duchon {
identifiability_transform,
..
} => {
let z = identifiability_transform
.as_ref()
.unwrap_or_else(|| panic!("{} failed", "term collection should store frozen Duchon transform"));
assert_eq!(z.nrows(), 4);
assert_eq!(z.ncols(), 3);
}
other => panic!("expected Duchon metadata, got {other:?}"),
}
}
#[test]
fn frozen_joint_maternspec_rebuild_keeps_adaptive_cache_in_sync() {
let n = 12usize;
let d = 2usize;
let mut data = Array2::<f64>::zeros((n, d));
for i in 0..n {
data[[i, 0]] = i as f64 * 0.13;
data[[i, 1]] = (i as f64 * 0.17).sin();
}
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "matern_joint".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: (0..d).collect(),
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
length_scale: 1.0,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "base design", e));
let frozen = freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze spec", e));
let rebuilt = build_term_collection_design(data.view(), &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "rebuilt design", e));
let caches =
extract_spatial_operator_runtime_caches(&frozen, &rebuilt).unwrap_or_else(|e| panic!("{} failed: {:?}", "adaptive caches", e));
assert_eq!(caches.len(), 1);
assert_eq!(caches[0].termname, "matern_joint");
assert_eq!(rebuilt.smooth.terms.len(), 1);
assert!(!rebuilt.smooth.terms[0].coeff_range.is_empty());
}
#[test]
fn tensor_bspline_term_builds_te_style_design_and_penalties() {
let n = 10usize;
let mut data = Array2::<f64>::zeros((n, 2));
for i in 0..n {
data[[i, 0]] = i as f64 / (n as f64 - 1.0);
data[[i, 1]] = (i as f64 / (n as f64 - 1.0)).powi(2);
}
let spec_x = BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 3,
},
double_penalty: false,
identifiability: BSplineIdentifiability::default(),
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
};
let spec_y = BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 2,
},
double_penalty: false,
identifiability: BSplineIdentifiability::default(),
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
};
let terms = vec![SmoothTermSpec {
name: "te_xy".to_string(),
basis: SmoothBasisSpec::TensorBSpline {
feature_cols: vec![0, 1],
spec: TensorBSplineSpec {
periods: Vec::new(),
marginalspecs: vec![spec_x, spec_y],
double_penalty: true,
identifiability: TensorBSplineIdentifiability::default(),
penalty_decomposition: Default::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}];
let sd = build_smooth_design(data.view(), &terms).unwrap();
assert_eq!(sd.nrows(), n);
assert_eq!(sd.terms.len(), 1);
assert_eq!(sd.penalties.len(), 3);
assert_eq!(sd.nullspace_dims.len(), 3);
assert!(
sd.penalties
.iter()
.all(|bp| bp.local.nrows() == bp.block_size())
);
assert!(
sd.penalties
.iter()
.all(|bp| bp.col_range.end <= sd.total_smooth_cols())
);
}
#[test]
fn tensor_binary_margin_is_penalized_factor_smooth_not_unidentified_raw_tensor() {
let n = 18usize;
let mut data = Array2::<f64>::zeros((n, 2));
for i in 0..n {
data[[i, 0]] = i as f64 / (n as f64 - 1.0);
data[[i, 1]] = if i % 2 == 0 { 0.0 } else { 1.0 };
}
let age_margin = BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 1,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
};
let binary_margin = BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 1,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
};
let terms = vec![SmoothTermSpec {
name: "te_age_binary".to_string(),
basis: SmoothBasisSpec::TensorBSpline {
feature_cols: vec![0, 1],
spec: TensorBSplineSpec {
periods: Vec::new(),
marginalspecs: vec![age_margin, binary_margin],
double_penalty: false,
identifiability: TensorBSplineIdentifiability::None,
penalty_decomposition: Default::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}];
let design = build_smooth_design(data.view(), &terms).unwrap_or_else(|e| panic!("{} failed: {:?}", "binary-margin tensor", e));
let kron = design.terms[0]
.kronecker_factored
.as_ref()
.unwrap_or_else(|| panic!("{} failed", "tensor term should preserve metadata"));
assert_eq!(kron.marginal_dims, vec![5, 5]);
assert_eq!(
numerical_rank(&kron.marginal_designs[1]),
2,
"a binary tensor margin has two data-supported levels even when its raw spline margin has five columns"
);
assert_eq!(
numerical_rank(&kron.marginal_penalties[1]),
3,
"the binary margin's second-difference roughness penalty must shrink the three unsupported spline-range directions"
);
assert_eq!(design.penalties.len(), 2);
assert!(design.penalties.iter().all(|penalty| {
penalty.local.nrows() == 25
&& penalty.local.ncols() == 25
&& penalty.local.iter().all(|value| value.is_finite())
}));
}
#[test]
fn centered_tensor_penalties_canonicalize_in_transformed_basis_width() {
let n = 16usize;
let mut data = Array2::<f64>::zeros((n, 2));
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
data[[i, 1]] = 0.5 + 0.25 * (2.0 * std::f64::consts::PI * t).sin();
}
let tensor_term = SmoothTermSpec {
name: "te_centered".to_string(),
basis: SmoothBasisSpec::TensorBSpline {
feature_cols: vec![0, 1],
spec: TensorBSplineSpec {
marginalspecs: vec![
BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 3,
},
double_penalty: false,
identifiability: BSplineIdentifiability::default(),
boundary_conditions: Default::default(),
boundary: OneDimensionalBoundary::Open,
},
BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 2,
},
double_penalty: false,
identifiability: BSplineIdentifiability::default(),
boundary_conditions: Default::default(),
boundary: OneDimensionalBoundary::Open,
},
],
periods: Vec::new(),
double_penalty: false,
identifiability: TensorBSplineIdentifiability::default(),
penalty_decomposition: Default::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![tensor_term],
};
let design = build_term_collection_design(data.view(), &spec).unwrap();
let penalty_specs = design
.penalties
.iter()
.map(PenaltySpec::from_blockwise_ref)
.collect::<Vec<_>>();
let (canonical, _) = gam_terms::construction::canonicalize_penalty_specs(
&penalty_specs,
&design.nullspace_dims,
design.design.ncols(),
"centered tensor penalty regression",
)
.unwrap();
for cp in canonical {
assert_eq!(cp.root.ncols(), cp.col_range.len());
assert_eq!(cp.local.nrows(), cp.col_range.len());
assert_eq!(cp.local.ncols(), cp.col_range.len());
}
}
#[test]
fn periodic_bspline_margin_wraps_exactly_at_period() {
let x = array![0.0, 1.25, 2.5, 3.75, 7.0, 8.25];
let spec = BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::PeriodicUniform {
data_range: (0.0, 1.0),
num_basis: 8,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary_conditions: Default::default(),
boundary: OneDimensionalBoundary::Open,
};
let built = build_bspline_basis_1d(x.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "periodic bspline", e));
let dense = built.design.to_dense();
assert_eq!(dense.ncols(), 8);
for j in 0..dense.ncols() {
assert!(
(dense[[0, j]] - dense[[4, j]]).abs() < 1e-12,
"seam row differs at column {j}: {} vs {}",
dense[[0, j]],
dense[[4, j]]
);
assert!(
(dense[[1, j]] - dense[[5, j]]).abs() < 1e-12,
"wrapped row differs at column {j}: {} vs {}",
dense[[1, j]],
dense[[5, j]]
);
}
for row in dense.rows() {
assert!((row.sum() - 1.0).abs() < 1e-12);
}
assert_eq!(built.nullspace_dims[0], 1);
}
#[test]
fn tensor_bspline_supports_two_periodic_margins_as_torus() {
let data = array![[0.0, 0.0], [7.0, 0.0], [0.0, 24.0], [7.0, 24.0], [1.5, 6.0]];
let spec_day = BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::PeriodicUniform {
data_range: (0.0, 7.0),
num_basis: 7,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary_conditions: Default::default(),
boundary: OneDimensionalBoundary::Open,
};
let spec_hour = BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::PeriodicUniform {
data_range: (0.0, 24.0),
num_basis: 8,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary_conditions: Default::default(),
boundary: OneDimensionalBoundary::Open,
};
let spec_collection = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "te_day_hour".to_string(),
basis: SmoothBasisSpec::TensorBSpline {
feature_cols: vec![0, 1],
spec: TensorBSplineSpec {
marginalspecs: vec![spec_day, spec_hour],
periods: Vec::new(),
double_penalty: false,
identifiability: TensorBSplineIdentifiability::None,
penalty_decomposition: Default::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec_collection)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "periodic tensor design", e));
let sd = &design.smooth;
let dense = sd.term_designs[0].to_dense();
assert_eq!(dense.ncols(), 56);
for j in 0..dense.ncols() {
assert!((dense[[0, j]] - dense[[1, j]]).abs() < 1e-12);
assert!((dense[[0, j]] - dense[[2, j]]).abs() < 1e-12);
assert!((dense[[0, j]] - dense[[3, j]]).abs() < 1e-12);
}
assert_eq!(sd.penalties.len(), 2);
assert!(sd.penalties.iter().all(|p| p.local.nrows() == 56));
let frozen = freeze_term_collection_from_design(&spec_collection, &design)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze periodic tensor", e));
match &frozen.smooth_terms[0].basis {
SmoothBasisSpec::TensorBSpline { spec, .. } => {
assert!(matches!(
spec.marginalspecs[0].knotspec,
BSplineKnotSpec::PeriodicUniform { data_range, .. }
if (data_range.1 - data_range.0 - 7.0).abs() < 1e-9
));
assert!(matches!(
spec.marginalspecs[1].knotspec,
BSplineKnotSpec::PeriodicUniform { data_range, .. }
if (data_range.1 - data_range.0 - 24.0).abs() < 1e-9
));
}
_ => panic!("expected tensor spec"),
}
}
#[test]
fn tensor_bspline_design_matches_extended_marginal_kronecker_product() {
let data = array![[-0.2, 0.1], [0.2, 0.4], [0.7, 0.8], [1.2, 1.1],];
let spec_x = BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 3,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
};
let spec_y = BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 2,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
};
let mx = build_bspline_basis_1d(data.column(0), &spec_x)
.unwrap()
.design
.to_dense();
let my = build_bspline_basis_1d(data.column(1), &spec_y)
.unwrap()
.design
.to_dense();
let expected = tensor_product_design_from_marginals(&[mx.clone(), my.clone()]).unwrap();
let term = SmoothTermSpec {
name: "te_xy".to_string(),
basis: SmoothBasisSpec::TensorBSpline {
feature_cols: vec![0, 1],
spec: TensorBSplineSpec {
periods: Vec::new(),
marginalspecs: vec![spec_x, spec_y],
double_penalty: false,
identifiability: TensorBSplineIdentifiability::None,
penalty_decomposition: Default::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let got = build_smooth_design(data.view(), &[term])
.unwrap()
.term_designs
.into_iter()
.next()
.unwrap()
.to_dense();
assert_eq!(got.dim(), expected.dim());
for i in 0..got.nrows() {
for j in 0..got.ncols() {
assert!((got[[i, j]] - expected[[i, j]]).abs() < 1e-10);
}
}
}
#[test]
fn tensor_bspline_periodic_margins_wrap_as_torus() {
let data = array![[1.25, 3.5], [8.25, 27.5], [-5.75, -20.5], [1.25, 27.5]];
let periodic_day = BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 7.0),
num_internal_knots: 4,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary_conditions: BSplineBoundaryConditions::default(),
boundary: OneDimensionalBoundary::Cyclic {
start: 0.0,
end: 7.0,
},
};
let periodic_hour = BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 24.0),
num_internal_knots: 4,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary_conditions: BSplineBoundaryConditions::default(),
boundary: OneDimensionalBoundary::Cyclic {
start: 0.0,
end: 24.0,
},
};
let term = SmoothTermSpec {
name: "te_day_hour".to_string(),
basis: SmoothBasisSpec::TensorBSpline {
feature_cols: vec![0, 1],
spec: TensorBSplineSpec {
marginalspecs: vec![periodic_day, periodic_hour],
periods: Vec::new(),
double_penalty: false,
identifiability: TensorBSplineIdentifiability::None,
penalty_decomposition: Default::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let design = build_smooth_design(data.view(), &[term])
.unwrap()
.term_designs
.into_iter()
.next()
.unwrap()
.to_dense();
for j in 0..design.ncols() {
assert!((design[[0, j]] - design[[1, j]]).abs() < 1e-12);
assert!((design[[0, j]] - design[[2, j]]).abs() < 1e-12);
}
}
#[test]
fn tensor_bspline_design_is_identifiable_against_global_intercept() {
let n = 120usize;
let mut data = Array2::<f64>::zeros((n, 2));
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
data[[i, 1]] = (3.0 * t).sin();
}
let tensor_term = SmoothTermSpec {
name: "te_xy".to_string(),
basis: SmoothBasisSpec::TensorBSpline {
feature_cols: vec![0, 1],
spec: TensorBSplineSpec {
periods: Vec::new(),
marginalspecs: vec![
BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 6,
},
double_penalty: false,
identifiability: BSplineIdentifiability::default(),
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
},
BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (-1.0, 1.0),
num_internal_knots: 6,
},
double_penalty: false,
identifiability: BSplineIdentifiability::default(),
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
},
],
double_penalty: false,
identifiability: TensorBSplineIdentifiability::default(),
penalty_decomposition: Default::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let sd = build_smooth_design(data.view(), &[tensor_term.clone()]).unwrap();
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![tensor_term],
};
let full = build_term_collection_design(data.view(), &spec).unwrap();
let ones = Array1::<f64>::ones(n);
let sd_dense_terms = sd
.term_designs
.iter()
.map(|d| d.to_dense())
.collect::<Vec<_>>();
let sd_assembled = ndarray::concatenate(
ndarray::Axis(1),
&sd_dense_terms.iter().map(|d| d.view()).collect::<Vec<_>>(),
)
.unwrap();
let residualvs_tensor = residual_norm_to_column_space(&sd_assembled, &ones);
let full_design_dense = full.design.to_dense();
let residualvs_full = residual_norm_to_column_space(&full_design_dense, &ones);
assert!(residualvs_tensor > 1e-6);
assert!(residualvs_full < 1e-8);
}
fn assert_matern_spatial_length_scale_optimization_monotone(
data: ArrayView2<'_, f64>,
y: &Array1<f64>,
weights: &Array1<f64>,
offset: &Array1<f64>,
spec: &TermCollectionSpec,
fit_opts: &FitOptions,
) {
let baseline = fit_term_collection_forspec(
data,
y.view(),
weights.view(),
offset.view(),
spec,
LikelihoodSpec::gaussian_identity(),
fit_opts,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "baseline fit should succeed", e));
let baseline_score = fit_score(&baseline.fit);
let optimized = fit_term_collectionwith_spatial_length_scale_optimization(
data,
y.clone(),
weights.clone(),
offset.clone(),
spec,
LikelihoodSpec::gaussian_identity(),
fit_opts,
&SpatialLengthScaleOptimizationOptions {
max_outer_iter: 16,
rel_tol: 1e-5,
pilot_subsample_threshold: 0,
outer_wall_clock_budget_secs: None,
..SpatialLengthScaleOptimizationOptions::default()
},
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "optimized fit should succeed", e));
let optimized_score = fit_score(&optimized.fit);
assert!(optimized_score <= baseline_score + 1e-10);
let ls = match &optimized.resolvedspec.smooth_terms[0].basis {
SmoothBasisSpec::Matern { spec, .. } => spec.length_scale,
_ => panic!("expected Matérn term"),
};
assert!(ls.is_finite() && (1e-3..=1e3).contains(&ls));
match &optimized.resolvedspec.smooth_terms[0].basis {
SmoothBasisSpec::Matern { spec, .. } => {
assert!(matches!(
spec.center_strategy,
CenterStrategy::UserProvided(_)
));
assert!(matches!(
spec.identifiability,
MaternIdentifiability::FrozenTransform { .. }
));
}
_ => panic!("expected Matérn term"),
}
}
#[test]
fn spatial_length_scale_optimization_monotone_improves_or_keeps_score_for_matern_two_feature() {
let n = 60usize;
let d = 3usize;
let mut data = Array2::<f64>::zeros((n, d));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let x0 = i as f64 / (n as f64 - 1.0);
let x1 = (i as f64 * 0.13).sin();
let x2 = (i as f64 * 0.07).cos();
data[[i, 0]] = x0;
data[[i, 1]] = x1;
data[[i, 2]] = x2;
y[i] = (2.5 * x0).sin() + 0.4 * x1 - 0.2 * x2;
}
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "matern".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1, 2],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 12 },
length_scale: 20.0,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let fit_opts = FitOptions {
max_iter: 40,
..FitOptions::default()
};
let weights = Array1::ones(n);
let offset = Array1::zeros(n);
assert_matern_spatial_length_scale_optimization_monotone(
data.view(),
&y,
&weights,
&offset,
&spec,
&fit_opts,
);
}
pub(super) fn run_two_block_exact_joint_optimize(
data: ArrayView2<'_, f64>,
meanspec: &TermCollectionSpec,
noisespec: &TermCollectionSpec,
expect_msg: &str,
) -> SpatialLengthScaleOptimizationResult<f64> {
let kappa_options = SpatialLengthScaleOptimizationOptions {
max_outer_iter: 1,
rel_tol: 1e-6,
pilot_subsample_threshold: 0,
outer_wall_clock_budget_secs: None,
..SpatialLengthScaleOptimizationOptions::default()
};
let joint_setup = two_block_exact_joint_hyper_setup(meanspec, noisespec, &kappa_options);
let theta_dim = joint_setup.theta0().len();
let mean_terms = spatial_length_scale_term_indices(meanspec);
let noise_terms = spatial_length_scale_term_indices(noisespec);
let policy = gam_model_api::families::custom_family::OuterDerivativePolicy {
capability: gam_problem::ExactOuterDerivativeOrder::Second,
predicted_hessian_work: 0,
predicted_gradient_work: 0,
subsample_capable: false,
};
optimize_spatial_length_scale_exact_joint(
data,
&[meanspec.clone(), noisespec.clone()],
&[mean_terms, noise_terms],
&kappa_options,
&joint_setup,
gam_problem::SeedRiskProfile::Gaussian,
true,
true,
false,
None,
policy,
|theta, specs, designs| {
assert_eq!(theta.len(), theta_dim);
assert_eq!(specs.len(), 2);
Ok(designs[0].design.ncols() as f64
+ designs[1].design.ncols() as f64
+ designs[0].penalties.len() as f64
+ designs[1].penalties.len() as f64)
},
|theta, specs, designs, eval_mode, _| {
assert_eq!(theta.len(), theta_dim);
assert_eq!(specs.len(), 2);
assert!(!designs.is_empty());
Ok((
0.0,
Array1::zeros(theta_dim),
if matches!(
eval_mode,
gam_solve::estimate::reml::reml_outer_engine::EvalMode::ValueGradientHessian
) {
gam_problem::HessianResult::Analytic(Array2::zeros((
theta_dim, theta_dim,
)))
} else {
gam_problem::HessianResult::Unavailable
},
))
},
|theta, specs, designs| {
assert_eq!(theta.len(), theta_dim);
assert_eq!(specs.len(), 2);
assert!(!designs.is_empty());
Ok(gam_problem::EfsEval {
cost: 0.0,
steps: vec![0.0; theta_dim],
beta: None,
psi_gradient: None,
psi_indices: None,
inner_hessian_scale: None,
logdet_enclosure_gap: None,
})
},
|_beta: &Array1<f64>| Ok(gam_solve::rho_optimizer::SeedOutcome::NoSlot),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", expect_msg, e))
}
#[test]
fn exact_joint_two_block_spatial_length_scale_freezes_matern_centers() {
let n = 40usize;
let mut data = Array2::<f64>::zeros((n, 2));
for i in 0..n {
let x0 = i as f64 / (n as f64 - 1.0);
let x1 = (i as f64 * 0.21).sin();
data[[i, 0]] = x0;
data[[i, 1]] = x1;
}
let matern_term = |name: &str, length_scale: f64| SmoothTermSpec {
name: name.to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
length_scale,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: Some(vec![0.0, 0.0]),
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let meanspec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![matern_term("mean_matern", 0.8)],
};
let noisespec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![matern_term("noise_matern", 1.1)],
};
let solved = run_two_block_exact_joint_optimize(
data.view(),
&meanspec,
&noisespec,
"exact joint two-block κ optimization should succeed",
);
for resolved in [&solved.resolved_specs[0], &solved.resolved_specs[1]] {
match &resolved.smooth_terms[0].basis {
SmoothBasisSpec::Matern { spec, .. } => {
assert!(matches!(
spec.center_strategy,
CenterStrategy::UserProvided(_)
));
assert!(matches!(
spec.identifiability,
MaternIdentifiability::FrozenTransform { .. }
));
}
_ => panic!("expected Matérn term"),
}
}
}
#[test]
fn exact_joint_spatial_outer_hessian_available_for_dense_non_gaussian_designs() {
let n = 5_001usize;
let p = 100usize;
let data = Array2::from_shape_fn((n, p), |(i, j)| ((i + j + 1) as f64).sin());
let spec = TermCollectionSpec {
linear_terms: (0..p)
.map(|j| LinearTermSpec {
name: format!("x{j}"),
feature_col: j,
feature_cols: vec![j],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
})
.collect(),
random_effect_terms: vec![],
smooth_terms: vec![],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
assert!(exact_joint_spatial_outer_hessian_available(
&LikelihoodSpec::binomial_logit(),
&design,
));
assert!(exact_joint_spatial_outer_hessian_available(
&LikelihoodSpec::gaussian_identity(),
&design,
));
}
#[test]
fn spatial_aniso_joint_exact_hessian_materializes_small_case() {
let n = 18usize;
let mut data = Array2::<f64>::zeros((n, 2));
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
data[[i, 1]] = (0.41 * i as f64).sin();
}
let y = Array1::from_iter((0..n).map(|i| {
let t = i as f64 / (n as f64 - 1.0);
0.4 + (2.0 * std::f64::consts::PI * t).sin()
}));
let weights = Array1::ones(n);
let offset = Array1::zeros(n);
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "matern_aniso".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 5 },
length_scale: 0.85,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: Some(vec![0.2, -0.2]),
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let fit_opts = FitOptions {
compute_inference: false,
max_iter: 120,
tol: 1e-10,
penalty_shrinkage_floor: None,
..FitOptions::default()
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let frozen = freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze", e));
let spatial_terms = spatial_length_scale_term_indices(&frozen);
let dims_per_term = spatial_dims_per_term(&frozen, &spatial_terms);
assert_eq!(dims_per_term, vec![2]);
let rho_dim = design.penalties.len();
let log_kappa0 = SpatialLogKappaCoords::from_length_scales_aniso(
&frozen,
&spatial_terms,
&SpatialLengthScaleOptimizationOptions::default(),
);
let mut theta = Array1::<f64>::zeros(rho_dim + log_kappa0.as_array().len());
for j in 0..rho_dim {
theta[j] = -0.15 + 0.07 * j as f64;
}
theta.slice_mut(s![rho_dim..]).assign(log_kappa0.as_array());
let external_opts =
external_opts_for_design(&LikelihoodSpec::gaussian_identity(), &design, &fit_opts);
let mut cache = SingleBlockExactJointDesignCache::new(
data.view(),
frozen,
design.clone(),
spatial_terms,
rho_dim,
dims_per_term,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "single-block cache", e));
let mut evaluator = gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&design.design,
offset.view(),
&design.penalties,
&external_opts,
"small aniso Hessian finite-difference evaluator",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "evaluator", e));
let eval_at = |theta: &Array1<f64>,
cache: &mut SingleBlockExactJointDesignCache<'_>,
evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>,
order: gam_model_api::OuterEvalOrder| {
cache.ensure_theta(theta).unwrap_or_else(|e| panic!("{} failed: {:?}", "theta applied", e));
let hyper_dirs = try_build_spatial_log_kappa_hyper_dirs(
data.view(),
cache.spec(),
cache.design(),
&cache.spatial_terms,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "hyper dirs build", e))
.expect("hyper dirs present");
evaluate_joint_reml_outer_eval_at_theta(
evaluator,
cache.design(),
theta,
rho_dim,
hyper_dirs,
None,
order,
None,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "outer eval", e))
};
let (_, gradient, hessian_result) = eval_at(
&theta,
&mut cache,
&mut evaluator,
gam_model_api::OuterEvalOrder::ValueGradientHessian,
);
let hessian = hessian_result
.materialize_dense()
.unwrap_or_else(|e| panic!("{} failed: {:?}", "hessian materializes", e))
.expect("hessian present");
assert_eq!(hessian.nrows(), theta.len());
assert_eq!(hessian.ncols(), theta.len());
assert!(hessian.iter().all(|value| value.is_finite()));
assert!(gradient.iter().all(|value| value.is_finite()));
let symmetry_diff = max_abs_diff_matrix(&hessian, &hessian.t().to_owned());
assert!(
symmetry_diff <= 1e-10,
"small aniso exact Hessian should be symmetric, max diff={symmetry_diff}"
);
let psi_block = hessian.slice(s![rho_dim.., rho_dim..]).to_owned();
assert!(
psi_block.iter().any(|value| value.abs() > 1e-10),
"small aniso exact Hessian should carry non-zero ψ curvature"
);
}
#[test]
fn iso_kappa_duchon_binomial_probit_joint_gradient_matches_finite_difference() {
let (pass, worst, violations) = iso_kappa_fd_variant_driver(
"duchon_probit_n80",
80,
LikelihoodSpec::binomial_probit(),
false,
false,
);
assert!(
pass,
"Duchon BinomialProbit n=80 FD failed; worst_psi_rel={worst:.3e}\n {}",
violations.join("\n ")
);
}
fn iso_kappa_fd_variant_driver(
label: &str,
n: usize,
family: LikelihoodSpec,
skip_psi: bool,
well_conditioned: bool,
) -> (bool, f64, Vec<String>) {
let mut data = Array2::<f64>::zeros((n, 1));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
let eta = 1.4 * (2.0 * std::f64::consts::PI * t).sin() + 0.5 * (t - 0.5);
let raw = eta + 0.7 * (3.7 * (i as f64) + 1.0).sin();
y[i] = if family.is_gaussian_identity() {
raw
} else if well_conditioned {
let p = 1.0 / (1.0 + (-0.6 * (2.0 * std::f64::consts::PI * t).sin()).exp());
let u = 0.5 * ((5.0 * (i as f64) + 0.5).sin() + 1.0);
if u < p { 1.0 } else { 0.0 }
} else if raw > 0.0 {
1.0
} else {
0.0
};
}
let weights = Array1::ones(n);
let offset = Array1::zeros(n);
let basis = if label.starts_with("matern") {
SmoothBasisSpec::Matern {
feature_cols: vec![0],
spec: MaternBasisSpec {
center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
periodic: None,
length_scale: 1.0,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: false,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
}
} else {
SmoothBasisSpec::Duchon {
feature_cols: vec![0],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
length_scale: Some(1.0),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
}
};
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "variant_1d".to_string(),
basis,
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let fit_opts = FitOptions {
compute_inference: false,
max_iter: 200,
tol: 1e-12,
penalty_shrinkage_floor: None,
..FitOptions::default()
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let frozen = freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze", e));
let frozen_design = build_term_collection_design(data.view(), &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "frozen design", e));
let spatial_terms = spatial_length_scale_term_indices(&frozen);
let dims_per_term = spatial_dims_per_term(&frozen, &spatial_terms);
assert_eq!(dims_per_term, vec![1], "{label}: expect one log-κ axis");
let rho_dim = frozen_design.penalties.len();
let psi_dim: usize = dims_per_term.iter().sum();
assert!(psi_dim >= 1);
let external_opts = external_opts_for_design(&family, &frozen_design, &fit_opts);
let mut cache = SingleBlockExactJointDesignCache::new(
data.view(),
frozen.clone(),
frozen_design.clone(),
spatial_terms.clone(),
rho_dim,
dims_per_term.clone(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "single-block cache", e));
let mut evaluator = gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&frozen_design.design,
offset.view(),
&frozen_design.penalties,
&external_opts,
"iso-κ variant FD evaluator",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "evaluator", e));
let cost_at = |theta: &Array1<f64>,
cache: &mut SingleBlockExactJointDesignCache<'_>,
evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>|
-> f64 {
cache.ensure_theta(theta).unwrap_or_else(|e| panic!("{} failed: {:?}", "ensure_theta", e));
let design = cache.design();
evaluator
.evaluate_cost_only(
&design.design,
&design.penalties,
&design.nullspace_dims,
design.linear_constraints.clone(),
theta,
rho_dim,
None,
"iso-κ variant FD cost-only",
None,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "cost-only eval", e))
};
let analytic_at = |theta: &Array1<f64>,
cache: &mut SingleBlockExactJointDesignCache<'_>,
evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>|
-> (f64, Array1<f64>) {
cache.ensure_theta(theta).expect("ensure_theta");
let hyper_dirs = try_build_spatial_log_kappa_hyper_dirs(
data.view(),
cache.spec(),
cache.design(),
&cache.spatial_terms,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "hyper dirs build", e))
.expect("hyper dirs present");
let (cost, grad, _hess) = evaluate_joint_reml_outer_eval_at_theta(
evaluator,
cache.design(),
theta,
rho_dim,
hyper_dirs,
None,
gam_model_api::OuterEvalOrder::ValueAndGradient,
None,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "outer eval", e));
(cost, grad)
};
let theta_dim = rho_dim + psi_dim;
let theta_zero = Array1::<f64>::zeros(theta_dim);
let mut theta_base = Array1::<f64>::zeros(theta_dim);
for j in 0..rho_dim {
theta_base[j] = 0.2 - 0.1 * j as f64;
}
let mut theta_psi_only = Array1::<f64>::zeros(theta_dim);
for k in 0..psi_dim {
theta_psi_only[rho_dim + k] = 0.4;
}
let mut theta_alt = theta_base.clone();
for j in 0..rho_dim {
theta_alt[j] = 1.0 + 0.05 * j as f64;
}
for k in 0..psi_dim {
theta_alt[rho_dim + k] = 0.4;
}
let h = 1e-5_f64;
let rel_tol = 5e-3_f64;
let mut violations: Vec<String> = Vec::new();
let mut worst_psi_rel = 0.0_f64;
for (probe, theta) in [
("zero", &theta_zero),
("psi_only", &theta_psi_only),
("base", &theta_base),
("alt", &theta_alt),
] {
let (cost_an, grad_an) = analytic_at(theta, &mut cache, &mut evaluator);
assert!(cost_an.is_finite(), "{label} {probe}: cost not finite");
let cost_via_fd_path = cost_at(theta, &mut cache, &mut evaluator);
eprintln!(
"[{label} {probe}] COST an={:+.10e} fd_path={:+.10e} diff={:.3e}",
cost_an,
cost_via_fd_path,
(cost_an - cost_via_fd_path).abs()
);
for j in 0..theta_dim {
let is_psi = j >= rho_dim;
if skip_psi && is_psi {
continue;
}
let mut plus = theta.clone();
plus[j] += h;
let mut minus = theta.clone();
minus[j] -= h;
let cp = cost_at(&plus, &mut cache, &mut evaluator);
let cm = cost_at(&minus, &mut cache, &mut evaluator);
let fd = (cp - cm) / (2.0 * h);
let denom = fd.abs().max(grad_an[j].abs()).max(1e-3);
let rel = (grad_an[j] - fd).abs() / denom;
let kind = if is_psi { "psi" } else { "rho" };
eprintln!(
"[{label} {probe}] {kind} j={j} an={:+.4e} fd={:+.4e} rel={:.3e}",
grad_an[j], fd, rel
);
if is_psi && rel > worst_psi_rel {
worst_psi_rel = rel;
}
if rel >= rel_tol {
violations.push(format!(
"{probe} {kind} j={j}: analytic={:+.6e} fd={:+.6e} rel={:.3e}",
grad_an[j], fd, rel
));
}
}
}
let pass = violations.is_empty();
eprintln!(
"[{label} SUMMARY] pass={pass} worst_psi_rel={worst_psi_rel:.3e} \
violations={}",
violations.len()
);
(pass, worst_psi_rel, violations)
}
#[test]
fn iso_kappa_duchon_gaussian_identity_fd() {
let (pass, worst, violations) = iso_kappa_fd_variant_driver(
"duchon_gaussian",
80,
LikelihoodSpec::gaussian_identity(),
false,
false,
);
assert!(
pass,
"Gaussian Identity FD failed; worst_psi_rel={worst:.3e}\n {}",
violations.join("\n ")
);
}
#[test]
fn iso_kappa_matern_gaussian_identity_fd() {
let (pass, worst, violations) = iso_kappa_fd_variant_driver(
"matern_gaussian",
80,
LikelihoodSpec::gaussian_identity(),
false,
false,
);
assert!(
pass,
"Matérn iso-κ Gaussian-identity outer-gradient FD failed; \
worst_psi_rel={worst:.3e}\n {}",
violations.join("\n ")
);
}
#[test]
fn exact_spatial_joint_engine_aniso_iso_parity_1d() {
let n = 80usize;
let mut data = Array2::<f64>::zeros((n, 1));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
let eta = 1.4 * (2.0 * std::f64::consts::PI * t).sin() + 0.5 * (t - 0.5);
y[i] = eta + 0.7 * (3.7 * (i as f64) + 1.0).sin();
}
let weights = Array1::<f64>::ones(n);
let offset = Array1::<f64>::zeros(n);
let family = LikelihoodSpec::gaussian_identity();
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "parity_1d".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
length_scale: Some(1.0),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::all_active(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let fit_opts = FitOptions {
compute_inference: false,
max_iter: 200,
tol: 1e-12,
penalty_shrinkage_floor: None,
..FitOptions::default()
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let frozen = freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze", e));
let frozen_design = build_term_collection_design(data.view(), &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "frozen design", e));
let spatial_terms = spatial_length_scale_term_indices(&frozen);
assert_eq!(spatial_terms.len(), 1, "expect a single spatial term");
let dims_per_term = spatial_dims_per_term(&frozen, &spatial_terms);
assert_eq!(dims_per_term, vec![1], "expect one log-scale axis");
let rho_dim = frozen_design.penalties.len();
assert!(rho_dim >= 1, "expect at least one penalty block");
const JOINT_RHO_BOUND: f64 = 12.0;
let kappa_options = SpatialLengthScaleOptimizationOptions::default();
let log_kappa0 =
SpatialLogKappaCoords::from_length_scales(&frozen, &spatial_terms, &kappa_options);
let log_kappa_lower = SpatialLogKappaCoords::lower_bounds_from_data(
data.view(),
&frozen,
&spatial_terms,
&kappa_options,
);
let log_kappa_upper = SpatialLogKappaCoords::upper_bounds_from_data(
data.view(),
&frozen,
&spatial_terms,
&kappa_options,
);
let log_kappa0 = log_kappa0.clamp_to_bounds(&log_kappa_lower, &log_kappa_upper);
let setup = ExactJointHyperSetup::new(
Array1::<f64>::zeros(rho_dim), Array1::<f64>::from_elem(rho_dim, -JOINT_RHO_BOUND),
Array1::<f64>::from_elem(rho_dim, JOINT_RHO_BOUND),
log_kappa0,
log_kappa_lower,
log_kappa_upper,
);
let theta0 = setup.theta0();
let lower = setup.lower();
let upper = setup.upper();
let run = |kind: SpatialHyperKind| -> (Array1<f64>, f64) {
run_exact_joint_spatial_optimization(
kind,
data.view(),
y.view(),
weights.view(),
offset.view(),
&frozen,
&frozen_design,
family.clone(),
&fit_opts,
&spatial_terms,
&dims_per_term,
&theta0,
&lower,
&upper,
rho_dim,
&kappa_options,
)
.map(|(outcome, _timing)| match outcome {
SpatialJointOutcome::Optimized {
theta_star,
final_value,
} => (theta_star, final_value),
SpatialJointOutcome::NonConverged { final_value, .. } => panic!(
"exact joint spatial optimization did not converge (final_value={final_value})"
),
})
.unwrap_or_else(|e| panic!("{} failed: {:?}", "exact joint spatial optimization", e))
};
let (theta_aniso, value_aniso) = run(SpatialHyperKind::Anisotropic);
let (theta_iso, value_iso) = run(SpatialHyperKind::Isotropic);
assert_eq!(
theta_aniso.len(),
theta_iso.len(),
"converged θ dimension must match across coordinate kinds"
);
for j in 0..theta_aniso.len() {
let diff = (theta_aniso[j] - theta_iso[j]).abs();
assert!(
diff <= 1e-9 * (1.0 + theta_aniso[j].abs()),
"θ[{j}] differs across kinds: aniso={:+.12e} iso={:+.12e} diff={:.3e}",
theta_aniso[j],
theta_iso[j],
diff,
);
}
let value_diff = (value_aniso - value_iso).abs();
assert!(
value_diff <= 1e-9 * (1.0 + value_aniso.abs()),
"final REML value differs across kinds: aniso={:+.12e} iso={:+.12e} diff={:.3e}",
value_aniso,
value_iso,
value_diff,
);
assert!(
value_aniso.is_finite() && value_iso.is_finite(),
"both kinds must produce a finite certified REML cost"
);
}
#[test]
fn psi_gram_tensor_lane_matches_streamed_reml_cost_and_gradient() {
use gam_model_api::OuterEvalOrder;
let n = 600usize;
let mut data = Array2::<f64>::zeros((n, 1));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
let signal = 1.2 * (2.0 * std::f64::consts::PI * t).sin() + 0.4 * (t - 0.5);
let noise = 0.15 * (((i as f64) * 12.9898).sin() * 43758.547).fract();
y[i] = signal + noise;
}
let weights = Array1::<f64>::ones(n);
let offset = Array1::<f64>::zeros(n);
let family = LikelihoodSpec::gaussian_identity();
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "psi_tensor_invariance".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 12 },
length_scale: Some(1.0),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::all_active(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let fit_opts = FitOptions {
compute_inference: false,
max_iter: 200,
tol: 1e-12,
penalty_shrinkage_floor: None,
..FitOptions::default()
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let frozen = freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze", e));
let frozen_design = build_term_collection_design(data.view(), &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "frozen design", e));
let spatial_terms = spatial_length_scale_term_indices(&frozen);
assert_eq!(spatial_terms.len(), 1, "expect a single spatial term");
let dims_per_term = spatial_dims_per_term(&frozen, &spatial_terms);
assert_eq!(
dims_per_term,
vec![1],
"expect one log-scale axis (coord_dim == 1)"
);
let rho_dim = frozen_design.penalties.len();
assert!(rho_dim >= 1, "expect at least one penalty block");
let kappa_options = SpatialLengthScaleOptimizationOptions::default();
let log_kappa0 =
SpatialLogKappaCoords::from_length_scales(&frozen, &spatial_terms, &kappa_options);
let log_kappa_lower = SpatialLogKappaCoords::lower_bounds_from_data(
data.view(),
&frozen,
&spatial_terms,
&kappa_options,
);
let log_kappa_upper = SpatialLogKappaCoords::upper_bounds_from_data(
data.view(),
&frozen,
&spatial_terms,
&kappa_options,
);
let log_kappa0 = log_kappa0.clamp_to_bounds(&log_kappa_lower, &log_kappa_upper);
const JOINT_RHO_BOUND: f64 = 12.0;
let setup = ExactJointHyperSetup::new(
Array1::<f64>::zeros(rho_dim),
Array1::<f64>::from_elem(rho_dim, -JOINT_RHO_BOUND),
Array1::<f64>::from_elem(rho_dim, JOINT_RHO_BOUND),
log_kappa0,
log_kappa_lower,
log_kappa_upper,
);
let theta0 = setup.theta0();
let lower = setup.lower();
let upper = setup.upper();
let psi_lo = lower[rho_dim];
let psi_hi = upper[rho_dim];
assert!(psi_hi > psi_lo, "ψ window must be non-degenerate");
let make_cache = || {
SingleBlockExactJointDesignCache::new(
data.view(),
frozen.clone(),
frozen_design.clone(),
spatial_terms.clone(),
rho_dim,
dims_per_term.clone(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "design cache", e))
};
let external_opts = external_opts_for_design(&family, &frozen_design, &fit_opts);
let mut streamed_eval = gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&frozen_design.design,
offset.view(),
&frozen_design.penalties,
&external_opts,
"psi_tensor_invariance/streamed",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "streamed evaluator", e));
let mut tensor_eval = gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&frozen_design.design,
offset.view(),
&frozen_design.penalties,
&external_opts,
"psi_tensor_invariance/tensor",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "tensor evaluator", e));
let z = Array1::from_iter(y.iter().zip(offset.iter()).map(|(yi, oi)| yi - oi));
let attached = {
let mut build_cache = make_cache();
let theta_probe_base = theta0.clone();
tensor_eval.build_and_set_psi_gram_tensor(
|psi| {
let mut theta_probe = theta_probe_base.clone();
theta_probe[rho_dim] = psi;
build_cache.ensure_theta(&theta_probe)?;
Ok(build_cache.design().design.clone())
},
weights.view(),
z.view(),
psi_lo,
psi_hi,
)
};
assert!(
attached,
"ψ-gram tensor failed to certify over the production window \
[{psi_lo:.3}, {psi_hi:.3}]; the invariance test would be vacuous"
);
let mut stream_cache = make_cache();
let mut tensor_cache = make_cache();
let psi_samples = [
psi_lo + 0.10 * (psi_hi - psi_lo),
psi_lo + 0.37 * (psi_hi - psi_lo),
0.5 * (psi_lo + psi_hi),
psi_lo + 0.78 * (psi_hi - psi_lo),
psi_hi - 0.05 * (psi_hi - psi_lo),
];
let rho_samples = [
Array1::<f64>::from_elem(rho_dim, -1.5),
Array1::<f64>::from_elem(rho_dim, 0.5),
];
let eval_one = |evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>,
cache: &mut SingleBlockExactJointDesignCache<'_>,
theta: &Array1<f64>,
with_hessian: bool|
-> (f64, Array1<f64>, Option<Array2<f64>>) {
use gam_problem::HessianResult;
cache.ensure_theta(theta).unwrap_or_else(|e| panic!("{} failed: {:?}", "ensure_theta", e));
let hyper_dirs = try_build_spatial_log_kappa_hyper_dirs(
data.view(),
cache.spec(),
cache.design(),
&spatial_terms,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "hyper_dirs build", e))
.expect("hyper_dirs present");
let design_revision = Some(cache.design_revision());
let order = if with_hessian {
OuterEvalOrder::ValueGradientHessian
} else {
OuterEvalOrder::ValueAndGradient
};
let (cost, grad, hess) = evaluate_joint_reml_outer_eval_at_theta(
evaluator,
cache.design(),
theta,
rho_dim,
hyper_dirs,
None,
order,
design_revision,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "evaluate_with_order", e));
let hess_mat = if with_hessian {
match hess {
HessianResult::Analytic(h) => Some(h),
_ => None,
}
} else {
None
};
(cost, grad, hess_mat)
};
let mut worst_cost_rel = 0.0_f64;
let mut worst_grad_abs = 0.0_f64;
let mut worst_hess_abs = 0.0_f64;
for with_hessian in [false, true] {
for rho in &rho_samples {
for &psi in &psi_samples {
assert!(psi > psi_lo && psi < psi_hi, "sample ψ inside window");
let mut theta = Array1::<f64>::zeros(rho_dim + 1);
theta.slice_mut(s![..rho_dim]).assign(rho);
theta[rho_dim] = psi;
let (cost_s, grad_s, hess_s) =
eval_one(&mut streamed_eval, &mut stream_cache, &theta, with_hessian);
let (cost_t, grad_t, hess_t) =
eval_one(&mut tensor_eval, &mut tensor_cache, &theta, with_hessian);
assert!(
cost_s.is_finite() && cost_t.is_finite(),
"non-finite REML cost at ψ={psi:.4} hessian={with_hessian}: \
streamed={cost_s}, tensor={cost_t}"
);
let cost_rel = (cost_s - cost_t).abs() / (1.0 + cost_s.abs());
worst_cost_rel = worst_cost_rel.max(cost_rel);
assert!(
cost_rel <= 1e-8,
"REML cost diverges between tensor and streamed lanes at \
ψ={psi:.4}, ρ={:+.2} hessian={with_hessian}: \
streamed={cost_s:.12e}, tensor={cost_t:.12e}, rel={cost_rel:.3e}",
rho[0],
);
assert_eq!(grad_s.len(), grad_t.len(), "gradient dimension mismatch");
for j in 0..grad_s.len() {
let gabs = (grad_s[j] - grad_t[j]).abs();
let grel = gabs / (1.0 + grad_s[j].abs());
worst_grad_abs = worst_grad_abs.max(gabs);
assert!(
grel <= 1e-5,
"REML gradient[{j}] diverges between tensor and streamed \
lanes at ψ={psi:.4}, ρ={:+.2} hessian={with_hessian}: \
streamed={:+.12e}, tensor={:+.12e}, |Δ|={gabs:.3e}, \
rel={grel:.3e} (far above summation-order roundoff ⇒ \
∂(XᵀWX)/∂ψ install has a frame/scaling bug)",
rho[0],
grad_s[j],
grad_t[j],
);
}
if let (Some(hs), Some(ht)) = (hess_s, hess_t) {
assert_eq!(
hs.shape(),
ht.shape(),
"Hessian shape mismatch at ψ={psi:.4} ρ={:+.2}",
rho[0],
);
for (((r, c), vs), (_, vt)) in hs.indexed_iter().zip(ht.indexed_iter()) {
let habs = (vs - vt).abs();
let hrel = habs / (1.0 + vs.abs());
worst_hess_abs = worst_hess_abs.max(habs);
assert!(
hrel <= 1e-6,
"REML Hessian[{r},{c}] diverges between tensor and \
streamed lanes at ψ={psi:.4}, ρ={:+.2}: \
streamed={vs:+.12e}, tensor={vt:+.12e}, \
|Δ|={habs:.3e}, rel={hrel:.3e} (both lanes use \
the slab for τ-τ Hessian — divergence is a bug)",
rho[0],
);
}
}
}
}
}
eprintln!(
"[psi-gram-tensor invariance] worst cost rel={worst_cost_rel:.3e}, \
worst grad |Δ|={worst_grad_abs:.3e}, worst hess |Δ|={worst_hess_abs:.3e} \
over {} (ρ,ψ) points × 2 orders",
rho_samples.len() * psi_samples.len(),
);
}
#[test]
fn psi_gram_tensor_e2e_kappa_optimum_matches_streamed() {
let n = 600usize;
let mut data = Array2::<f64>::zeros((n, 1));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
let signal = 1.2 * (2.0 * std::f64::consts::PI * t).sin() + 0.4 * (t - 0.5);
let noise = 0.15 * (((i as f64) * 12.9898).sin() * 43758.547).fract();
y[i] = signal + noise;
}
let weights = Array1::<f64>::ones(n);
let offset = Array1::<f64>::zeros(n);
let family = LikelihoodSpec::gaussian_identity();
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "e2e_kappa_optimum".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 12 },
length_scale: Some(1.0),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::all_active(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let frozen = freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze", e));
let frozen_design = build_term_collection_design(data.view(), &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "frozen design", e));
let spatial_terms = spatial_length_scale_term_indices(&frozen);
let dims_per_term = spatial_dims_per_term(&frozen, &spatial_terms);
let rho_dim = frozen_design.penalties.len();
let kappa_options = SpatialLengthScaleOptimizationOptions::default();
let log_kappa0 =
SpatialLogKappaCoords::from_length_scales(&frozen, &spatial_terms, &kappa_options);
let log_kappa_lower = SpatialLogKappaCoords::lower_bounds_from_data(
data.view(),
&frozen,
&spatial_terms,
&kappa_options,
);
let log_kappa_upper = SpatialLogKappaCoords::upper_bounds_from_data(
data.view(),
&frozen,
&spatial_terms,
&kappa_options,
);
let log_kappa0 = log_kappa0.clamp_to_bounds(&log_kappa_lower, &log_kappa_upper);
const JOINT_RHO_BOUND: f64 = 12.0;
let setup = ExactJointHyperSetup::new(
Array1::<f64>::zeros(rho_dim),
Array1::<f64>::from_elem(rho_dim, -JOINT_RHO_BOUND),
Array1::<f64>::from_elem(rho_dim, JOINT_RHO_BOUND),
log_kappa0.clone(),
log_kappa_lower.clone(),
log_kappa_upper.clone(),
);
let theta0 = setup.theta0();
let lower = setup.lower();
let upper = setup.upper();
let psi_lo = lower[rho_dim];
let psi_hi = upper[rho_dim];
let z = Array1::from_iter(y.iter().zip(offset.iter()).map(|(yi, oi)| yi - oi));
let external_opts = external_opts_for_design(
&family,
&frozen_design,
&FitOptions {
compute_inference: false,
max_iter: 200,
tol: 1e-12,
penalty_shrinkage_floor: None,
..FitOptions::default()
},
);
let make_eval = || {
gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&frozen_design.design,
offset.view(),
&frozen_design.penalties,
&external_opts,
"e2e_kappa_optimum",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "evaluator", e))
};
let make_cache = || {
SingleBlockExactJointDesignCache::new(
data.view(),
frozen.clone(),
frozen_design.clone(),
spatial_terms.clone(),
rho_dim,
dims_per_term.clone(),
)
.expect("design cache")
};
let mut streamed_eval = make_eval();
let mut stream_cache = make_cache();
let mut tensor_eval = make_eval();
let mut tensor_cache = make_cache();
let attached = {
let mut build_cache = make_cache();
let theta_probe_base = theta0.clone();
tensor_eval.build_and_set_psi_gram_tensor(
|psi| {
let mut theta_probe = theta_probe_base.clone();
theta_probe[rho_dim] = psi;
build_cache.ensure_theta(&theta_probe)?;
Ok(build_cache.design().design.clone())
},
weights.view(),
z.view(),
psi_lo,
psi_hi,
)
};
assert!(
attached,
"tensor must certify on this fixture for a non-vacuous gate"
);
let check_theta = theta0.clone();
stream_cache.ensure_theta(&check_theta).unwrap();
tensor_cache.ensure_theta(&check_theta).unwrap();
let hyper_s = try_build_spatial_log_kappa_hyper_dirs(
data.view(),
stream_cache.spec(),
stream_cache.design(),
&spatial_terms,
)
.unwrap()
.unwrap();
let hyper_t = try_build_spatial_log_kappa_hyper_dirs(
data.view(),
tensor_cache.spec(),
tensor_cache.design(),
&spatial_terms,
)
.unwrap()
.unwrap();
let (c_s, g_s, _) = evaluate_joint_reml_outer_eval_at_theta(
&mut streamed_eval,
stream_cache.design(),
&check_theta,
rho_dim,
hyper_s,
None,
OuterEvalOrder::ValueAndGradient,
Some(stream_cache.design_revision()),
)
.unwrap();
let (c_t, g_t, _) = evaluate_joint_reml_outer_eval_at_theta(
&mut tensor_eval,
tensor_cache.design(),
&check_theta,
rho_dim,
hyper_t,
None,
OuterEvalOrder::ValueAndGradient,
Some(tensor_cache.design_revision()),
)
.unwrap();
let cost_rel = (c_s - c_t).abs() / (1.0 + c_s.abs());
assert!(
cost_rel <= 1e-8,
"e2e smoke-check: cost diverges at θ₀: streamed={c_s:.10e} tensor={c_t:.10e} rel={cost_rel:.3e}"
);
for j in 0..g_s.len() {
let grel = (g_s[j] - g_t[j]).abs() / (1.0 + g_s[j].abs());
assert!(
grel <= 1e-5,
"e2e smoke-check: gradient[{j}] diverges at θ₀: \
streamed={:+.10e} tensor={:+.10e} rel={grel:.3e}",
g_s[j],
g_t[j],
);
}
eprintln!(
"[psi-gram-tensor e2e] θ₀ smoke-check: cost rel={cost_rel:.3e}, \
max grad rel={:.3e} — tensor lane bit-tight at the optimizer entry point",
g_s.iter()
.zip(g_t.iter())
.map(|(a, b)| (a - b).abs() / (1.0 + a.abs()))
.fold(0.0_f64, f64::max),
);
let psi_sweep = [
psi_lo + 0.12 * (psi_hi - psi_lo),
psi_lo + 0.40 * (psi_hi - psi_lo),
0.5 * (psi_lo + psi_hi),
psi_lo + 0.71 * (psi_hi - psi_lo),
psi_hi - 0.08 * (psi_hi - psi_lo),
];
let rho_sweep = [
Array1::<f64>::from_elem(rho_dim, -2.0),
Array1::<f64>::from_elem(rho_dim, 0.0),
Array1::<f64>::from_elem(rho_dim, 1.5),
];
let mut worst_beta_abs = 0.0_f64;
let beta_one = |evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>,
cache: &mut SingleBlockExactJointDesignCache<'_>,
theta: &Array1<f64>|
-> Array1<f64> {
cache.ensure_theta(theta).unwrap_or_else(|e| panic!("{} failed: {:?}", "ensure_theta", e));
let hyper = try_build_spatial_log_kappa_hyper_dirs(
data.view(),
cache.spec(),
cache.design(),
&spatial_terms,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "hyper_dirs build", e))
.expect("hyper_dirs present");
evaluate_joint_reml_outer_eval_at_theta(
evaluator,
cache.design(),
theta,
rho_dim,
hyper,
None,
OuterEvalOrder::ValueAndGradient,
Some(cache.design_revision()),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "evaluate_with_order", e));
evaluator
.current_beta()
.unwrap_or_else(|| panic!("{} failed", "converged inner β̂ available after the PIRLS solve"))
};
for rho in &rho_sweep {
for &psi in &psi_sweep {
assert!(psi > psi_lo && psi < psi_hi, "sweep ψ inside window");
let mut theta = Array1::<f64>::zeros(rho_dim + 1);
theta.slice_mut(s![..rho_dim]).assign(rho);
theta[rho_dim] = psi;
let beta_s = beta_one(&mut streamed_eval, &mut stream_cache, &theta);
let beta_t = beta_one(&mut tensor_eval, &mut tensor_cache, &theta);
assert_eq!(
beta_s.len(),
beta_t.len(),
"coefficient dimension mismatch at ψ={psi:.4}"
);
for j in 0..beta_s.len() {
assert!(
beta_s[j].is_finite() && beta_t[j].is_finite(),
"non-finite β̂[{j}] at ψ={psi:.4}: streamed={}, tensor={}",
beta_s[j],
beta_t[j],
);
let babs = (beta_s[j] - beta_t[j]).abs();
let brel = babs / (1.0 + beta_s[j].abs());
worst_beta_abs = worst_beta_abs.max(babs);
assert!(
brel <= 1e-6,
"converged β̂[{j}] diverges between tensor and streamed lanes \
at ψ={psi:.4}, ρ={:+.2}: streamed={:+.12e}, tensor={:+.12e}, \
|Δ|={babs:.3e}, rel={brel:.3e} — the assembled-Gram handoff \
changed the inner solution (EDF/κ-optimum would diverge)",
rho[0],
beta_s[j],
beta_t[j],
);
}
}
}
eprintln!(
"[psi-gram-tensor e2e] coefficient bit-tightness: worst |Δβ̂|={worst_beta_abs:.3e} \
over {} (ρ,ψ) window points — converged inner solution (⇒ EDF, κ-optimum) \
is lane-invariant end-to-end",
rho_sweep.len() * psi_sweep.len(),
);
}
#[test]
fn psi_gram_tensor_fast_path_skips_n_row_lane_and_matches_streamed() {
use gam_model_api::OuterEvalOrder;
let n = 600usize;
let mut data = Array2::<f64>::zeros((n, 1));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
let signal = 1.2 * (2.0 * std::f64::consts::PI * t).sin() + 0.4 * (t - 0.5);
let noise = 0.15 * (((i as f64) * 12.9898).sin() * 43758.547).fract();
y[i] = signal + noise;
}
let weights = Array1::<f64>::ones(n);
let offset = Array1::<f64>::zeros(n);
let family = LikelihoodSpec::gaussian_identity();
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "fast_path_skip".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 12 },
length_scale: Some(1.0),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::all_active(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let frozen = freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze", e));
let frozen_design = build_term_collection_design(data.view(), &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "frozen design", e));
let spatial_terms = spatial_length_scale_term_indices(&frozen);
let dims_per_term = spatial_dims_per_term(&frozen, &spatial_terms);
let rho_dim = frozen_design.penalties.len();
let kappa_options = SpatialLengthScaleOptimizationOptions::default();
let log_kappa0 =
SpatialLogKappaCoords::from_length_scales(&frozen, &spatial_terms, &kappa_options);
let log_kappa_lower = SpatialLogKappaCoords::lower_bounds_from_data(
data.view(),
&frozen,
&spatial_terms,
&kappa_options,
);
let log_kappa_upper = SpatialLogKappaCoords::upper_bounds_from_data(
data.view(),
&frozen,
&spatial_terms,
&kappa_options,
);
let log_kappa0 = log_kappa0.clamp_to_bounds(&log_kappa_lower, &log_kappa_upper);
const JOINT_RHO_BOUND: f64 = 12.0;
let setup = ExactJointHyperSetup::new(
Array1::<f64>::zeros(rho_dim),
Array1::<f64>::from_elem(rho_dim, -JOINT_RHO_BOUND),
Array1::<f64>::from_elem(rho_dim, JOINT_RHO_BOUND),
log_kappa0.clone(),
log_kappa_lower.clone(),
log_kappa_upper.clone(),
);
let theta0 = setup.theta0();
let lower = setup.lower();
let upper = setup.upper();
let psi_lo = lower[rho_dim];
let psi_hi = upper[rho_dim];
let z = Array1::from_iter(y.iter().zip(offset.iter()).map(|(yi, oi)| yi - oi));
let external_opts = external_opts_for_design(
&family,
&frozen_design,
&FitOptions {
compute_inference: false,
max_iter: 200,
tol: 1e-12,
penalty_shrinkage_floor: None,
..FitOptions::default()
},
);
let make_eval = || {
gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&frozen_design.design,
offset.view(),
&frozen_design.penalties,
&external_opts,
"fast_path_skip",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "evaluator", e))
};
let make_cache = || {
SingleBlockExactJointDesignCache::new(
data.view(),
frozen.clone(),
frozen_design.clone(),
spatial_terms.clone(),
rho_dim,
dims_per_term.clone(),
)
.expect("design cache")
};
let mut tensor_eval = make_eval();
let mut tensor_cache = make_cache();
let attached = {
let mut build_cache = make_cache();
let theta_probe_base = theta0.clone();
tensor_eval.build_and_set_psi_gram_tensor(
|psi| {
let mut theta_probe = theta_probe_base.clone();
theta_probe[rho_dim] = psi;
build_cache.ensure_theta(&theta_probe)?;
Ok(build_cache.design().design.clone())
},
weights.view(),
z.view(),
psi_lo,
psi_hi,
)
};
assert!(
attached,
"tensor must certify on this fixture for a non-vacuous gate"
);
let nfree_penalty = tensor_cache.supports_nfree_penalty_rekey();
assert!(
nfree_penalty,
"single frozen Duchon term must admit the exact n-free S(ψ) re-key"
);
tensor_eval.set_supports_nfree_penalty_rekey(nfree_penalty);
let span = psi_hi - psi_lo;
let psi_a = psi_lo + 0.5 * span;
let psi_b = psi_a + span / 4096.0;
let psi_c = psi_lo + 0.02 * span;
assert!(
tensor_eval.psi_gram_tensor_covers(psi_a)
&& tensor_eval.psi_gram_tensor_covers(psi_b)
&& tensor_eval.psi_gram_tensor_covers(psi_c),
"all probes must lie inside the certified value window so the #1033 \
production skip can stay on the n-free lane"
);
let rho = Array1::<f64>::from_elem(rho_dim, 0.5);
let theta_at = |psi: f64| {
let mut theta = Array1::<f64>::zeros(rho_dim + 1);
theta.slice_mut(s![..rho_dim]).assign(&rho);
theta[rho_dim] = psi;
theta
};
let eval_tensor = |evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>,
cache: &mut SingleBlockExactJointDesignCache<'_>,
theta: &Array1<f64>,
realize: bool|
-> (f64, Array1<f64>, Array1<f64>) {
if realize {
cache.ensure_theta(theta).unwrap_or_else(|e| panic!("{} failed: {:?}", "ensure_theta", e));
}
let hyper_dirs = cache
.hyper_dirs_for_current_design(data.view(), SpatialHyperKind::Isotropic)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "hyper_dirs", e));
let penalty = cache
.canonical_penalties_at(theta)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "exact n-free S(ψ) rebuild", e));
evaluator.stage_fast_path_penalty(Some(penalty));
let (cost, grad, _h) = evaluate_joint_reml_outer_eval_at_theta(
evaluator,
cache.design(),
theta,
rho_dim,
hyper_dirs,
None,
OuterEvalOrder::ValueAndGradient,
Some(cache.design_revision()),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "tensor eval", e));
let beta = evaluator
.current_beta()
.unwrap_or_else(|| panic!("{} failed", "converged inner β̂ available"));
(cost, grad, beta)
};
assert_eq!(
tensor_eval.slow_path_reset_count(),
0,
"no slow-path reset before the first eval"
);
let (c_a, _g_a, beta_a) =
eval_tensor(&mut tensor_eval, &mut tensor_cache, &theta_at(psi_a), true);
assert_eq!(
tensor_eval.slow_path_reset_count(),
1,
"first eval at a fresh revision must take the slow path exactly once"
);
let covers_b = tensor_eval.psi_gram_tensor_covers(psi_b);
let covers_c = tensor_eval.psi_gram_tensor_covers(psi_c);
let skip_b = tensor_eval.psi_gram_tensor_covers_skip(psi_b);
let skip_c = tensor_eval.psi_gram_tensor_covers_skip(psi_c);
eprintln!(
"[DIAG1264-FP] pinned ψ_a={psi_a:.5} ψ_b={psi_b:.5} value_covers={covers_b} \
covers_skip={skip_b} ψ_c={psi_c:.5} value_covers={covers_c} \
covers_skip={skip_c}"
);
assert!(
covers_b && covers_c,
"both moved trials must lie inside the certified value window (the necessary \
condition for the skip gate to even consider firing); ψ_b={psi_b:.6} \
(covers={covers_b}) ψ_c={psi_c:.6} (covers={covers_c})"
);
assert!(
!skip_c,
"ψ_C={psi_c:.6} must GENUINELY rotate the reduced basis so the restored \
#1264 `reduced_basis_equal` precondition REFUSES the skip and routes it to \
the exact β̂-sound path; got covers_skip={skip_c}"
);
let before_b = tensor_eval.slow_path_reset_count();
let (c_b, _g_b, beta_b) =
eval_tensor(&mut tensor_eval, &mut tensor_cache, &theta_at(psi_b), !skip_b);
let after_b = tensor_eval.slow_path_reset_count();
let reset_b = after_b - before_b;
assert_eq!(
reset_b == 0,
skip_b,
"ψ_B fast-path engagement must match the restored `covers_skip` verdict: \
covers_skip={skip_b} but slow_path_resets={reset_b} (0 ⇒ skip fired)"
);
let before_c = tensor_eval.slow_path_reset_count();
let (c_c, _g_c, beta_c) =
eval_tensor(&mut tensor_eval, &mut tensor_cache, &theta_at(psi_c), !skip_c);
let after_c = tensor_eval.slow_path_reset_count();
let reset_c = after_c - before_c;
assert_eq!(
reset_c, 1,
"ψ_C ROTATES the reduced basis (covers_skip REFUSED it), so the restored \
#1264 precondition MUST take the exact slow path (exactly one reset); got \
resets={reset_c}"
);
assert!(
c_a.is_finite() && c_b.is_finite() && c_c.is_finite(),
"all fast/slow path costs must be finite"
);
let mut streamed_eval = make_eval();
let mut stream_cache = make_cache();
let beta_streamed = |evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>,
cache: &mut SingleBlockExactJointDesignCache<'_>,
theta: &Array1<f64>|
-> Array1<f64> {
cache.ensure_theta(theta).unwrap_or_else(|e| panic!("{} failed: {:?}", "ensure_theta", e));
let hyper = try_build_spatial_log_kappa_hyper_dirs(
data.view(),
cache.spec(),
cache.design(),
&spatial_terms,
)
.unwrap()
.unwrap();
evaluate_joint_reml_outer_eval_at_theta(
evaluator,
cache.design(),
theta,
rho_dim,
hyper,
None,
OuterEvalOrder::ValueAndGradient,
Some(cache.design_revision()),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "streamed eval", e));
evaluator.current_beta().unwrap_or_else(|| panic!("{} failed", "streamed β̂"))
};
let mut worst = 0.0_f64;
for (label, theta, beta_tensor) in [
("psi_b", theta_at(psi_b), &beta_b),
("psi_c", theta_at(psi_c), &beta_c),
] {
let beta_slow = beta_streamed(&mut streamed_eval, &mut stream_cache, &theta);
let r = beta_tensor
.iter()
.zip(beta_slow.iter())
.fold(0.0_f64, |a, (f, s)| a.max((f - s).abs() / (1.0 + s.abs())));
eprintln!(
"[DIAG1264-FP] {label} ψ={:.4} β̂rel={r:.3e} issue-bar=1e-6 \
β̂tensor[0]={:+.6e} β̂streamed[0]={:+.6e}",
theta[rho_dim], beta_tensor[0], beta_slow[0]
);
assert_eq!(beta_tensor.len(), beta_slow.len(), "β̂ dim mismatch @ {label}");
for j in 0..beta_tensor.len() {
assert!(
beta_tensor[j].is_finite() && beta_slow[j].is_finite(),
"non-finite β̂[{j}] @ {label}"
);
let babs = (beta_tensor[j] - beta_slow[j]).abs();
let brel = babs / (1.0 + beta_slow[j].abs());
worst = worst.max(babs);
const ISSUE_BETA_BAR: f64 = 1.0e-6;
assert!(
brel <= ISSUE_BETA_BAR,
"β̂[{j}] @ {label} diverges from the streamed exact solve beyond the \
#1264 issue-mandated 1e-6 bar: tensor={:+.12e} streamed={:+.12e} \
|Δ|={babs:.3e} rel={brel:.3e} > 1e-6 — the restored `reduced_basis_equal` \
skip precondition failed to route this ψ to the exact path, so an \
interpolated-Gram κ-amplified κ-optimum shipped",
beta_tensor[j],
beta_slow[j],
);
}
}
assert!(
beta_a.iter().all(|v| v.is_finite()),
"ψ_A pinning β̂ must be finite"
);
eprintln!(
"[psi-gram-tensor #1033/#1264] pairwise reduced-basis-equality witness: \
skip(ψ_b)={skip_b} skip(ψ_c={psi_c:.4} ROTATED)={skip_c}, \
slow_path_reset_count={}, worst |Δβ̂| vs streamed slow path = {worst:.3e} \
(asserted ≤ the #1264 issue-mandated 1e-6 β̂-rel bar across the basis rotation \
the witness refuses)",
tensor_eval.slow_path_reset_count()
);
}
#[test]
fn psi_gram_skip_forced_rotation_beta_error_ladder_diag() {
use gam_model_api::OuterEvalOrder;
let n = 600usize;
let mut data = Array2::<f64>::zeros((n, 1));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
let signal = 1.2 * (2.0 * std::f64::consts::PI * t).sin() + 0.4 * (t - 0.5);
let noise = 0.15 * (((i as f64) * 12.9898).sin() * 43758.547).fract();
y[i] = signal + noise;
}
let weights = Array1::<f64>::ones(n);
let offset = Array1::<f64>::zeros(n);
let family = LikelihoodSpec::gaussian_identity();
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "forced_skip_diag".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 12 },
length_scale: Some(1.0),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::all_active(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None, },
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let frozen = freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze", e));
let frozen_design = build_term_collection_design(data.view(), &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "frozen design", e));
let spatial_terms = spatial_length_scale_term_indices(&frozen);
let dims_per_term = spatial_dims_per_term(&frozen, &spatial_terms);
let rho_dim = frozen_design.penalties.len();
let kappa_options = SpatialLengthScaleOptimizationOptions::default();
let log_kappa0 =
SpatialLogKappaCoords::from_length_scales(&frozen, &spatial_terms, &kappa_options);
let log_kappa_lower = SpatialLogKappaCoords::lower_bounds_from_data(
data.view(),
&frozen,
&spatial_terms,
&kappa_options,
);
let log_kappa_upper = SpatialLogKappaCoords::upper_bounds_from_data(
data.view(),
&frozen,
&spatial_terms,
&kappa_options,
);
let log_kappa0 = log_kappa0.clamp_to_bounds(&log_kappa_lower, &log_kappa_upper);
const JOINT_RHO_BOUND: f64 = 12.0;
let setup = ExactJointHyperSetup::new(
Array1::<f64>::zeros(rho_dim),
Array1::<f64>::from_elem(rho_dim, -JOINT_RHO_BOUND),
Array1::<f64>::from_elem(rho_dim, JOINT_RHO_BOUND),
log_kappa0.clone(),
log_kappa_lower.clone(),
log_kappa_upper.clone(),
);
let theta0 = setup.theta0();
let lower = setup.lower();
let upper = setup.upper();
let psi_lo = lower[rho_dim];
let psi_hi = upper[rho_dim];
let z = Array1::from_iter(y.iter().zip(offset.iter()).map(|(yi, oi)| yi - oi));
let external_opts = external_opts_for_design(
&family,
&frozen_design,
&FitOptions {
compute_inference: false,
max_iter: 200,
tol: 1e-12,
penalty_shrinkage_floor: None,
..FitOptions::default()
},
);
let make_eval = || {
gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&frozen_design.design,
offset.view(),
&frozen_design.penalties,
&external_opts,
"forced_skip_diag",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "evaluator", e))
};
let make_cache = || {
SingleBlockExactJointDesignCache::new(
data.view(),
frozen.clone(),
frozen_design.clone(),
spatial_terms.clone(),
rho_dim,
dims_per_term.clone(),
)
.expect("design cache")
};
let mut tensor_eval = make_eval();
let mut tensor_cache = make_cache();
let attached = {
let mut build_cache = make_cache();
let theta_probe_base = theta0.clone();
tensor_eval.build_and_set_psi_gram_tensor(
|psi| {
let mut theta_probe = theta_probe_base.clone();
theta_probe[rho_dim] = psi;
build_cache.ensure_theta(&theta_probe)?;
Ok(build_cache.design().design.clone())
},
weights.view(),
z.view(),
psi_lo,
psi_hi,
)
};
assert!(attached, "tensor must certify for a non-vacuous diagnostic");
let nfree_penalty = tensor_cache.supports_nfree_penalty_rekey();
assert!(nfree_penalty, "frozen Duchon term must admit exact n-free S(ψ)");
tensor_eval.set_supports_nfree_penalty_rekey(nfree_penalty);
let span = psi_hi - psi_lo;
let psi_a = psi_lo + 0.5 * span; let rho = Array1::<f64>::from_elem(rho_dim, 0.5);
let theta_at = |psi: f64| {
let mut theta = Array1::<f64>::zeros(rho_dim + 1);
theta.slice_mut(s![..rho_dim]).assign(&rho);
theta[rho_dim] = psi;
theta
};
let eval_tensor = |evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>,
cache: &mut SingleBlockExactJointDesignCache<'_>,
theta: &Array1<f64>,
realize: bool|
-> Array1<f64> {
if realize {
cache.ensure_theta(theta).unwrap_or_else(|e| panic!("{} failed: {:?}", "ensure_theta", e));
}
let hyper_dirs = cache
.hyper_dirs_for_current_design(data.view(), SpatialHyperKind::Isotropic)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "hyper_dirs", e));
let penalty = cache
.canonical_penalties_at(theta)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "exact n-free S(ψ) rebuild", e));
evaluator.stage_fast_path_penalty(Some(penalty));
evaluate_joint_reml_outer_eval_at_theta(
evaluator,
cache.design(),
theta,
rho_dim,
hyper_dirs,
None,
OuterEvalOrder::ValueAndGradient,
Some(cache.design_revision()),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "tensor eval", e));
evaluator.current_beta().unwrap_or_else(|| panic!("{} failed", "converged inner β̂"))
};
let beta_a = eval_tensor(&mut tensor_eval, &mut tensor_cache, &theta_at(psi_a), true);
assert!(beta_a.iter().all(|v| v.is_finite()), "pin β̂ finite");
assert_eq!(tensor_eval.slow_path_reset_count(), 1, "one pinning reset");
let mut streamed_eval = make_eval();
let mut stream_cache = make_cache();
let beta_streamed = |evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>,
cache: &mut SingleBlockExactJointDesignCache<'_>,
theta: &Array1<f64>|
-> Array1<f64> {
cache.ensure_theta(theta).unwrap_or_else(|e| panic!("{} failed: {:?}", "ensure_theta", e));
let hyper = try_build_spatial_log_kappa_hyper_dirs(
data.view(),
cache.spec(),
cache.design(),
&spatial_terms,
)
.unwrap()
.unwrap();
evaluate_joint_reml_outer_eval_at_theta(
evaluator,
cache.design(),
theta,
rho_dim,
hyper,
None,
OuterEvalOrder::ValueAndGradient,
Some(cache.design_revision()),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "streamed eval", e));
evaluator.current_beta().unwrap_or_else(|| panic!("{} failed", "streamed β̂"))
};
eprintln!(
"[DIAG1033-FORCE] 513-node tensor; pin ψ_a={psi_a:.5} window=[{psi_lo:.4},{psi_hi:.4}] span={span:.4}"
);
eprintln!("[DIAG1033-FORCE] frac_of_span psi covers_skip(gate) forced_skip_β̂rel_vs_streamed");
let fracs = [
1.0 / 4096.0,
1.0 / 512.0,
1.0 / 64.0,
1.0 / 16.0,
1.0 / 8.0,
1.0 / 4.0,
-0.48, ];
let frob_rel = |a: &Array2<f64>, b: &Array2<f64>| -> f64 {
let num: f64 = a
.iter()
.zip(b.iter())
.map(|(x, y)| (x - y).powi(2))
.sum::<f64>()
.sqrt();
let den: f64 = b.iter().map(|y| y * y).sum::<f64>().sqrt().max(1e-300);
num / den
};
let exact_gram_at = |psi: f64| -> Array2<f64> {
let mut c = make_cache();
let mut th = theta0.clone();
th[rho_dim] = psi;
c.ensure_theta(&th).unwrap_or_else(|e| panic!("{} failed: {:?}", "ensure_theta exact", e));
let x = c.design().design.to_dense();
let mut wx = x.clone();
for (mut row, &wi) in wx.outer_iter_mut().zip(weights.iter()) {
row.mapv_inplace(|v| v * wi);
}
x.t().dot(&wx)
};
let exact_rhs_at = |psi: f64| -> Array1<f64> {
let mut c = make_cache();
let mut th = theta0.clone();
th[rho_dim] = psi;
c.ensure_theta(&th).unwrap_or_else(|e| panic!("{} failed: {:?}", "ensure_theta exact", e));
let x = c.design().design.to_dense();
let mut wz = Array1::<f64>::zeros(x.nrows());
for ((s, &wi), (&yi, &oi)) in wz
.iter_mut()
.zip(weights.iter())
.zip(y.iter().zip(offset.iter()))
{
*s = wi * (yi - oi);
}
x.t().dot(&wz)
};
let vec_rel = |a: &Array1<f64>, b: &Array1<f64>| -> f64 {
let num: f64 = a
.iter()
.zip(b.iter())
.map(|(x, y)| (x - y).powi(2))
.sum::<f64>()
.sqrt();
let den: f64 = b.iter().map(|y| y * y).sum::<f64>().sqrt().max(1e-300);
num / den
};
let tensor_arc = tensor_eval
.psi_gram_tensor_arc()
.unwrap_or_else(|| panic!("{} failed", "tensor attached"));
let tensor_gram_at = |psi: f64| -> Array2<f64> { tensor_arc.gram_at(psi) };
let tensor_rhs_at = |psi: f64| -> Array1<f64> { tensor_arc.rhs_at(psi) };
let g_interp_a = tensor_gram_at(psi_a);
let g_exact_a = exact_gram_at(psi_a);
let r_interp_a = tensor_rhs_at(psi_a);
let r_exact_a = exact_rhs_at(psi_a);
let mut worst_forced = 0.0_f64;
for &f in &fracs {
let psi = psi_a + f * span;
if !tensor_eval.psi_gram_tensor_covers(psi) {
eprintln!("[DIAG1033-FORCE] frac={f:+.5} psi={psi:.5} OUT-OF-VALUE-WINDOW (skipped)");
continue;
}
let gate = tensor_eval.psi_gram_tensor_covers_skip(psi);
let beta_forced =
eval_tensor(&mut tensor_eval, &mut tensor_cache, &theta_at(psi), false);
let beta_slow = beta_streamed(&mut streamed_eval, &mut stream_cache, &theta_at(psi));
let rel = beta_forced
.iter()
.zip(beta_slow.iter())
.fold(0.0_f64, |a, (g, s)| a.max((g - s).abs() / (1.0 + s.abs())));
worst_forced = worst_forced.max(rel);
let g_interp = tensor_gram_at(psi);
let g_exact = exact_gram_at(psi);
let interp_res = frob_rel(&g_interp, &g_exact);
let mut g_anchored = g_interp.clone();
g_anchored += &g_exact_a;
g_anchored -= &g_interp_a;
let anchored_res = frob_rel(&g_anchored, &g_exact);
let r_interp = tensor_rhs_at(psi);
let r_exact = exact_rhs_at(psi);
let rhs_interp_res = vec_rel(&r_interp, &r_exact);
let mut r_anchored = r_interp.clone();
r_anchored += &r_exact_a;
r_anchored -= &r_interp_a;
let rhs_anchored_res = vec_rel(&r_anchored, &r_exact);
let p_cols = stream_cache.design().design.ncols();
let to_dense_canonical =
|pens: &[gam_terms::construction::CanonicalPenalty]| -> Array2<f64> {
let mut s = Array2::<f64>::zeros((p_cols, p_cols));
for cp in pens {
let r = &cp.col_range;
s.slice_mut(ndarray::s![r.start..r.end, r.start..r.end])
.scaled_add(1.0, &cp.local);
}
s
};
let (staged_pens, _) = tensor_cache
.canonical_penalties_at(&theta_at(psi))
.unwrap_or_else(|e| panic!("{} failed: {:?}", "staged canonical penalties", e));
let s_staged = to_dense_canonical(&staged_pens);
let mut s_streamed = Array2::<f64>::zeros((p_cols, p_cols));
for bp in &stream_cache.design().penalties {
let r = &bp.col_range;
s_streamed
.slice_mut(ndarray::s![r.start..r.end, r.start..r.end])
.scaled_add(1.0, &bp.local);
}
let pen_res = frob_rel(&s_staged, &s_streamed);
let pen_blocks_staged = staged_pens.len();
let pen_blocks_streamed = stream_cache.design().penalties.len();
eprintln!(
"[DIAG1033-FORCE] frac={f:+.5} psi={psi:.5} covers_skip={gate} forced_β̂rel={rel:.3e} \
gram_interp_res={interp_res:.3e} gram_anchored_res={anchored_res:.3e} \
rhs_interp_res={rhs_interp_res:.3e} rhs_anchored_res={rhs_anchored_res:.3e} \
pen_res={pen_res:.3e} pen_blocks={pen_blocks_staged}/{pen_blocks_streamed} {}",
if rel <= 1e-6 { "SOUND(≤1e-6)" } else { "UNSOUND(>1e-6)" }
);
drop(eval_tensor(&mut tensor_eval, &mut tensor_cache, &theta_at(psi_a), true));
}
eprintln!(
"[DIAG1033-FORCE] worst forced-skip β̂rel across the move ladder = {worst_forced:.3e} \
(≤1e-6 ⇒ gate over-conservative, replace with residual self-check; \
growing-to-~1e-5 ⇒ interpolation×conditioning wall, needs deeper tensor)"
);
assert!(
worst_forced.is_finite(),
"forced-skip β̂rel ladder must be finite to be a usable measurement"
);
}
#[test]
fn iso_kappa_duchon_binomial_logit_fd() {
let (pass, worst, violations) =
iso_kappa_fd_variant_driver("duchon_logit", 80, LikelihoodSpec::binomial_logit(), false, false);
assert!(
pass,
"BinomialLogit FD failed; worst_psi_rel={worst:.3e}\n {}",
violations.join("\n ")
);
}
#[test]
fn iso_kappa_duchon_n_smaller_fd() {
let (pass, worst, violations) = iso_kappa_fd_variant_driver(
"duchon_probit_n20",
20,
LikelihoodSpec::binomial_probit(),
false,
true,
);
assert!(
pass,
"Duchon Probit n=20 FD failed; worst_psi_rel={worst:.3e}\n {}",
violations.join("\n ")
);
}
#[test]
fn iso_kappa_duchon_no_psi_fd() {
let (pass, _worst, violations) = iso_kappa_fd_variant_driver(
"duchon_probit_rho_only",
80,
LikelihoodSpec::binomial_probit(),
true,
false,
);
assert!(
pass,
"Duchon Probit ρ-only FD failed:\n {}",
violations.join("\n ")
);
}
struct DuchonProbitSetup {
data: Array2<f64>,
y: Array1<f64>,
weights: Array1<f64>,
offset: Array1<f64>,
frozen: TermCollectionSpec,
frozen_design: TermCollectionDesign,
spatial_terms: Vec<usize>,
dims_per_term: Vec<usize>,
rho_dim: usize,
psi_dim: usize,
}
fn build_duchon_probit_setup() -> DuchonProbitSetup {
let n = 80usize;
let mut data = Array2::<f64>::zeros((n, 1));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
let eta = 1.4 * (2.0 * std::f64::consts::PI * t).sin() + 0.5 * (t - 0.5);
y[i] = if eta + 0.7 * (3.7 * (i as f64) + 1.0).sin() > 0.0 {
1.0
} else {
0.0
};
}
let weights = Array1::ones(n);
let offset = Array1::zeros(n);
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "duchon_1d".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
length_scale: Some(1.0),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::all_active(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let frozen = freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze", e));
let frozen_design = build_term_collection_design(data.view(), &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "frozen design", e));
let spatial_terms = spatial_length_scale_term_indices(&frozen);
let dims_per_term = spatial_dims_per_term(&frozen, &spatial_terms);
let rho_dim = frozen_design.penalties.len();
let psi_dim: usize = dims_per_term.iter().sum();
DuchonProbitSetup {
data,
y,
weights,
offset,
frozen,
frozen_design,
spatial_terms,
dims_per_term,
rho_dim,
psi_dim,
}
}
#[test]
fn iso_kappa_duchon_outer_gradient_matches_centered_fd() {
let DuchonProbitSetup {
data,
y,
weights,
offset,
frozen,
frozen_design,
spatial_terms,
dims_per_term,
rho_dim,
psi_dim,
} = build_duchon_probit_setup();
let fit_opts = FitOptions {
compute_inference: false,
max_iter: 200,
tol: 1e-12,
penalty_shrinkage_floor: None,
..FitOptions::default()
};
let external_opts = external_opts_for_design(
&LikelihoodSpec::binomial_probit(),
&frozen_design,
&fit_opts,
);
let mut cache = SingleBlockExactJointDesignCache::new(
data.view(),
frozen.clone(),
frozen_design.clone(),
spatial_terms.clone(),
rho_dim,
dims_per_term.clone(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "cache", e));
let mut evaluator = gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&frozen_design.design,
offset.view(),
&frozen_design.penalties,
&external_opts,
"iso-kappa Duchon gradient FD pin",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "evaluator", e));
let theta_dim = rho_dim + psi_dim;
let theta_zero = Array1::<f64>::zeros(theta_dim);
let eval_at =
|theta: &Array1<f64>,
order: gam_model_api::OuterEvalOrder,
cache: &mut SingleBlockExactJointDesignCache<'_>,
evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>| {
cache.ensure_theta(theta).unwrap_or_else(|e| panic!("{} failed: {:?}", "ensure_theta", e));
let hyper_dirs = try_build_spatial_log_kappa_hyper_dirs(
data.view(),
cache.spec(),
cache.design(),
&cache.spatial_terms,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "hyper dirs build", e))
.expect("hyper dirs present");
evaluate_joint_reml_outer_eval_at_theta(
evaluator,
cache.design(),
theta,
rho_dim,
hyper_dirs,
None,
order,
None,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "outer eval", e))
};
let (cost_at_zero, grad_at_zero, _hess) = eval_at(
&theta_zero,
gam_model_api::OuterEvalOrder::ValueAndGradient,
&mut cache,
&mut evaluator,
);
let h = 1e-5_f64;
let psi_idx = rho_dim;
let mut theta_p = theta_zero.clone();
theta_p[psi_idx] += h;
let mut theta_m = theta_zero.clone();
theta_m[psi_idx] -= h;
let (cost_p, _, _) = eval_at(
&theta_p,
gam_model_api::OuterEvalOrder::Value,
&mut cache,
&mut evaluator,
);
let (cost_m, _, _) = eval_at(
&theta_m,
gam_model_api::OuterEvalOrder::Value,
&mut cache,
&mut evaluator,
);
let fd_psi_gradient = (cost_p - cost_m) / (2.0 * h);
let analytic_psi_gradient = grad_at_zero[psi_idx];
let scale = 1.0 + analytic_psi_gradient.abs().max(fd_psi_gradient.abs());
let rel = (analytic_psi_gradient - fd_psi_gradient).abs() / scale;
assert!(
rel < 1e-3,
"Duchon ψ outer gradient must match centered FD of the production objective: \
analytic={:+.4e}, fd={:+.4e}, rel={:+.3e}",
analytic_psi_gradient,
fd_psi_gradient,
rel
);
assert!(
cost_at_zero.is_finite() && grad_at_zero.iter().all(|v| v.is_finite()),
"ψ-gradient and cost must be finite at θ=0"
);
}
#[test]
fn duchon_probit_pirls_determinism_at_zero() {
let DuchonProbitSetup {
data,
y,
weights,
offset,
frozen,
frozen_design,
spatial_terms,
dims_per_term,
rho_dim,
psi_dim,
} = build_duchon_probit_setup();
let fit_opts = FitOptions {
compute_inference: false,
max_iter: 200,
tol: 1e-6,
penalty_shrinkage_floor: None,
..FitOptions::default()
};
let external_opts = external_opts_for_design(
&LikelihoodSpec::binomial_probit(),
&frozen_design,
&fit_opts,
);
let mut cache = SingleBlockExactJointDesignCache::new(
data.view(),
frozen.clone(),
frozen_design.clone(),
spatial_terms.clone(),
rho_dim,
dims_per_term.clone(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "cache", e));
let mut evaluator = gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&frozen_design.design,
offset.view(),
&frozen_design.penalties,
&external_opts,
"PIRLS-determinism",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "evaluator", e));
let theta_dim = rho_dim + psi_dim;
let theta_zero = Array1::<f64>::zeros(theta_dim);
let mut h_calls = Vec::new();
for trial in 0..3 {
cache.ensure_theta(&theta_zero).unwrap_or_else(|e| panic!("{} failed: {:?}", "ensure_theta", e));
let d = cache.design().clone();
let h_i = evaluator
.debug_full_h(
&d.design,
&d.penalties,
&d.nullspace_dims,
d.linear_constraints.clone(),
&theta_zero,
rho_dim,
&format!("determinism trial {}", trial),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "debug_full_h", e));
h_calls.push(h_i);
}
let h0 = &h_calls[0];
let norm_h0 = h0.iter().map(|v| v.abs()).fold(0.0_f64, f64::max);
let abs_tol = (1e-5_f64) * norm_h0 + 1e-12_f64;
for trial in 1..3 {
let diff = &h_calls[trial] - h0;
let max_abs = diff.iter().map(|v| v.abs()).fold(0.0_f64, f64::max);
assert!(
max_abs <= abs_tol,
"PIRLS non-deterministic at fixed θ (trial {trial} vs 0): \
max|Δ|={max_abs:+.6e} > tol={abs_tol:+.6e} \
(‖H‖_∞={norm_h0:+.6e})"
);
}
}
#[test]
fn spatial_aniso_joint_large_psi_dim_keeps_second_order_route() {
let cap = gam_solve::rho_optimizer::OuterCapability {
gradient: gam_problem::Derivative::Analytic,
hessian: gam_problem::DeclaredHessianForm::Either,
n_params: 40,
psi_dim: 31,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let route = gam_solve::rho_optimizer::plan(&cap);
assert_eq!(route.solver, gam_solve::rho_optimizer::Solver::Arc);
assert_eq!(
route.hessian_source,
gam_solve::rho_optimizer::HessianSource::Analytic
);
assert!(route.routing_log_line().contains("matrix-free=false"));
}
#[test]
fn exact_joint_spatial_outer_hessian_available_for_sparse_designs() {
let n = 96usize;
let x = Array1::linspace(0.0, 1.0, n);
let mut data = Array2::<f64>::zeros((n, 1));
data.column_mut(0).assign(&x);
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "s_x".to_string(),
basis: SmoothBasisSpec::BSpline1D {
feature_col: 0,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 32,
},
double_penalty: true,
identifiability: BSplineIdentifiability::None,
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
assert!(matches!(design.design, DesignMatrix::Sparse(_)));
assert!(exact_joint_spatial_outer_hessian_available(
&LikelihoodSpec::binomial_logit(),
&design,
));
}
#[test]
fn iso_kappa_duchon_dx_dpsi_matches_fd() {
let n = 80usize;
let mut data = Array2::<f64>::zeros((n, 1));
for i in 0..n {
data[[i, 0]] = i as f64 / (n as f64 - 1.0);
}
let spec_orig = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "duchon_1d".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
length_scale: Some(1.0),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec_orig).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let frozen = freeze_term_collection_from_design(&spec_orig, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze", e));
let build_design_at = |psi: f64| -> Array2<f64> {
let mut s = frozen.clone();
if let SmoothBasisSpec::Duchon {
spec: ref mut duchon,
..
} = s.smooth_terms[0].basis
{
duchon.length_scale = Some((-psi).exp());
}
let d = build_term_collection_design(data.view(), &s).unwrap_or_else(|e| panic!("{} failed: {:?}", "rebuild", e));
d.design.to_dense()
};
let psi_eval = 0.0_f64;
let derivative_bundle =
try_build_spatial_term_log_kappa_derivative(data.view(), &frozen, &design, 0)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "formula Duchon derivative should build", e))
.expect("Duchon derivative should be available");
let global_range = derivative_bundle.0;
let p_total = derivative_bundle.1;
let implicit_operator = derivative_bundle.8;
let op = implicit_operator.unwrap_or_else(|| panic!("{} failed", "Duchon derivative should expose implicit operator"));
let p = op.p_out();
assert_eq!(p_total, design.design.ncols());
assert_eq!(global_range.end - global_range.start, p);
let h = 1e-4_f64;
let x_plus = build_design_at(psi_eval + h);
let x_minus = build_design_at(psi_eval - h);
eprintln!(
"[DXDPSI_FD] X(+h)[0,0..3]={:?} X(-h)[0,0..3]={:?}",
x_plus.row(0).iter().take(3).copied().collect::<Vec<_>>(),
x_minus.row(0).iter().take(3).copied().collect::<Vec<_>>(),
);
eprintln!(
"[DXDPSI_FD] X(+h) shape={:?} X(-h) shape={:?} p_out={}",
x_plus.shape(),
x_minus.shape(),
p,
);
let x_at = build_design_at(psi_eval);
let orig_design = build_term_collection_design(data.view(), &spec_orig).unwrap_or_else(|e| panic!("{} failed: {:?}", "rebuild orig", e));
eprintln!(
"[DXDPSI_FD] X(psi_eval) shape={:?} orig_design.ncols={}",
x_at.shape(),
orig_design.design.ncols(),
);
let mut analytic = Array2::<f64>::zeros((n, p));
let mut basisv = Array1::<f64>::zeros(p);
for j in 0..p {
basisv[j] = 1.0;
let col = op.forward_mul(0, &basisv.view()).unwrap_or_else(|e| panic!("{} failed: {:?}", "forward_mul", e));
analytic.column_mut(j).assign(&col);
basisv[j] = 0.0;
}
let smooth_start = global_range.start;
let v_test = Array1::<f64>::from_shape_fn(n, |i| (i as f64 * 0.07).sin());
let analytic_tv = op.transpose_mul(0, &v_test.view()).unwrap_or_else(|e| panic!("{} failed: {:?}", "transpose_mul", e));
let fd_tv_full = (&x_plus.t() - &x_minus.t()) / (2.0 * h);
let fd_tv = fd_tv_full.dot(&v_test);
let fd_tv_smooth = fd_tv.slice(s![smooth_start..(smooth_start + p)]).to_owned();
let max_tv_diff = analytic_tv
.iter()
.zip(fd_tv_smooth.iter())
.map(|(a, b)| (a - b).abs())
.fold(0.0f64, f64::max);
let max_tv_abs = analytic_tv.iter().map(|v| v.abs()).fold(0.0f64, f64::max);
eprintln!(
"[DXDPSI_TV] max|analytic_tv - fd_tv|={:.3e} max|analytic_tv|={:.3e}",
max_tv_diff, max_tv_abs
);
eprintln!(
"[DXDPSI_TV] analytic_tv={:?}",
analytic_tv.iter().take(p).copied().collect::<Vec<_>>()
);
eprintln!(
"[DXDPSI_TV] fd_tv_smooth={:?}",
fd_tv_smooth.iter().take(p).copied().collect::<Vec<_>>()
);
let fd_full = (&x_plus - &x_minus) / (2.0 * h);
let fd = fd_full
.slice(s![.., smooth_start..(smooth_start + p)])
.to_owned();
let mut max_diff = 0.0_f64;
let mut max_abs = 0.0_f64;
for i in 0..n {
for j in 0..p {
let d = (analytic[[i, j]] - fd[[i, j]]).abs();
if d > max_diff {
max_diff = d;
}
if analytic[[i, j]].abs() > max_abs {
max_abs = analytic[[i, j]].abs();
}
}
}
eprintln!(
"[DXDPSI_FD] max|analytic - fd|={:.3e} max|analytic|={:.3e}",
max_diff, max_abs
);
eprintln!(
"[DXDPSI_FD] analytic[0,..]={:?}",
analytic.row(0).iter().take(p).copied().collect::<Vec<_>>(),
);
eprintln!(
"[DXDPSI_FD] fd[0,..]={:?}",
fd.row(0).iter().take(p).copied().collect::<Vec<_>>(),
);
assert!(max_diff < 5e-3 * max_abs.max(1e-3), "dX/dψ mismatch");
}
#[test]
fn joint_build_and_freeze_shares_auto_spatial_centers_across_blocks() {
let n = 400usize;
let mut data = Array2::<f64>::zeros((n, 2));
for i in 0..n {
data[[i, 0]] = i as f64 / (n as f64 - 1.0);
data[[i, 1]] = (i as f64 * 0.19).sin();
}
let matern_term = |name: &str| SmoothTermSpec {
name: name.to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::Auto(Box::new(CenterStrategy::FarthestPoint {
num_centers: 8,
})),
length_scale: 0.8,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let marginalspec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![matern_term("marginal")],
};
let logslopespec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![matern_term("logslope")],
};
let (designs, resolved_specs) = build_term_collection_designs_and_freeze_joint(
data.view(),
&[marginalspec.clone(), logslopespec.clone()],
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "joint build and freeze should succeed", e));
assert_eq!(designs.len(), 2);
assert_eq!(resolved_specs.len(), 2);
let extract_centers = |spec: &TermCollectionSpec| match &spec.smooth_terms[0].basis {
SmoothBasisSpec::Matern { spec, .. } => match &spec.center_strategy {
CenterStrategy::UserProvided(centers) => centers.clone(),
other => panic!("expected frozen user-provided centers, got {other:?}"),
},
other => panic!("expected Matérn term, got {other:?}"),
};
let marginal_centers = extract_centers(&resolved_specs[0]);
let logslope_centers = extract_centers(&resolved_specs[1]);
let separate_marginal_design =
build_term_collection_design(data.view(), &marginalspec).unwrap_or_else(|e| panic!("{} failed: {:?}", "separate marginal", e));
let separate_marginal =
freeze_term_collection_from_design(&marginalspec, &separate_marginal_design)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze separate marginal", e));
let separate_marginal_centers = extract_centers(&separate_marginal);
assert_eq!(marginal_centers, logslope_centers);
assert_eq!(marginal_centers.ncols(), 2);
assert_eq!(marginal_centers.nrows(), separate_marginal_centers.nrows());
}
#[test]
fn exact_joint_two_block_no_spatial_fast_path_returns_fully_frozen_specs() {
let n = 24usize;
let mut data = Array2::<f64>::zeros((n, 2));
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
data[[i, 1]] = (i % 3) as f64;
}
let pspline_term = |name: &str| SmoothTermSpec {
name: name.to_string(),
basis: SmoothBasisSpec::BSpline1D {
feature_col: 0,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 3,
},
double_penalty: true,
identifiability: BSplineIdentifiability::None,
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let random_effect = RandomEffectTermSpec {
name: "grp".to_string(),
feature_col: 1,
drop_first_level: false,
penalized: true,
frozen_levels: None,
};
let meanspec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![random_effect.clone()],
smooth_terms: vec![pspline_term("mean_ps")],
};
let noisespec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![random_effect],
smooth_terms: vec![pspline_term("noise_ps")],
};
let kappa_options = SpatialLengthScaleOptimizationOptions {
max_outer_iter: 1,
rel_tol: 1e-6,
pilot_subsample_threshold: 0,
outer_wall_clock_budget_secs: None,
..SpatialLengthScaleOptimizationOptions::default()
};
let joint_setup = two_block_exact_joint_hyper_setup(&meanspec, &noisespec, &kappa_options);
let theta_dim = joint_setup.theta0().len();
let mean_terms = spatial_length_scale_term_indices(&meanspec);
let noise_terms = spatial_length_scale_term_indices(&noisespec);
assert!(mean_terms.is_empty());
assert!(noise_terms.is_empty());
let policy = gam_model_api::families::custom_family::OuterDerivativePolicy {
capability: gam_problem::ExactOuterDerivativeOrder::Second,
predicted_hessian_work: 0,
predicted_gradient_work: 0,
subsample_capable: false,
};
let solved = optimize_spatial_length_scale_exact_joint(
data.view(),
&[meanspec.clone(), noisespec.clone()],
&[mean_terms, noise_terms],
&kappa_options,
&joint_setup,
gam_problem::SeedRiskProfile::Gaussian,
true,
true,
false,
None,
policy,
|theta, specs, designs| {
assert_eq!(theta.len(), theta_dim);
assert_eq!(specs.len(), 2);
assert_eq!(designs.len(), 2);
Ok(designs[0].design.ncols() as f64 + designs[1].design.ncols() as f64)
},
|theta, specs, designs, eval_mode, _| {
assert_eq!(theta.len(), theta_dim);
assert_eq!(specs.len(), 2);
assert_eq!(designs.len(), 2);
Ok((
0.0,
Array1::zeros(theta_dim),
if matches!(
eval_mode,
gam_solve::estimate::reml::reml_outer_engine::EvalMode::ValueGradientHessian
) {
gam_problem::HessianResult::Analytic(Array2::zeros((
theta_dim, theta_dim,
)))
} else {
gam_problem::HessianResult::Unavailable
},
))
},
|theta, specs, designs| {
assert_eq!(theta.len(), theta_dim);
assert_eq!(specs.len(), 2);
assert_eq!(designs.len(), 2);
Ok(gam_problem::EfsEval {
cost: 0.0,
steps: vec![0.0; theta_dim],
beta: None,
psi_gradient: None,
psi_indices: None,
inner_hessian_scale: None,
logdet_enclosure_gap: None,
})
},
|_beta: &Array1<f64>| Ok(gam_solve::rho_optimizer::SeedOutcome::NoSlot),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "exact joint no-spatial fast path should succeed", e));
for resolved in [&solved.resolved_specs[0], &solved.resolved_specs[1]] {
resolved
.validate_frozen("resolvedspec")
.unwrap_or_else(|e| panic!("{} failed: {:?}", "exact joint no-spatial fast path should fully freeze specs", e));
match &resolved.smooth_terms[0].basis {
SmoothBasisSpec::BSpline1D { spec, .. } => {
assert!(matches!(spec.knotspec, BSplineKnotSpec::Provided(_)));
}
_ => panic!("expected P-spline term"),
}
assert!(
resolved.random_effect_terms[0].frozen_levels.is_some(),
"random-effect levels should be frozen in exact joint no-spatial fast path"
);
}
}
#[test]
fn incremental_frozen_realizer_matches_unified_full_rebuild() {
let n = 24usize;
let mut data = Array2::<f64>::zeros((n, 4));
for i in 0..n {
let t = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = t;
data[[i, 1]] = (0.35 * i as f64).sin();
data[[i, 2]] = (i % 3) as f64;
data[[i, 3]] = t * t;
}
let spec = TermCollectionSpec {
linear_terms: vec![LinearTermSpec {
name: "lin".to_string(),
feature_col: 1,
feature_cols: vec![1],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: Some(-0.5),
coefficient_max: None,
}],
random_effect_terms: vec![RandomEffectTermSpec {
name: "grp".to_string(),
feature_col: 2,
drop_first_level: false,
penalized: true,
frozen_levels: None,
}],
smooth_terms: vec![
SmoothTermSpec {
name: "spatial".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
length_scale: 0.8,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: Some(vec![0.15, -0.15]),
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
SmoothTermSpec {
name: "mono".to_string(),
basis: SmoothBasisSpec::BSpline1D {
feature_col: 3,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 3,
},
double_penalty: false,
identifiability: BSplineIdentifiability::None,
boundary: OneDimensionalBoundary::Open,
boundary_conditions: BSplineBoundaryConditions::default(),
},
},
shape: ShapeConstraint::MonotoneIncreasing,
joint_null_rotation: None,
},
],
};
let base_design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "base design", e));
let frozen = freeze_term_collection_from_design(&spec, &base_design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze", e));
let frozen_design = build_term_collection_design(data.view(), &frozen).unwrap_or_else(|e| panic!("{} failed: {:?}", "frozen design", e));
let spatial_terms = spatial_length_scale_term_indices(&frozen);
assert_eq!(spatial_terms, vec![0]);
let smooth_start = frozen_design.design.ncols() - frozen_design.smooth.total_smooth_cols();
let fixed_before = frozen_design.design.clone();
let nonspatial_range = frozen_design.smooth.terms[1].coeff_range.clone();
let full_nonspatial_range =
(smooth_start + nonspatial_range.start)..(smooth_start + nonspatial_range.end);
let mut realizer = FrozenTermCollectionIncrementalRealizer::new(
data.view(),
frozen.clone(),
frozen_design.clone(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "incremental realizer", e));
let updated_log_kappa = SpatialLogKappaCoords::new_with_dims(array![0.30, -0.20], vec![2]);
let updated_spec = updated_log_kappa
.apply_tospec(&frozen, &spatial_terms)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "updated spec", e));
realizer
.apply_log_kappa(&updated_log_kappa, &spatial_terms)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "incremental update", e));
let rebuilt = build_term_collection_design(data.view(), &updated_spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "rebuilt design", e));
assert_term_collection_designs_match(realizer.design(), &rebuilt, "incremental realizer");
let linear_range = frozen_design.linear_ranges[0].1.clone();
let random_range = frozen_design.random_effect_ranges[0].1.clone();
let fixed_before_dense = fixed_before.to_dense();
let updated_full_dense = realizer.design().design.to_dense();
let linear_diff = max_abs_diff_matrix(
&fixed_before_dense
.slice(s![.., linear_range.clone()])
.to_owned(),
&updated_full_dense.slice(s![.., linear_range]).to_owned(),
);
let random_diff = max_abs_diff_matrix(
&fixed_before_dense
.slice(s![.., random_range.clone()])
.to_owned(),
&updated_full_dense.slice(s![.., random_range]).to_owned(),
);
let nonspatial_diff = max_abs_diff_matrix(
&fixed_before_dense
.slice(s![.., full_nonspatial_range.clone()])
.to_owned(),
&updated_full_dense
.slice(s![.., full_nonspatial_range.clone()])
.to_owned(),
);
let spatial_range = frozen_design.smooth.terms[0].coeff_range.clone();
let full_spatial_range =
(smooth_start + spatial_range.start)..(smooth_start + spatial_range.end);
let spatial_change = max_abs_diff_matrix(
&fixed_before_dense
.slice(s![.., full_spatial_range.clone()])
.to_owned(),
&updated_full_dense
.slice(s![.., full_spatial_range])
.to_owned(),
);
assert!(
linear_diff <= 1e-12,
"linear block changed max_abs={linear_diff}"
);
assert!(
random_diff <= 1e-12,
"random-effect block changed max_abs={random_diff}"
);
assert!(
nonspatial_diff <= 1e-12,
"unchanged smooth block changed max_abs={nonspatial_diff}"
);
assert!(
spatial_change > 1e-8,
"spatial block did not update max_abs={spatial_change}"
);
}
#[test]
fn two_block_exact_joint_design_cache_clears_memo_on_theta_change() {
let n = 20usize;
let mut data = Array2::<f64>::zeros((n, 2));
for i in 0..n {
let x0 = i as f64 / (n as f64 - 1.0);
let x1 = (0.19 * i as f64).sin();
data[[i, 0]] = x0;
data[[i, 1]] = x1;
}
let matern_term = |name: &str, length_scale: f64| SmoothTermSpec {
name: name.to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 5 },
length_scale,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: Some(vec![0.0, 0.0]),
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
};
let meanspec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![matern_term("mean", 0.7)],
};
let noisespec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![matern_term("noise", 1.1)],
};
let kappa_options = SpatialLengthScaleOptimizationOptions {
max_outer_iter: 1,
rel_tol: 1e-6,
pilot_subsample_threshold: 0,
outer_wall_clock_budget_secs: None,
..SpatialLengthScaleOptimizationOptions::default()
};
let joint_setup = two_block_exact_joint_hyper_setup(&meanspec, &noisespec, &kappa_options);
let theta0 = joint_setup.theta0();
let mean_design = build_term_collection_design(data.view(), &meanspec).unwrap_or_else(|e| panic!("{} failed: {:?}", "mean", e));
let noise_design = build_term_collection_design(data.view(), &noisespec).unwrap_or_else(|e| panic!("{} failed: {:?}", "noise", e));
let mean_frozen =
freeze_term_collection_from_design(&meanspec, &mean_design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze mean", e));
let noise_frozen =
freeze_term_collection_from_design(&noisespec, &noise_design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze noise", e));
let mean_term_indices = spatial_length_scale_term_indices(&mean_frozen);
let noise_term_indices = spatial_length_scale_term_indices(&noise_frozen);
let mut cache = ExactJointDesignCache::new(
data.view(),
vec![
(
mean_frozen.clone(),
mean_design.clone(),
mean_term_indices.clone(),
),
(
noise_frozen.clone(),
noise_design.clone(),
noise_term_indices.clone(),
),
],
joint_setup.rho_dim(),
joint_setup.log_kappa_dims_per_term(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "n-block cache", e));
cache.ensure_theta(&theta0).unwrap_or_else(|e| panic!("{} failed: {:?}", "initial theta", e));
assert!(cache.memoized_cost(&theta0).is_none());
assert!(cache.memoized_eval(&theta0).is_none());
let eval = (
2.25,
Array1::<f64>::ones(theta0.len()),
gam_problem::HessianResult::Analytic(Array2::<f64>::eye(theta0.len())),
);
cache.store_eval(eval.clone());
let cached_eval = cache.memoized_eval(&theta0).unwrap_or_else(|| panic!("{} failed", "cached eval"));
assert!((cached_eval.0 - eval.0).abs() <= 1e-12);
assert_eq!(cached_eval.1, eval.1);
assert_eq!(
cached_eval
.2
.materialize_dense()
.unwrap_or_else(|e| panic!("{} failed: {:?}", "materialize cached hessian", e)),
eval.2
.materialize_dense()
.expect("materialize eval hessian"),
);
let mut theta1 = theta0.clone();
theta1[joint_setup.rho_dim()] += 0.25;
cache.ensure_theta(&theta1).unwrap_or_else(|e| panic!("{} failed: {:?}", "updated theta", e));
assert!(cache.memoized_cost(&theta1).is_none());
assert!(cache.memoized_eval(&theta1).is_none());
let log_kappa = SpatialLogKappaCoords::from_theta_tail_with_dims(
&theta1,
joint_setup.rho_dim(),
joint_setup.log_kappa_dims_per_term(),
);
let mean_terms = spatial_length_scale_term_indices(&mean_frozen);
let noise_terms = spatial_length_scale_term_indices(&noise_frozen);
let (mean_lk, noise_lk) = log_kappa.split_at(mean_terms.len());
let mean_updated = mean_lk
.apply_tospec(&mean_frozen, &mean_terms)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "mean updated spec", e));
let noise_updated = noise_lk
.apply_tospec(&noise_frozen, &noise_terms)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "noise updated spec", e));
let mean_rebuilt =
build_term_collection_design(data.view(), &mean_updated).unwrap_or_else(|e| panic!("{} failed: {:?}", "mean rebuilt", e));
let noise_rebuilt =
build_term_collection_design(data.view(), &noise_updated).unwrap_or_else(|e| panic!("{} failed: {:?}", "noise rebuilt", e));
let cache_designs = cache.designs();
assert_term_collection_designs_match(cache_designs[0], &mean_rebuilt, "mean cache");
assert_term_collection_designs_match(cache_designs[1], &noise_rebuilt, "noise cache");
}
#[test]
fn single_block_exact_joint_design_cache_clears_memo_on_theta_change() {
let n = 22usize;
let mut data = Array2::<f64>::zeros((n, 2));
for i in 0..n {
let x0 = i as f64 / (n as f64 - 1.0);
let x1 = (0.23 * i as f64).cos();
data[[i, 0]] = x0;
data[[i, 1]] = x1;
}
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "duchon_hybrid".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0, 1],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
length_scale: Some(0.9),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let frozen = freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze spec", e));
let spatial_terms = spatial_length_scale_term_indices(&frozen);
let rho_dim = design.penalties.len();
let dims_per_term = vec![1];
let mut theta0 = Array1::<f64>::zeros(rho_dim + 1);
theta0[rho_dim] = -get_spatial_length_scale(&frozen, spatial_terms[0])
.unwrap_or_else(|| panic!("{} failed", "length scale"))
.ln();
let mut cache = SingleBlockExactJointDesignCache::new(
data.view(),
frozen.clone(),
design.clone(),
spatial_terms.clone(),
rho_dim,
dims_per_term.clone(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "single-block cache", e));
cache.ensure_theta(&theta0).unwrap_or_else(|e| panic!("{} failed: {:?}", "initial theta", e));
assert!(cache.memoized_cost(&theta0).is_none());
assert!(cache.memoized_eval(&theta0).is_none());
let eval = (
0.5,
Array1::<f64>::ones(theta0.len()),
gam_problem::HessianResult::Analytic(Array2::<f64>::eye(theta0.len())),
);
cache.store_eval_at(&theta0, eval.clone());
let cached_eval = cache.memoized_eval(&theta0).unwrap_or_else(|| panic!("{} failed", "cached eval"));
assert!((cached_eval.0 - eval.0).abs() <= 1e-12);
assert_eq!(cached_eval.1, eval.1);
assert_eq!(
cached_eval
.2
.materialize_dense()
.unwrap_or_else(|e| panic!("{} failed: {:?}", "materialize cached hessian", e)),
eval.2
.materialize_dense()
.expect("materialize eval hessian"),
);
let mut theta1 = theta0.clone();
theta1[rho_dim] += 0.35;
cache.ensure_theta(&theta1).unwrap_or_else(|e| panic!("{} failed: {:?}", "updated theta", e));
assert!(cache.memoized_cost(&theta1).is_none());
assert!(cache.memoized_eval(&theta1).is_none());
let updated_log_kappa =
SpatialLogKappaCoords::from_theta_tail_with_dims(&theta1, rho_dim, dims_per_term);
let updated_spec = updated_log_kappa
.apply_tospec(&frozen, &spatial_terms)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "updated spec", e));
let rebuilt = build_term_collection_design(data.view(), &updated_spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "rebuilt design", e));
assert_term_collection_designs_match(cache.design(), &rebuilt, "single-block cache");
}
#[test]
fn single_block_latent_coord_design_cache_invalidates_memo_on_outer_iter_advance() {
use gam_solve::estimate::reml::outer_eval::{
current_outer_iter, record_current_outer_iter_for_ift,
};
use gam_terms::latent::{LatentCoordValues, LatentIdMode, LatentManifold};
let n = 16usize;
let latent_dim = 1usize;
let mut data = Array2::<f64>::zeros((n, 1));
for i in 0..n {
data[[i, 0]] = i as f64 / (n as f64 - 1.0);
}
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "s(x0)".to_string(),
basis: SmoothBasisSpec::BSpline1D {
feature_col: 0,
spec: BSplineBasisSpec {
degree: 3,
penalty_order: 2,
knotspec: BSplineKnotSpec::Generate {
data_range: (0.0, 1.0),
num_internal_knots: 3,
},
double_penalty: false,
identifiability: BSplineIdentifiability::WeightedSumToZero { weights: None },
boundary_conditions: BSplineBoundaryConditions::default(),
boundary: OneDimensionalBoundary::Open,
},
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let rho_dim = design.penalties.len();
let flat = Array1::<f64>::from_iter((0..n).map(|i| i as f64 / (n as f64 - 1.0)));
let latent_values = std::sync::Arc::new(LatentCoordValues::from_flat_with_manifold(
flat,
n,
latent_dim,
LatentIdMode::None,
LatentManifold::Euclidean,
));
let latent = StandardLatentCoordConfig {
values: latent_values,
term_index: gam_problem::SmoothTermIdx::new(0),
feature_cols: vec![0],
manifold: LatentManifold::Euclidean,
manifold_auto: false,
retraction_registry: gam_problem::LatentRetractionRegistry::default(),
analytic_penalties: None,
};
let mut cache = SingleBlockLatentCoordDesignCache::new(
data.clone(),
spec.clone(),
design.clone(),
&latent,
rho_dim,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "latent-coord cache", e));
let theta = Array1::<f64>::zeros(rho_dim + n * latent_dim);
cache.current_theta = Some(theta.clone());
let initial_iter = current_outer_iter();
let eval = (
1.25_f64,
Array1::<f64>::from_elem(theta.len(), 0.5),
gam_problem::HessianResult::Analytic(Array2::<f64>::eye(theta.len())),
);
cache.store_eval(eval.clone());
let hit = cache
.memoized_eval(&theta)
.unwrap_or_else(|| panic!("{} failed", "memo should hit at the same outer iter"));
assert!((hit.0 - eval.0).abs() <= 1e-12);
assert_eq!(cache.memoized_cost(&theta), Some(eval.0));
record_current_outer_iter_for_ift(initial_iter.wrapping_add(1));
assert!(
cache.memoized_eval(&theta).is_none(),
"memoized_eval returned a stale cached eval after current_outer_iter advanced"
);
assert!(
cache.memoized_cost(&theta).is_none(),
"memoized_cost returned a stale cached cost after current_outer_iter advanced"
);
record_current_outer_iter_for_ift(initial_iter);
}
#[test]
fn external_joint_evaluator_reuse_matches_fresh_state_after_theta_update() {
let n = 26usize;
let mut data = Array2::<f64>::zeros((n, 2));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let x0 = i as f64 / (n as f64 - 1.0);
let x1 = (0.21 * i as f64).sin();
data[[i, 0]] = x0;
data[[i, 1]] = x1;
y[i] = (2.0 * std::f64::consts::PI * x0).sin() + 0.35 * x1;
}
let spec = TermCollectionSpec {
linear_terms: vec![LinearTermSpec {
name: "x0".to_string(),
feature_col: 0,
feature_cols: vec![0],
categorical_levels: vec![],
double_penalty: false,
coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
coefficient_min: None,
coefficient_max: None,
}],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "matern".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
length_scale: 0.85,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let weights = Array1::ones(n);
let offset = Array1::zeros(n);
let fit_opts = FitOptions {
compute_inference: false,
max_iter: 40,
tol: 1e-7,
penalty_shrinkage_floor: None,
..FitOptions::default()
};
let design = build_term_collection_design(data.view(), &spec).unwrap_or_else(|e| panic!("{} failed: {:?}", "design", e));
let frozen = freeze_term_collection_from_design(&spec, &design).unwrap_or_else(|e| panic!("{} failed: {:?}", "freeze", e));
let spatial_terms = spatial_length_scale_term_indices(&frozen);
let dims_per_term = spatial_dims_per_term(&frozen, &spatial_terms);
let rho_dim = design.penalties.len();
let mut theta0 = Array1::<f64>::zeros(rho_dim + dims_per_term.iter().sum::<usize>());
for j in 0..rho_dim {
theta0[j] = 0.2 - 0.1 * j as f64;
}
theta0[rho_dim] = -get_spatial_length_scale(&frozen, spatial_terms[0])
.unwrap_or_else(|| panic!("{} failed", "length scale"))
.ln();
let mut theta1 = theta0.clone();
theta1[rho_dim] += 0.3;
let external_opts =
external_opts_for_design(&LikelihoodSpec::gaussian_identity(), &design, &fit_opts);
let mut cache = SingleBlockExactJointDesignCache::new(
data.view(),
frozen,
design.clone(),
spatial_terms,
rho_dim,
dims_per_term,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "single-block cache", e));
let mut reused = gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&design.design,
offset.view(),
&design.penalties,
&external_opts,
"reused evaluator",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "reused evaluator", e));
let compare_eval =
|theta: &Array1<f64>,
cache: &mut SingleBlockExactJointDesignCache<'_>,
reused: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>| {
cache.ensure_theta(theta).unwrap_or_else(|e| panic!("{} failed: {:?}", "theta applied", e));
let build_hyper_dirs = || {
try_build_spatial_log_kappa_hyper_dirs(
data.view(),
cache.spec(),
cache.design(),
&cache.spatial_terms,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "hyper dirs build", e))
.expect("hyper dirs present")
};
let reused_eval = evaluate_joint_reml_outer_eval_at_theta(
reused,
cache.design(),
theta,
rho_dim,
build_hyper_dirs(),
None,
gam_model_api::OuterEvalOrder::ValueGradientHessian,
None,
)
.expect("reused eval");
let fresh_opts = external_opts_for_design(
&LikelihoodSpec::gaussian_identity(),
cache.design(),
&fit_opts,
);
let mut fresh = gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&cache.design().design,
offset.view(),
&cache.design().penalties,
&fresh_opts,
"fresh evaluator",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "fresh evaluator", e));
let fresh_eval = evaluate_joint_reml_outer_eval_at_theta(
&mut fresh,
cache.design(),
theta,
rho_dim,
build_hyper_dirs(),
None,
gam_model_api::OuterEvalOrder::ValueGradientHessian,
None,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "fresh eval", e));
let cost_diff = (reused_eval.0 - fresh_eval.0).abs();
assert!(cost_diff <= 1e-10, "cost mismatch: {cost_diff}");
let grad_diff = reused_eval
.1
.iter()
.zip(fresh_eval.1.iter())
.map(|(left, right)| (left - right).abs())
.fold(0.0_f64, f64::max);
assert!(grad_diff <= 1e-9, "gradient mismatch: {grad_diff}");
let reused_hess = reused_eval
.2
.materialize_dense()
.unwrap_or_else(|e| panic!("{} failed: {:?}", "reused hessian materializes", e))
.expect("reused hessian present");
let fresh_hess = fresh_eval
.2
.materialize_dense()
.unwrap_or_else(|e| panic!("{} failed: {:?}", "fresh hessian materializes", e))
.expect("fresh hessian present");
let hess_diff = max_abs_diff_matrix(&reused_hess, &fresh_hess);
assert!(hess_diff <= 1e-9, "hessian mismatch: {hess_diff}");
let reused_efs = evaluate_joint_reml_efs_at_theta(
reused,
cache.design(),
theta,
rho_dim,
build_hyper_dirs(),
None,
None,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "reused EFS eval", e));
let mut fresh_efs_eval = gam_solve::estimate::ExternalJointHyperEvaluator::new(
y.view(),
weights.view(),
&cache.design().design,
offset.view(),
&cache.design().penalties,
&fresh_opts,
"fresh EFS evaluator",
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "fresh EFS evaluator", e));
let fresh_efs = evaluate_joint_reml_efs_at_theta(
&mut fresh_efs_eval,
cache.design(),
theta,
rho_dim,
build_hyper_dirs(),
None,
None,
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "fresh EFS eval", e));
let efs_cost_diff = (reused_efs.cost - fresh_efs.cost).abs();
assert!(efs_cost_diff <= 1e-10, "EFS cost mismatch: {efs_cost_diff}");
assert_eq!(reused_efs.steps.len(), fresh_efs.steps.len());
let efs_step_diff = reused_efs
.steps
.iter()
.zip(fresh_efs.steps.iter())
.map(|(left, right)| (left - right).abs())
.fold(0.0_f64, f64::max);
assert!(efs_step_diff <= 1e-9, "EFS step mismatch: {efs_step_diff}");
};
compare_eval(&theta0, &mut cache, &mut reused);
compare_eval(&theta1, &mut cache, &mut reused);
}
#[test]
fn exact_matern_log_kappa_derivative_uses_feature_columns_only() {
let n = 24usize;
let p = 17usize;
let mut data = Array2::<f64>::zeros((n, p));
for i in 0..n {
let x = i as f64 / (n as f64 - 1.0);
data[[i, 0]] = x;
for j in 1..p {
data[[i, j]] = ((i + j) as f64 * 0.13).sin();
}
}
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "matern".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
length_scale: 0.4,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "baseline Matérn design should build", e));
let frozenspec = freeze_term_collection_from_design(&spec, &design)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "freezing Matérn centers from design should succeed", e));
match &frozenspec.smooth_terms[0].basis {
SmoothBasisSpec::Matern { spec, .. } => match &spec.center_strategy {
CenterStrategy::UserProvided(centers) => {
assert_eq!(centers.ncols(), 1, "frozen centers should stay term-local");
}
_ => panic!("expected frozen user-provided centers"),
},
_ => panic!("expected Matérn term"),
}
let derivative =
try_build_spatial_term_log_kappa_derivative(data.view(), &frozenspec, &design, 0);
assert!(
derivative.is_ok(),
"exact Matérn log-kappa derivative should use only feature_cols; got {derivative:?}"
);
assert!(
derivative
.unwrap_or_else(|e| panic!("{} failed: {:?}", "derivative call should succeed", e))
.is_some(),
"Matérn term should expose an exact derivative"
);
}
#[test]
fn exact_thin_plate_log_kappa_derivative_uses_feature_columns_only() {
let n = 28usize;
let p = 15usize;
let mut data = Array2::<f64>::zeros((n, p));
for i in 0..n {
let x0 = i as f64 / (n as f64 - 1.0);
let x1 = (0.17 * i as f64).sin();
data[[i, 0]] = x0;
data[[i, 1]] = x1;
for j in 2..p {
data[[i, j]] = ((i + 3 * j) as f64 * 0.07).cos();
}
}
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "thinplate".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 7 },
length_scale: 0.7,
double_penalty: true,
identifiability: SpatialIdentifiability::default(),
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "baseline ThinPlate design should build", e));
let frozenspec = freeze_term_collection_from_design(&spec, &design)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "freezing ThinPlate centers from design should succeed", e));
match &frozenspec.smooth_terms[0].basis {
SmoothBasisSpec::ThinPlate { spec, .. } => match &spec.center_strategy {
CenterStrategy::UserProvided(centers) => {
assert_eq!(centers.ncols(), 2, "frozen centers should stay term-local");
}
_ => panic!("expected frozen user-provided centers"),
},
_ => panic!("expected ThinPlate term"),
}
let smooth_term = &design.smooth.terms[0];
let termspec = &frozenspec.smooth_terms[0];
let BasisPsiDerivativeResult {
design_derivative: local_x_psi,
penalties_derivative: local_s_psi,
..
} = match &termspec.basis {
SmoothBasisSpec::ThinPlate {
feature_cols, spec, ..
} => {
let x =
select_columns(data.view(), feature_cols).unwrap_or_else(|e| panic!("{} failed: {:?}", "select ThinPlate feature cols", e));
gam_terms::basis::build_thin_plate_basis_log_kappa_derivative(x.view(), spec)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "direct ThinPlate derivative should build", e))
}
_ => panic!("expected ThinPlate term"),
};
let BasisPsiSecondDerivativeResult {
designsecond_derivative: local_x_psi_psi,
penaltiessecond_derivative: local_s_psi_psi,
..
} = match &termspec.basis {
SmoothBasisSpec::ThinPlate {
feature_cols, spec, ..
} => {
let x =
select_columns(data.view(), feature_cols).expect("select ThinPlate feature cols");
gam_terms::basis::build_thin_plate_basis_log_kappasecond_derivative(x.view(), spec)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "direct ThinPlate second derivative should build", e))
}
_ => panic!("expected ThinPlate term"),
};
assert_eq!(local_x_psi.ncols(), smooth_term.coeff_range.len());
assert_eq!(local_x_psi_psi.ncols(), smooth_term.coeff_range.len());
assert!(!local_s_psi.is_empty());
assert_eq!(local_s_psi.len(), local_s_psi_psi.len());
assert!(local_s_psi.iter().all(|s| {
s.nrows() == smooth_term.coeff_range.len() && s.ncols() == smooth_term.coeff_range.len()
}));
assert!(local_s_psi_psi.iter().all(|s| {
s.nrows() == smooth_term.coeff_range.len() && s.ncols() == smooth_term.coeff_range.len()
}));
let derivative =
try_build_spatial_term_log_kappa_derivative(data.view(), &frozenspec, &design, 0);
assert!(
derivative.is_ok(),
"exact ThinPlate log-kappa derivative should use only feature_cols; got {derivative:?}"
);
let derivative = derivative.unwrap_or_else(|e| panic!("{} failed: {:?}", "derivative call should succeed", e));
assert!(
derivative.is_some(),
"ThinPlate term should expose an exact derivative"
);
}
#[test]
fn exact_duchon_log_kappa_derivative_uses_feature_columns_only() {
let n = 28usize;
let p = 15usize;
let mut data = Array2::<f64>::zeros((n, p));
for i in 0..n {
let x0 = i as f64 / (n as f64 - 1.0);
let x1 = (0.21 * i as f64).cos();
data[[i, 0]] = x0;
data[[i, 1]] = x1;
for j in 2..p {
data[[i, j]] = ((i + 2 * j) as f64 * 0.09).sin();
}
}
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "duchon".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0, 1],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 7 },
length_scale: Some(0.7),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let design = build_term_collection_design(data.view(), &spec)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "baseline Duchon design should build", e));
let frozenspec = freeze_term_collection_from_design(&spec, &design)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "freezing Duchon centers from design should succeed", e));
match &frozenspec.smooth_terms[0].basis {
SmoothBasisSpec::Duchon { spec, .. } => match &spec.center_strategy {
CenterStrategy::UserProvided(centers) => {
assert_eq!(centers.ncols(), 2, "frozen centers should stay term-local");
}
_ => panic!("expected frozen user-provided centers"),
},
_ => panic!("expected Duchon term"),
}
let smooth_term = &design.smooth.terms[0];
let termspec = &frozenspec.smooth_terms[0];
let derivative_bundle = match &termspec.basis {
SmoothBasisSpec::Duchon {
feature_cols, spec, ..
} => {
let x = select_columns(data.view(), feature_cols).unwrap_or_else(|e| panic!("{} failed: {:?}", "select Duchon feature cols", e));
gam_terms::basis::build_duchon_basis_log_kappa_derivatives(x.view(), spec)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "direct Duchon derivative bundle should build", e))
}
_ => panic!("expected Duchon term"),
};
let local_implicit = derivative_bundle.implicit_operator;
let BasisPsiDerivativeResult {
design_derivative: local_x_psi,
penalties_derivative: local_s_psi,
implicit_operator: local_implicit_psi_unused,
} = derivative_bundle.first;
let BasisPsiSecondDerivativeResult {
designsecond_derivative: local_x_psi_psi,
penaltiessecond_derivative: local_s_psi_psi,
implicit_operator: local_implicit_psi_psi_unused,
} = derivative_bundle.second;
assert!(local_implicit_psi_unused.is_none());
assert!(local_implicit_psi_psi_unused.is_none());
assert_spatial_derivative_width(
"Duchon first log-kappa",
&local_x_psi,
local_implicit.as_ref(),
smooth_term.coeff_range.len(),
);
assert_spatial_derivative_width(
"Duchon second log-kappa",
&local_x_psi_psi,
local_implicit.as_ref(),
smooth_term.coeff_range.len(),
);
assert!(!local_s_psi.is_empty());
assert_eq!(local_s_psi.len(), local_s_psi_psi.len());
assert!(local_s_psi.iter().all(|s| {
s.nrows() == smooth_term.coeff_range.len() && s.ncols() == smooth_term.coeff_range.len()
}));
assert!(local_s_psi_psi.iter().all(|s| {
s.nrows() == smooth_term.coeff_range.len() && s.ncols() == smooth_term.coeff_range.len()
}));
let derivative =
try_build_spatial_term_log_kappa_derivative(data.view(), &frozenspec, &design, 0);
assert!(
derivative.is_ok(),
"exact Duchon log-kappa derivative should use only feature_cols; got {derivative:?}"
);
let derivative = derivative.unwrap_or_else(|e| panic!("{} failed: {:?}", "derivative call should succeed", e));
assert!(
derivative.is_some(),
"Duchon term should expose an exact derivative"
);
}
#[test]
fn spatial_length_scale_optimization_monotone_improves_or_keeps_score_for_matern() {
let n = 60usize;
let d = 2usize;
let mut data = Array2::<f64>::zeros((n, d));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let x0 = i as f64 / (n as f64 - 1.0);
let x1 = (i as f64 * 0.17).sin();
data[[i, 0]] = x0;
data[[i, 1]] = x1;
y[i] = (3.0 * x0).cos() + 0.35 * x1;
}
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "matern".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 12 },
length_scale: 12.0,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let fit_opts = FitOptions {
max_iter: 40,
penalty_shrinkage_floor: None,
..FitOptions::default()
};
let weights = Array1::ones(n);
let offset = Array1::zeros(n);
assert_matern_spatial_length_scale_optimization_monotone(
data.view(),
&y,
&weights,
&offset,
&spec,
&fit_opts,
);
}
#[test]
fn spatial_length_scale_optimization_runs_binomial_logit_matern_with_exact_laml_derivatives() {
let n = 80usize;
let d = 2usize;
let mut data = Array2::<f64>::zeros((n, d));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let x0 = i as f64 / (n as f64 - 1.0);
let x1 = (i as f64 * 0.19).cos();
data[[i, 0]] = x0;
data[[i, 1]] = x1;
let eta = -0.15 + 0.45 * x0 - 0.25 * x1 + 0.10 * (6.0 * x0).sin();
let mu = 1.0 / (1.0 + (-eta).exp());
let u = (((i * 37 + 17) % 101) as f64 + 0.5) / 101.0;
y[i] = if u < mu { 1.0 } else { 0.0 };
}
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "matern".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 10 },
length_scale: 1.8,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::CenterSumToZero,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let fit_opts = FitOptions {
max_iter: 60,
..FitOptions::default()
};
let weights = Array1::ones(n);
let offset = Array1::zeros(n);
fit_term_collectionwith_spatial_length_scale_optimization(
data.view(),
y,
weights,
offset,
&spec,
LikelihoodSpec::binomial_logit(),
&fit_opts,
&SpatialLengthScaleOptimizationOptions {
max_outer_iter: 2,
rel_tol: 1e-5,
pilot_subsample_threshold: 0,
outer_wall_clock_budget_secs: None,
..SpatialLengthScaleOptimizationOptions::default()
},
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "standard binomial-logit spatial kappa optimization should use exact non-TK LAML derivatives", e));
}
#[test]
fn spatial_kappa_result_requires_exact_availability() {
let err = require_successful_spatial_optimization_result::<()>(0.0, Ok(None))
.expect_err("missing exact spatial result must be surfaced");
let msg = err.to_string();
assert!(msg.contains("unavailable"), "unexpected error: {msg}");
}
#[test]
fn spatial_kappa_result_rejects_worse_exact_score() {
let err = require_successful_spatial_optimization_result(1.0, Ok(Some(((), 1.5))))
.expect_err("worse exact spatial result must be rejected");
let msg = err.to_string();
assert!(
msg.contains("made REML score worse"),
"unexpected error: {msg}"
);
assert!(msg.contains("1.000000e0"), "unexpected error: {msg}");
assert!(msg.contains("1.500000e0"), "unexpected error: {msg}");
}
#[test]
fn spatial_kappa_result_surfaces_optimizer_failure() {
let err = require_successful_spatial_optimization_result::<()>(
0.0,
Err(EstimationError::InvalidInput("boom".to_string())),
)
.expect_err("exact spatial optimizer failure must be surfaced");
let msg = err.to_string();
assert!(
msg.contains("spatial kappa optimization failed"),
"unexpected error: {msg}"
);
assert!(msg.contains("boom"), "unexpected error: {msg}");
}
#[test]
fn duchon_terms_participate_in_kappa_optimization() {
let data = array![
[0.0, 0.1, 0.2],
[0.2, 0.0, 0.4],
[0.4, 0.3, 0.1],
[0.6, 0.5, 0.7],
[0.8, 0.7, 0.3],
[1.0, 0.9, 0.8],
];
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "duchon".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0, 1],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: Some(0.9),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
assert_eq!(spatial_length_scale_term_indices(&spec), vec![0]);
let fit_opts = FitOptions {
max_iter: 40,
penalty_shrinkage_floor: None,
..FitOptions::default()
};
let y = Array1::linspace(0.0, 1.0, data.nrows());
let weights = Array1::ones(data.nrows());
let offset = Array1::zeros(data.nrows());
let design = build_term_collection_design(data.view(), &spec)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "baseline Duchon design should build", e));
let frozenspec = freeze_term_collection_from_design(&spec, &design)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "freezing Duchon centers from design should succeed", e));
let derivative =
try_build_spatial_term_log_kappa_derivative(data.view(), &frozenspec, &design, 0);
assert!(
derivative
.unwrap_or_else(|e| panic!("{} failed: {:?}", "Duchon exact derivative call should succeed", e))
.is_some(),
"Duchon term should expose an exact derivative"
);
let optimized = fit_term_collectionwith_spatial_length_scale_optimization(
data.view(),
y,
weights,
offset,
&spec,
LikelihoodSpec::gaussian_identity(),
&fit_opts,
&SpatialLengthScaleOptimizationOptions::default(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "Duchon fit should use exact κ optimization", e));
let optimized_ls = match &optimized.resolvedspec.smooth_terms[0].basis {
SmoothBasisSpec::Duchon { spec, .. } => spec.length_scale,
_ => panic!("expected Duchon term"),
};
assert!(optimized_ls.is_some());
match &optimized.resolvedspec.smooth_terms[0].basis {
SmoothBasisSpec::Duchon { spec, .. } => {
assert!(matches!(
spec.center_strategy,
CenterStrategy::UserProvided(_)
));
assert!(matches!(
spec.identifiability,
SpatialIdentifiability::FrozenTransform { .. }
));
}
_ => panic!("expected Duchon term"),
}
}
#[test]
fn pure_duchon_scale_dimensions_seed_geometry_but_enroll_no_hyper_axis() {
let mut spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "pure_duchon".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0, 1],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 4 },
length_scale: None,
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: None,
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
gam_terms::term_builder::enable_scale_dimensions(&mut spec);
assert!(
spatial_length_scale_term_indices(&spec).is_empty(),
"pure Duchon must enroll no outer hyper axis even with scale_dims on"
);
match &spec.smooth_terms[0].basis {
SmoothBasisSpec::Duchon { spec, .. } => {
assert_eq!(spec.length_scale, None);
assert_eq!(spec.aniso_log_scales.as_deref(), Some(&[0.0, 0.0][..]));
}
_ => panic!("expected Duchon term"),
}
}
#[test]
fn thin_plate_terms_anchor_length_scale_and_enroll_no_kappa_axis() {
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "thin_plate".to_string(),
basis: SmoothBasisSpec::ThinPlate {
feature_cols: vec![0, 1],
spec: ThinPlateBasisSpec {
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
length_scale: 0.75,
double_penalty: false,
identifiability: SpatialIdentifiability::default(),
radial_reparam: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
assert!(
spatial_length_scale_term_indices(&spec).is_empty(),
"penalized thin-plate regression splines must not contribute a redundant isotropic kappa axis"
);
assert!(
all_spatial_terms_kappa_fixed(&spec),
"with no TPS kappa axis, all spatial terms are effectively fixed-geometry"
);
}
#[test]
fn pure_duchon_from_length_scales_aniso_is_isotropic_single_psi() {
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "pure_duchon".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0, 1, 2],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::UserProvided(array![
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0],
]),
length_scale: None,
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::None,
aniso_log_scales: Some(vec![0.7, 0.2, 0.1]),
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let opts = SpatialLengthScaleOptimizationOptions::default();
let coords = SpatialLogKappaCoords::from_length_scales_aniso(&spec, &[0], &opts);
assert_eq!(coords.dims_per_term(), &[1]);
assert_eq!(coords.as_array().len(), 1);
let expected_psi = -opts.min_length_scale.ln();
assert!((coords.as_array()[0] - expected_psi).abs() <= 1e-12);
}
#[test]
fn explicit_duchon_aniso_length_scale_is_locked_kappa() {
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "duchon_fixed_geometry".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0, 1, 2],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::UserProvided(array![
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0],
]),
length_scale: Some(1.0),
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::None,
aniso_log_scales: Some(vec![0.7, 0.2, 0.1]),
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
assert!(
spatial_term_has_locked_kappa(&spec, 0),
"Duchon anisotropy is fixed geometry and must not force ψ optimization"
);
assert!(
all_spatial_terms_kappa_fixed(&spec),
"a Duchon term with explicit length_scale and fixed anisotropy has no REML κ/ψ axis"
);
}
#[test]
fn from_length_scales_aniso_keeps_nonaniso_spatial_terms_scalar() {
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![
SmoothTermSpec {
name: "matern_aniso".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::UserProvided(array![
[0.0, 0.0],
[1.0, 0.0],
[0.0, 1.0],
]),
length_scale: 0.5,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: false,
identifiability: MaternIdentifiability::None,
aniso_log_scales: Some(vec![0.3, -0.3]),
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
SmoothTermSpec {
name: "matern_iso".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::UserProvided(array![
[0.0, 0.0],
[1.0, 0.0],
[0.0, 1.0],
]),
length_scale: 0.25,
nu: MaternNu::ThreeHalves,
include_intercept: false,
double_penalty: false,
identifiability: MaternIdentifiability::None,
aniso_log_scales: None,
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
},
],
};
let term_indices = [0usize, 1usize];
let coords = SpatialLogKappaCoords::from_length_scales_aniso(
&spec,
&term_indices,
&SpatialLengthScaleOptimizationOptions::default(),
);
assert_eq!(spatial_dims_per_term(&spec, &term_indices), vec![2, 1]);
assert_eq!(coords.dims_per_term(), &[2, 1]);
let expected = [-0.5_f64.ln() + 0.3, -0.5_f64.ln() - 0.3, -0.25_f64.ln()];
for (got, want) in coords.as_array().iter().zip(expected.iter()) {
assert!((got - want).abs() <= 1e-12);
}
}
#[test]
fn aniso_bounds_clamp_preserves_in_range_global_length_scale_and_eta() {
let data = array![[0.0, 0.0], [1.0, 0.2], [0.1, 1.0], [1.1, 1.2]];
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "matern_aniso".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::UserProvided(array![
[0.0, 0.0],
[1.0, 0.0],
[0.0, 1.0],
[1.0, 1.0],
]),
length_scale: 1.0,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::None,
aniso_log_scales: Some(vec![3.0, -3.0]),
nullspace_shrinkage_survived: None,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let options = SpatialLengthScaleOptimizationOptions {
max_outer_iter: 1,
rel_tol: 1e-6,
min_length_scale: (-2.0_f64).exp(),
max_length_scale: 1.0_f64.exp(),
pilot_subsample_threshold: 0,
outer_wall_clock_budget_secs: None,
..SpatialLengthScaleOptimizationOptions::default()
};
let spatial_terms = vec![0];
let dims_per_term = spatial_dims_per_term(&spec, &spatial_terms);
let seed = SpatialLogKappaCoords::from_length_scales_aniso(&spec, &spatial_terms, &options);
let lower = SpatialLogKappaCoords::lower_bounds_aniso_from_data(
data.view(),
&spec,
&spatial_terms,
&dims_per_term,
&options,
);
let upper = SpatialLogKappaCoords::upper_bounds_aniso_from_data(
data.view(),
&spec,
&spatial_terms,
&dims_per_term,
&options,
);
let projected = seed.clone().clamp_to_bounds(&lower, &upper);
assert_eq!(projected.as_array(), seed.as_array());
let updated = projected
.apply_tospec(&spec, &spatial_terms)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "aniso projection should decode", e));
match &updated.smooth_terms[0].basis {
SmoothBasisSpec::Matern { spec, .. } => {
assert!((spec.length_scale - 1.0).abs() <= 1e-12);
let eta = spec
.aniso_log_scales
.as_ref()
.unwrap_or_else(|| panic!("{} failed", "anisotropy should be preserved"));
assert!((eta[0] - 3.0).abs() <= 1e-12);
assert!((eta[1] + 3.0).abs() <= 1e-12);
}
_ => panic!("expected Matérn term"),
}
}
#[test]
fn pure_duchon_aniso_fit_optimizes_without_introducing_hybrid_scale() {
let data = array![
[0.0, 0.1, 0.2],
[0.2, 0.0, 0.4],
[0.4, 0.3, 0.1],
[0.6, 0.5, 0.7],
[0.8, 0.7, 0.3],
[1.0, 0.9, 0.8],
];
let spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "pure_duchon".to_string(),
basis: SmoothBasisSpec::Duchon {
feature_cols: vec![0, 1, 2],
spec: DuchonBasisSpec {
radial_reparam: None,
periodic: None,
center_strategy: CenterStrategy::FarthestPoint { num_centers: 5 },
length_scale: None,
power: 1.0,
nullspace_order: DuchonNullspaceOrder::Linear,
identifiability: SpatialIdentifiability::default(),
aniso_log_scales: Some(vec![0.0, 0.0, 0.0]),
operator_penalties: DuchonOperatorPenaltySpec::default(),
boundary: OneDimensionalBoundary::Open,
},
input_scales: None,
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let fit_opts = FitOptions {
max_iter: 40,
penalty_shrinkage_floor: None,
..FitOptions::default()
};
let optimized = fit_term_collectionwith_spatial_length_scale_optimization(
data.view(),
Array1::linspace(0.0, 1.0, data.nrows()),
Array1::ones(data.nrows()),
Array1::zeros(data.nrows()),
&spec,
LikelihoodSpec::gaussian_identity(),
&fit_opts,
&SpatialLengthScaleOptimizationOptions::default(),
)
.unwrap_or_else(|e| panic!("{} failed: {:?}", "pure Duchon anisotropic fit should optimize", e));
match &optimized.resolvedspec.smooth_terms[0].basis {
SmoothBasisSpec::Duchon { spec, .. } => {
assert_eq!(spec.length_scale, None);
assert!(
spec.aniso_log_scales.is_some(),
"pure Duchon anisotropy should remain enabled"
);
}
_ => panic!("expected Duchon term"),
}
}
#[test]
fn spatial_anisotropy_pilot_initializer_seeds_geometry_without_fit() {
let data = Array2::from_shape_fn((32, 2), |(i, j)| {
if j == 0 {
i as f64 / 31.0
} else {
((i % 8) as f64) * 0.03
}
});
let mut spec = TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![SmoothTermSpec {
name: "pc_matern".to_string(),
basis: SmoothBasisSpec::Matern {
feature_cols: vec![0, 1],
spec: MaternBasisSpec {
periodic: None,
center_strategy: CenterStrategy::UserProvided(array![
[0.0, 0.0],
[1.0, 0.0],
[0.0, 0.05],
[1.0, 0.05],
]),
length_scale: 1.0,
nu: MaternNu::FiveHalves,
include_intercept: false,
double_penalty: true,
identifiability: MaternIdentifiability::None,
aniso_log_scales: Some(vec![0.0, 0.0]),
nullspace_shrinkage_survived: None,
},
input_scales: Some(vec![1.0, 1.0]),
},
shape: ShapeConstraint::None,
joint_null_rotation: None,
}],
};
let spatial_terms = spatial_length_scale_term_indices(&spec);
let updated = apply_spatial_anisotropy_pilot_initializer(
data.view(),
&mut spec,
&spatial_terms,
8,
&SpatialLengthScaleOptimizationOptions::default(),
);
assert_eq!(updated, 1);
match &spec.smooth_terms[0].basis {
SmoothBasisSpec::Matern { spec, .. } => {
let eta = spec
.aniso_log_scales
.as_ref()
.unwrap_or_else(|| panic!("{} failed", "pilot initializer should preserve anisotropy"));
assert_eq!(eta.len(), 2);
assert!((eta[0] + eta[1]).abs() <= 1e-12);
assert!(
eta.iter().any(|value| value.abs() > 1e-6),
"pilot geometry should seed nonzero axis contrast"
);
assert!(spec.length_scale.is_finite() && spec.length_scale > 0.0);
}
_ => panic!("expected Matern term"),
}
}
}