1use std::collections::{BTreeMap, BTreeSet, HashMap};
8use std::path::PathBuf;
9
10use ndarray::{Array2, ArrayView1};
11
12use crate::basis::{
13 BSplineBasisSpec, BSplineBoundaryConditions, BSplineEndpointBoundaryCondition,
14 BSplineIdentifiability, BSplineKnotSpec, CenterCountRequest, CenterStrategy,
15 ConstantCurvatureBasisSpec, ConstantCurvatureIdentifiability, DuchonBasisSpec,
16 DuchonNullspaceOrder, DuchonOperatorPenaltySpec, DuchonSpectralBasis, MaternBasisSpec,
17 MaternIdentifiability, MaternLengthScale, MaternNu, MeasureJetBasisSpec,
18 MeasureJetIdentifiability, OneDimensionalBoundary, SpatialIdentifiability, SphereMethod,
19 SphereWahbaKernel, SphericalSplineBasisSpec, SphericalSplineIdentifiability,
20 ThinPlateBasisSpec, auto_spatial_center_strategy, count_unique_coordinate_rows,
21 default_num_centers, default_spatial_center_strategy, default_spherical_harmonic_degree,
22 plan_spatial_basis, select_r_uniform_subsample_centers, thin_plate_penalty_order,
23};
24use crate::inference::formula_dsl::{
25 ParsedTerm, SmoothKind, option_bool, option_f64, option_f64_strict, option_usize,
26 option_usize_any, option_usize_any_strict, option_usize_strict, strip_quotes,
27};
28use crate::smooth::{
29 BySmoothKind, ByVarKind, ByVariableSpec, FactorSmoothFlavour, FactorSmoothSpec,
30 LinearCoefficientGeometry, LinearTermSpec, RandomEffectTermSpec, ShapeConstraint,
31 SmoothBasisSpec, SmoothTermSpec, TensorBSplineIdentifiability,
32 TensorBSplinePenaltyDecomposition, TensorBSplineSpec, TermCollectionSpec,
33};
34use gam_data::{ColumnKindTag, DataError, EncodedDataset as Dataset};
35use gam_problem::types::ColIdx;
36use gam_runtime::resource::ResourcePolicy;
37
38const DEFAULT_BSPLINE_DEGREE: usize = 3;
42
43const DEFAULT_PENALTY_ORDER: usize = 2;
47
48const SPHERE_TRUNCATION_LMAX_RANGE: std::ops::RangeInclusive<usize> = 5..=200;
53
54const CYCLIC_DEFAULT_BASIS_DIM: usize = 12;
60
61const FACTOR_SMOOTH_DEFAULT_BASIS_DIM: usize = 10;
67
68const DEFAULT_PCA_CHUNK_SIZE: usize = 4096;
72
73#[derive(Clone, Debug)]
83pub enum TermBuilderError {
84 MissingColumn { reason: String },
91 ColumnNotFound {
97 name: String,
98 role: Option<String>,
99 available: Vec<String>,
100 similar: Vec<String>,
101 tsv_hint: bool,
102 },
103 IncompatibleConfig { reason: String },
107 InvalidOption { reason: String },
110 UnsupportedFeature { reason: String },
114 DegenerateData { reason: String },
117 MalformedFormula { reason: String },
120}
121
122impl std::fmt::Display for TermBuilderError {
123 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
124 match self {
125 TermBuilderError::MissingColumn { reason }
126 | TermBuilderError::IncompatibleConfig { reason }
127 | TermBuilderError::InvalidOption { reason }
128 | TermBuilderError::UnsupportedFeature { reason }
129 | TermBuilderError::DegenerateData { reason }
130 | TermBuilderError::MalformedFormula { reason } => f.write_str(reason),
131 TermBuilderError::ColumnNotFound {
137 name,
138 role,
139 available,
140 similar,
141 tsv_hint,
142 } => {
143 let canonical = DataError::ColumnNotFound {
144 name: name.clone(),
145 role: role.clone(),
146 available: available.clone(),
147 similar: similar.clone(),
148 tsv_hint: *tsv_hint,
149 };
150 std::fmt::Display::fmt(&canonical, f)
151 }
152 }
153 }
154}
155
156impl From<TermBuilderError> for String {
157 fn from(err: TermBuilderError) -> String {
158 err.to_string()
159 }
160}
161
162impl From<String> for TermBuilderError {
169 fn from(reason: String) -> Self {
170 Self::IncompatibleConfig { reason }
171 }
172}
173
174impl From<DataError> for TermBuilderError {
181 fn from(err: DataError) -> Self {
182 match err {
183 DataError::ColumnNotFound {
184 name,
185 role,
186 available,
187 similar,
188 tsv_hint,
189 } => Self::ColumnNotFound {
190 name,
191 role,
192 available,
193 similar,
194 tsv_hint,
195 },
196 DataError::SchemaMismatch { reason }
197 | DataError::ParseError { reason }
198 | DataError::EncodingFailure { reason }
199 | DataError::EmptyInput { reason }
200 | DataError::InvalidValue { reason } => Self::MissingColumn { reason },
201 }
202 }
203}
204
205impl TermBuilderError {
207 #[inline]
208 fn missing_column(reason: impl Into<String>) -> Self {
209 TermBuilderError::MissingColumn {
210 reason: reason.into(),
211 }
212 }
213 #[inline]
214 fn incompatible_config(reason: impl Into<String>) -> Self {
215 TermBuilderError::IncompatibleConfig {
216 reason: reason.into(),
217 }
218 }
219 #[inline]
220 fn invalid_option(reason: impl Into<String>) -> Self {
221 TermBuilderError::InvalidOption {
222 reason: reason.into(),
223 }
224 }
225 #[inline]
226 fn unsupported_feature(reason: impl Into<String>) -> Self {
227 TermBuilderError::UnsupportedFeature {
228 reason: reason.into(),
229 }
230 }
231 #[inline]
232 fn degenerate_data(reason: impl Into<String>) -> Self {
233 TermBuilderError::DegenerateData {
234 reason: reason.into(),
235 }
236 }
237 #[inline]
238 fn malformed_formula(reason: impl Into<String>) -> Self {
239 TermBuilderError::MalformedFormula {
240 reason: reason.into(),
241 }
242 }
243}
244
245pub fn resolve_col(col_map: &HashMap<String, usize>, name: &str) -> Result<usize, DataError> {
256 col_map
257 .get(name)
258 .copied()
259 .ok_or_else(|| DataError::column_not_found(col_map, name, None))
260}
261
262pub fn resolve_role_col(
267 col_map: &HashMap<String, usize>,
268 name: &str,
269 role: &str,
270) -> Result<usize, DataError> {
271 col_map
272 .get(name)
273 .copied()
274 .ok_or_else(|| DataError::column_not_found(col_map, name, Some(role)))
275}
276
277fn encoded_levels_for_column(ds: &Dataset, col: ColIdx) -> Vec<(u64, String)> {
278 let mut seen = BTreeSet::<u64>::new();
279 for value in ds.values.column(col.get()) {
280 if value.is_finite() {
281 seen.insert(gam_data::canonical_level_bits(*value));
282 }
283 }
284 let schema_levels = ds
285 .schema
286 .columns
287 .get(col.get())
288 .map(|column| column.levels.as_slice())
289 .unwrap_or(&[]);
290 seen.into_iter()
291 .enumerate()
292 .map(|(idx, bits)| {
293 let fallback = format!("level{}", idx + 1);
294 let label = schema_levels.get(idx).cloned().unwrap_or(fallback);
295 (bits, label)
296 })
297 .collect()
298}
299
300const DEFAULT_SIZING_ROWS_OPTION: &str = "__default_sizing_rows";
312
313fn min_categorical_by_level_rows(ds: &Dataset, by_col: usize) -> Option<usize> {
317 let mut counts: BTreeMap<u64, usize> = BTreeMap::new();
318 for value in ds.values.column(by_col) {
319 if value.is_finite() {
320 *counts
321 .entry(gam_data::canonical_level_bits(*value))
322 .or_insert(0) += 1;
323 }
324 }
325 counts.values().copied().min()
326}
327
328fn inject_by_level_sizing_rows(
332 inner_options: &mut BTreeMap<String, String>,
333 ds: &Dataset,
334 by_col: usize,
335) {
336 if matches!(
337 ds.column_kinds.get(by_col).copied(),
338 Some(ColumnKindTag::Categorical)
339 ) && let Some(min_rows) = min_categorical_by_level_rows(ds, by_col)
340 {
341 inner_options.insert(DEFAULT_SIZING_ROWS_OPTION.to_string(), min_rows.to_string());
342 }
343}
344
345pub fn column_map_with_alias(
346 col_map: &HashMap<String, usize>,
347 alias: &str,
348 target_column: &str,
349) -> HashMap<String, usize> {
350 let mut aliased = col_map.clone();
351 if let Some(idx) = col_map.get(target_column).copied() {
352 aliased.entry(alias.to_string()).or_insert(idx);
353 }
354 aliased
355}
356
357pub const MARGINAL_SLOPE_Z_ALIAS: &str = "z";
360
361pub fn marginal_slope_z_alias_is_live(col_map: &HashMap<String, usize>, z_column: &str) -> bool {
368 col_map.contains_key(z_column) && !col_map.contains_key(MARGINAL_SLOPE_Z_ALIAS)
369}
370
371pub fn build_termspec(
376 terms: &[ParsedTerm],
377 ds: &Dataset,
378 col_map: &HashMap<String, usize>,
379 inference_notes: &mut Vec<String>,
380 policy: &ResourcePolicy,
381) -> Result<TermCollectionSpec, TermBuilderError> {
382 let mut linear_terms = Vec::<LinearTermSpec>::new();
383 let mut random_terms = Vec::<RandomEffectTermSpec>::new();
384 let mut smooth_terms = Vec::<SmoothTermSpec>::new();
385 let smooth_coordinate_count = terms
386 .iter()
387 .map(|term| match term {
388 ParsedTerm::Smooth { vars, .. } => vars.len(),
389 _ => 0,
390 })
391 .sum::<usize>();
392
393 for t in terms {
394 match t {
395 ParsedTerm::Linear {
396 name,
397 explicit,
398 double_penalty,
399 coefficient_min,
400 coefficient_max,
401 } => {
402 let col = resolve_col(col_map, name)?;
403 let auto_kind = ds.column_kinds.get(col).copied().ok_or_else(|| {
404 TermBuilderError::missing_column(format!(
405 "internal column-kind lookup failed for '{name}'"
406 ))
407 .to_string()
408 })?;
409 if *explicit {
410 linear_terms.push(LinearTermSpec {
411 name: name.clone(),
412 feature_col: col,
413 feature_cols: vec![col],
414 categorical_levels: vec![],
415 double_penalty: *double_penalty,
419 coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
420 coefficient_min: *coefficient_min,
421 coefficient_max: *coefficient_max,
422 frozen_function_mass: None,
423 });
424 } else {
425 match auto_kind {
426 ColumnKindTag::Continuous | ColumnKindTag::Binary => {
427 linear_terms.push(LinearTermSpec {
428 name: name.clone(),
429 feature_col: col,
430 feature_cols: vec![col],
431 categorical_levels: vec![],
432 double_penalty: *double_penalty,
435 coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
436 coefficient_min: *coefficient_min,
437 coefficient_max: *coefficient_max,
438 frozen_function_mass: None,
439 });
440 }
441 ColumnKindTag::Categorical => {
442 if coefficient_min.is_some() || coefficient_max.is_some() {
443 return Err(TermBuilderError::incompatible_config(format!(
444 "coefficient constraints are not supported for categorical auto-random-effect term '{name}'; use group({name}) or an unconstrained numeric term"
445 )));
446 }
447 random_terms.push(RandomEffectTermSpec {
448 name: name.clone(),
449 feature_col: col,
450 drop_first_level: false,
451 penalized: true,
452 frozen_levels: None,
453 lenient_unseen: false,
459 });
460 }
461 }
462 }
463 }
464 ParsedTerm::BoundedLinear {
465 name,
466 min,
467 max,
468 prior,
469 double_penalty,
470 } => {
471 let col = resolve_col(col_map, name)?;
472 let auto_kind = ds.column_kinds.get(col).copied().ok_or_else(|| {
473 TermBuilderError::missing_column(format!(
474 "internal column-kind lookup failed for '{name}'"
475 ))
476 .to_string()
477 })?;
478 if !matches!(auto_kind, ColumnKindTag::Continuous | ColumnKindTag::Binary) {
479 return Err(TermBuilderError::incompatible_config(format!(
480 "bounded() currently supports only numeric columns, got categorical '{name}'"
481 )));
482 }
483 linear_terms.push(LinearTermSpec {
484 name: name.clone(),
485 feature_col: col,
486 feature_cols: vec![col],
487 categorical_levels: vec![],
488 double_penalty: *double_penalty,
489 coefficient_geometry: LinearCoefficientGeometry::Bounded {
490 min: *min,
491 max: *max,
492 prior: prior.clone(),
493 },
494 coefficient_min: None,
495 coefficient_max: None,
496 frozen_function_mass: None,
497 });
498 }
499 ParsedTerm::RandomEffect {
500 name,
501 lenient_unseen,
502 } => {
503 let col = resolve_col(col_map, name)?;
504 random_terms.push(RandomEffectTermSpec {
505 name: name.clone(),
506 feature_col: col,
507 drop_first_level: false,
508 penalized: true,
509 frozen_levels: None,
510 lenient_unseen: *lenient_unseen,
518 });
519 }
520 ParsedTerm::Smooth {
521 label,
522 vars,
523 kind,
524 options,
525 } => {
526 let smooth_vars = vars.clone();
527 let by_name = options.get("by").cloned();
528 let cols = smooth_vars
538 .iter()
539 .map(|v| resolve_col(col_map, v))
540 .collect::<Result<Vec<_>, _>>()?;
541 let mut inner_options = options.clone();
542 inner_options.remove("by");
543 inner_options.remove("ordered");
547 let shape = match inner_options.remove("shape") {
553 None => ShapeConstraint::None,
554 Some(raw) => crate::smooth::parse_shape_constraint(&raw)
555 .map_err(TermBuilderError::invalid_option)?,
556 };
557 if let Some(by_name) = by_name.as_deref() {
562 let by_col = resolve_col(col_map, by_name)?;
563 inject_by_level_sizing_rows(&mut inner_options, ds, by_col);
564 }
565 let inner_basis = build_smooth_basis(
566 *kind,
567 &smooth_vars,
568 &cols,
569 &inner_options,
570 ds,
571 inference_notes,
572 policy,
573 smooth_coordinate_count,
574 )?;
575 if let Some(by_name) = by_name {
585 let by_col = resolve_col(col_map, &by_name)?;
586 match ds.column_kinds.get(by_col).copied().ok_or_else(|| {
587 format!("internal column-kind lookup failed for by variable '{by_name}'")
588 })? {
589 ColumnKindTag::Categorical => {
590 let levels = encoded_levels_for_column(ds, ColIdx::new(by_col));
591 let penalized_group_owner_present =
604 terms.iter().any(|other| match other {
605 ParsedTerm::RandomEffect { name, .. } => name == &by_name,
606 ParsedTerm::Linear {
607 name,
608 explicit: false,
609 ..
610 } if name == &by_name => col_map
611 .get(name)
612 .and_then(|c| ds.column_kinds.get(*c).copied())
613 .map(|kind| matches!(kind, ColumnKindTag::Categorical))
614 .unwrap_or(false),
615 _ => false,
616 });
617 if !random_terms.iter().any(|rt| rt.name == by_name)
628 && !penalized_group_owner_present
629 {
630 random_terms.push(RandomEffectTermSpec {
631 name: by_name.clone(),
632 feature_col: by_col,
633 drop_first_level: true,
634 penalized: false,
635 frozen_levels: None,
636 lenient_unseen: false,
641 });
642 }
643 for (level_bits, level_label) in levels {
650 smooth_terms.push(SmoothTermSpec {
651 frozen_parametric_residualization: None,
652 name: format!("{label}:by={by_name}[{level_label}]"),
653 basis: SmoothBasisSpec::ByVariable {
654 inner: Box::new(inner_basis.clone()),
655 by_col,
656 kind: BySmoothKind::Level { level_bits },
657 by: ByVariableSpec::Level {
658 value_bits: level_bits,
659 label: level_label,
660 },
661 },
662 shape: shape.clone(),
663 joint_null_rotation: None,
664 });
665 }
666 }
667 ColumnKindTag::Binary | ColumnKindTag::Continuous => {
668 smooth_terms.push(SmoothTermSpec {
669 frozen_parametric_residualization: None,
670 name: label.clone(),
671 basis: SmoothBasisSpec::ByVariable {
672 inner: Box::new(inner_basis),
673 by_col,
674 kind: BySmoothKind::Numeric,
675 by: ByVariableSpec::Numeric,
676 },
677 shape,
678 joint_null_rotation: None,
679 });
680 }
681 }
682 } else {
683 smooth_terms.push(SmoothTermSpec {
684 frozen_parametric_residualization: None,
685 name: label.clone(),
686 basis: inner_basis,
687 shape,
688 joint_null_rotation: None,
689 });
690 }
691 }
692 ParsedTerm::LinkWiggle { .. }
693 | ParsedTerm::TimeWiggle { .. }
694 | ParsedTerm::LinkConfig { .. }
695 | ParsedTerm::SurvivalConfig { .. } => {
696 }
698 ParsedTerm::LogSlopeSurface { .. } => {
699 return Err(TermBuilderError::malformed_formula(
700 "logslope(...) declarations must be resolved by the marginal-slope formula path before building a term spec",
701 ));
702 }
703 ParsedTerm::Interaction {
704 vars,
705 double_penalty,
706 } => {
707 let main_effect_present = |target: &str| -> bool {
740 terms.iter().any(|other| match other {
741 ParsedTerm::Linear { name, .. }
742 | ParsedTerm::BoundedLinear { name, .. }
743 | ParsedTerm::RandomEffect { name, .. } => name == target,
744 _ => false,
745 })
746 };
747 let parent_present = |drop_var: &str| -> bool {
753 vars.iter()
754 .filter(|v| v.as_str() != drop_var)
755 .all(|v| main_effect_present(v))
756 };
757
758 let mut numeric_cols = Vec::<usize>::new();
759 let mut categorical_factors =
762 Vec::<(String, usize, Vec<(u64, String)>, bool)>::new();
763 for var in vars {
764 let col = resolve_col(col_map, var)?;
765 let kind = ds.column_kinds.get(col).copied().ok_or_else(|| {
766 TermBuilderError::missing_column(format!(
767 "internal column-kind lookup failed for '{var}'"
768 ))
769 .to_string()
770 })?;
771 match kind {
772 ColumnKindTag::Continuous | ColumnKindTag::Binary => numeric_cols.push(col),
773 ColumnKindTag::Categorical => {
774 let mut levels = encoded_levels_for_column(ds, ColIdx::new(col));
775 let treatment_coded = parent_present(var);
779 if treatment_coded && levels.len() > 1 {
780 levels.remove(0);
781 }
782 if levels.is_empty() {
783 return Err(TermBuilderError::incompatible_config(format!(
784 "interaction `{}` references categorical column `{var}` with no usable levels",
785 vars.join(":")
786 )));
787 }
788 categorical_factors.push((var.clone(), col, levels, treatment_coded));
789 }
790 }
791 }
792
793 let label = vars.join(":");
794
795 if categorical_factors.is_empty() {
796 linear_terms.push(LinearTermSpec {
799 name: label,
800 feature_col: numeric_cols[0],
801 feature_cols: numeric_cols,
802 categorical_levels: vec![],
803 double_penalty: *double_penalty,
805 coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
806 coefficient_min: None,
807 coefficient_max: None,
808 frozen_function_mass: None,
809 });
810 inference_notes.push(format!(
811 "wired linear interaction `{}` as product of numeric columns",
812 vars.join(":")
813 ));
814 } else {
815 let mut cells: Vec<Vec<(usize, u64, String)>> = vec![Vec::new()];
820 for (_var, col, levels, _treatment_coded) in &categorical_factors {
821 let mut next = Vec::with_capacity(cells.len() * levels.len());
822 for cell in &cells {
823 for (bits, level_label) in levels {
824 let mut extended = cell.clone();
825 extended.push((*col, *bits, level_label.clone()));
826 next.push(extended);
827 }
828 }
829 cells = next;
830 }
831
832 let any_dummy_coded = categorical_factors
844 .iter()
845 .any(|(_, _, _, treatment_coded)| !*treatment_coded);
846 if numeric_cols.is_empty() && any_dummy_coded {
847 let reference_cell: Vec<(usize, u64)> = categorical_factors
850 .iter()
851 .map(|(_, col, _, _)| {
852 let levels = encoded_levels_for_column(ds, ColIdx::new(*col));
853 (*col, levels[0].0)
854 })
855 .collect();
856 cells.retain(|cell| {
857 !reference_cell.iter().all(|(rcol, rbits)| {
858 cell.iter()
859 .any(|(col, bits, _)| col == rcol && bits == rbits)
860 })
861 });
862 }
863
864 let n_cells = cells.len();
865 for cell in cells {
866 let cell_suffix = cell
867 .iter()
868 .map(|(_, _, level_label)| level_label.as_str())
869 .collect::<Vec<_>>()
870 .join(":");
871 let categorical_levels =
872 cell.iter().map(|(col, bits, _)| (*col, *bits)).collect();
873 let feature_col = numeric_cols
879 .first()
880 .copied()
881 .unwrap_or(categorical_factors[0].1);
882 linear_terms.push(LinearTermSpec {
883 name: format!("{label}:{cell_suffix}"),
884 feature_col,
885 feature_cols: numeric_cols.clone(),
886 categorical_levels,
887 double_penalty: *double_penalty,
888 coefficient_geometry: LinearCoefficientGeometry::Unconstrained,
889 coefficient_min: None,
890 coefficient_max: None,
891 frozen_function_mass: None,
892 });
893 }
894 let all_treatment_coded = !any_dummy_coded;
895 let coding = if all_treatment_coded {
896 "treatment-coded"
897 } else {
898 "marginality-aware (full dummy / saturated)"
899 };
900 inference_notes.push(format!(
901 "wired factor-aware linear interaction `{}` as {} {} cell column(s)",
902 vars.join(":"),
903 n_cells,
904 coding
905 ));
906 }
907 }
908 }
909 }
910
911 Ok(TermCollectionSpec {
912 linear_terms,
913 random_effect_terms: random_terms,
914 smooth_terms,
915 })
916}
917
918fn split_list_option(raw: &str) -> Vec<String> {
919 let t = raw.trim();
920 let inner = t
927 .strip_prefix('[')
928 .and_then(|u| u.strip_suffix(']'))
929 .or_else(|| {
930 t.strip_prefix("c(")
931 .or_else(|| t.strip_prefix("C("))
932 .or_else(|| t.strip_prefix('('))
933 .and_then(|u| u.strip_suffix(')'))
934 })
935 .unwrap_or(t);
936 inner
937 .split(',')
938 .map(|v| v.trim().to_string())
939 .filter(|v| !v.is_empty())
940 .collect()
941}
942
943fn parse_numeric_expr(raw: &str) -> Result<f64, String> {
944 let mut acc = 1.0f64;
945 let normalized = raw.replace(' ', "");
946 if normalized.eq_ignore_ascii_case("none") {
947 return Err("None is not numeric".to_string());
948 }
949 for factor in normalized.split('*') {
950 if factor.is_empty() {
951 return Err(format!("invalid numeric expression '{raw}'"));
952 }
953 let value = if factor.eq_ignore_ascii_case("pi") || factor == "π" {
954 std::f64::consts::PI
955 } else if factor.eq_ignore_ascii_case("tau") || factor == "τ" {
956 std::f64::consts::TAU
957 } else if let Some(prefix) = factor
958 .strip_suffix("pi")
959 .or_else(|| factor.strip_suffix("π"))
960 {
961 let coefficient = if prefix.is_empty() {
962 1.0
963 } else {
964 prefix
965 .parse::<f64>()
966 .map_err(|err| format!("invalid numeric expression '{raw}': {err}"))?
967 };
968 coefficient * std::f64::consts::PI
969 } else if let Some(prefix) = factor
970 .strip_suffix("tau")
971 .or_else(|| factor.strip_suffix("τ"))
972 {
973 let coefficient = if prefix.is_empty() {
974 1.0
975 } else {
976 prefix
977 .parse::<f64>()
978 .map_err(|err| format!("invalid numeric expression '{raw}': {err}"))?
979 };
980 coefficient * std::f64::consts::TAU
981 } else {
982 factor
983 .parse::<f64>()
984 .map_err(|err| format!("invalid numeric expression '{raw}': {err}"))?
985 };
986 acc *= value;
987 }
988 Ok(acc)
989}
990
991fn option_numeric_expr(
1001 options: &BTreeMap<String, String>,
1002 key: &str,
1003) -> Result<Option<f64>, String> {
1004 match options.get(key) {
1005 None => Ok(None),
1006 Some(raw) => parse_numeric_expr(raw)
1007 .map(Some)
1008 .map_err(|err| format!("option `{key}={raw}` is not a valid numeric value: {err}")),
1009 }
1010}
1011
1012fn parse_periods_option(
1013 options: &BTreeMap<String, String>,
1014 dim: usize,
1015) -> Result<Option<Vec<Option<f64>>>, String> {
1016 let Some(raw) = options.get("period") else {
1017 return Ok(None);
1018 };
1019 let values = split_list_option(raw);
1020 let mut periods = vec![None; dim];
1021 if values.len() == 1 && dim == 1 {
1022 periods[0] = Some(parse_numeric_expr(&values[0])?);
1023 } else {
1024 if values.len() != dim {
1025 return Err(format!(
1026 "period list length {} must match smooth dimension {}",
1027 values.len(),
1028 dim
1029 ));
1030 }
1031 for (i, v) in values.iter().enumerate() {
1032 if v.eq_ignore_ascii_case("none") {
1033 continue;
1034 }
1035 periods[i] = Some(parse_numeric_expr(v)?);
1036 }
1037 }
1038 Ok(Some(periods))
1039}
1040
1041fn parse_periodic_axes_option(
1042 options: &BTreeMap<String, String>,
1043 dim: usize,
1044) -> Result<Option<Vec<Option<f64>>>, String> {
1045 let Some(raw_axes) = options.get("periodic").or_else(|| options.get("cyclic")) else {
1048 let declared = parse_periods_option(options, dim)?;
1055 return Ok(match declared {
1056 Some(periods) if periods.iter().any(Option::is_some) => Some(periods),
1057 _ => None,
1058 });
1059 };
1060 let mut periods = parse_periods_option(options, dim)?.unwrap_or_else(|| vec![None; dim]);
1061 let lowered = raw_axes.trim().to_ascii_lowercase();
1070 if matches!(lowered.as_str(), "true" | "yes" | "y") {
1071 return Ok(Some(periods));
1072 }
1073 if matches!(lowered.as_str(), "false" | "no" | "n") {
1083 return Ok(None);
1084 }
1085 let axes = split_list_option(raw_axes);
1086 if axes.is_empty() {
1087 return Ok(Some(periods));
1088 }
1089
1090 let is_bool = |t: &str| {
1098 matches!(
1099 t.to_ascii_lowercase().as_str(),
1100 "true" | "yes" | "y" | "false" | "no" | "n"
1101 )
1102 };
1103 let is_truthy = |t: &str| matches!(t.to_ascii_lowercase().as_str(), "true" | "yes" | "y");
1104
1105 if axes.len() == 1 && is_bool(&axes[0]) {
1107 if !is_truthy(&axes[0]) {
1108 return Ok(None);
1111 }
1112 return Ok(Some(periods));
1115 }
1116
1117 if axes.iter().all(|a| is_bool(a)) {
1119 if axes.len() != dim {
1120 return Err(format!(
1121 "periodic flag list length {} must match smooth dimension {dim}",
1122 axes.len()
1123 ));
1124 }
1125 if !axes.iter().any(|a| is_truthy(a)) {
1126 return Ok(None);
1127 }
1128 for (i, a) in axes.iter().enumerate() {
1129 if !is_truthy(a) {
1130 periods[i] = None;
1131 }
1132 }
1133 return Ok(Some(periods));
1134 }
1135
1136 for a in &axes {
1139 let axis = a
1140 .parse::<usize>()
1141 .map_err(|err| format!("invalid periodic axis '{a}': {err}"))?;
1142 if axis >= dim {
1143 return Err(format!(
1144 "periodic axis {axis} out of range for {dim}D smooth"
1145 ));
1146 }
1147 if periods[axis].is_none() {
1148 return Err(format!(
1149 "periodic axis {axis} requires period[{axis}] to be finite"
1150 ));
1151 }
1152 }
1153 let listed: std::collections::BTreeSet<usize> = axes
1155 .iter()
1156 .filter_map(|a| a.parse::<usize>().ok())
1157 .collect();
1158 for i in 0..dim {
1159 if !listed.contains(&i) {
1160 periods[i] = None;
1161 }
1162 }
1163 Ok(Some(periods))
1164}
1165
1166fn parse_option_list(raw: &str) -> Vec<String> {
1171 let trimmed = raw.trim();
1172 let inner = trimmed
1178 .strip_prefix('[')
1179 .and_then(|v| v.strip_suffix(']'))
1180 .or_else(|| {
1181 trimmed
1182 .strip_prefix("c(")
1183 .or_else(|| trimmed.strip_prefix("C("))
1184 .or_else(|| trimmed.strip_prefix('('))
1185 .and_then(|v| v.strip_suffix(')'))
1186 })
1187 .unwrap_or(trimmed);
1188 inner
1189 .split(',')
1190 .map(|v| {
1191 v.trim()
1192 .trim_matches('"')
1193 .trim_matches('\'')
1194 .to_ascii_lowercase()
1195 })
1196 .filter(|v| !v.is_empty())
1197 .collect()
1198}
1199
1200fn axes_with_declared_period(
1217 options: &BTreeMap<String, String>,
1218 dim: usize,
1219) -> Result<Vec<bool>, String> {
1220 let mut axes = vec![false; dim];
1221 if let Some(raw) = options.get("period").or_else(|| options.get("periods")) {
1222 let values = split_list_option(raw);
1223 if values.len() == dim {
1224 for (axis, value) in values.iter().enumerate() {
1225 if !value.trim().eq_ignore_ascii_case("none") {
1226 axes[axis] = true;
1227 }
1228 }
1229 }
1230 }
1231 if dim == 1
1235 && PERIOD_ENDPOINT_OPTION_KEYS
1236 .iter()
1237 .any(|key| options.contains_key(*key))
1238 {
1239 axes[0] = true;
1240 }
1241 Ok(axes)
1242}
1243
1244const PERIOD_ENDPOINT_OPTION_KEYS: [&str; 4] = ["period_start", "period_end", "start", "end"];
1246
1247const PERIOD_LENGTH_OPTION_KEYS: [&str; 2] = ["period", "periods"];
1249
1250const PERIOD_ORIGIN_OPTION_KEYS: [&str; 5] = [
1252 "origin",
1253 "origins",
1254 "period_origin",
1255 "period-origin",
1256 "domain_origin",
1257];
1258
1259fn reject_unconsumable_period_declaration(
1271 term_name: &str,
1272 options: &BTreeMap<String, String>,
1273 periodic_axes: &[bool],
1274) -> Result<(), String> {
1275 if periodic_axes.iter().any(|periodic| *periodic) {
1276 return Ok(());
1277 }
1278 let dim = periodic_axes.len();
1279 if let Some(key) = PERIOD_LENGTH_OPTION_KEYS
1280 .iter()
1281 .find(|key| options.contains_key(**key))
1282 {
1283 let hint = if dim > 1 {
1284 format!(
1285 "a scalar `{key}=` does not say which of the {dim} margins wraps; write one entry \
1286 per margin (e.g. {key}=[<value>, None]) or name the axis with periodic=<axis>"
1287 )
1288 } else {
1289 "declare it on a periodic axis or drop it".to_string()
1290 };
1291 return Err(TermBuilderError::invalid_option(format!(
1292 "{term_name}(): `{key}=` declares a period, but no axis of this smooth is periodic — {hint}"
1293 ))
1294 .to_string());
1295 }
1296 if let Some(key) = PERIOD_ORIGIN_OPTION_KEYS
1297 .iter()
1298 .find(|key| options.contains_key(**key))
1299 {
1300 return Err(TermBuilderError::invalid_option(format!(
1301 "{term_name}(): `{key}=` places the start of a periodic domain, but this smooth \
1302 declares no period; add period=<value> or drop it"
1303 ))
1304 .to_string());
1305 }
1306 if let Some(key) = PERIOD_ENDPOINT_OPTION_KEYS
1307 .iter()
1308 .find(|key| options.contains_key(**key))
1309 {
1310 return Err(TermBuilderError::invalid_option(format!(
1311 "{term_name}(): `{key}=` declares a periodic domain endpoint, but no axis of this \
1312 smooth is periodic; on a tensor smooth use periods=[...] with origins=[...], which \
1313 name their margin"
1314 ))
1315 .to_string());
1316 }
1317 Ok(())
1318}
1319
1320fn reject_unconsumable_radial_period_declaration(
1342 term_name: &str,
1343 options: &BTreeMap<String, String>,
1344 dim: usize,
1345 periodic: Option<&[Option<f64>]>,
1346 boundary_is_cyclic: bool,
1347) -> Result<(), String> {
1348 if dim > 1
1349 && let Some(key) = PERIOD_ENDPOINT_OPTION_KEYS
1350 .iter()
1351 .find(|key| options.contains_key(**key))
1352 {
1353 return Err(TermBuilderError::invalid_option(format!(
1354 "{term_name}(): `{key}=` names one axis's periodic domain and is only read on a \
1355 one-dimensional radial smooth; this one has {dim} covariates, so give the wrap as \
1356 period=[…] with one entry per axis"
1357 ))
1358 .to_string());
1359 }
1360 let any_axis_wraps = boundary_is_cyclic
1361 || periodic.is_some_and(|axes| {
1362 (dim == 1 && !axes.is_empty()) || axes.iter().any(Option::is_some)
1363 });
1364 if any_axis_wraps {
1365 return Ok(());
1366 }
1367 let declared = ["periodic", "cyclic"]
1368 .iter()
1369 .chain(PERIOD_LENGTH_OPTION_KEYS.iter())
1370 .chain(PERIOD_ENDPOINT_OPTION_KEYS.iter())
1371 .find(|key| options.contains_key(**key));
1372 let Some(key) = declared else {
1373 return Ok(());
1374 };
1375 if matches!(*key, "periodic" | "cyclic")
1378 && options
1379 .get(*key)
1380 .map(|raw| raw.trim().to_ascii_lowercase())
1381 .is_some_and(|raw| matches!(raw.as_str(), "false" | "no" | "n"))
1382 {
1383 return Ok(());
1384 }
1385 Err(TermBuilderError::invalid_option(format!(
1386 "{term_name}(): `{key}=` declares periodicity, but no axis of this smooth ends up \
1387 periodic. A radial smooth derives its wrap from the center lattice only in one \
1388 dimension (this one has {dim}), so name the period per axis: \
1389 period=[<value>, None, …]"
1390 ))
1391 .to_string())
1392}
1393
1394fn parse_periodic_axes(
1395 options: &BTreeMap<String, String>,
1396 dim: usize,
1397) -> Result<Vec<bool>, String> {
1398 let mut axes = vec![false; dim];
1399 let mut explicitly_aperiodic = false;
1403 if let Some(raw) = options.get("periodic").or_else(|| options.get("cyclic")) {
1404 let lowered = raw.trim().to_ascii_lowercase();
1405 if matches!(lowered.as_str(), "true" | "yes" | "y") {
1406 axes.fill(true);
1407 } else if matches!(lowered.as_str(), "false" | "no" | "n") {
1408 explicitly_aperiodic = true;
1409 } else {
1410 for axis_raw in parse_option_list(raw) {
1411 let axis = axis_raw
1412 .parse::<usize>()
1413 .map_err(|err| format!("invalid periodic axis '{axis_raw}': {err}"))?;
1414 if axis >= dim {
1415 return Err(format!(
1416 "periodic axis {axis} out of range for {dim}D smooth"
1417 ));
1418 }
1419 axes[axis] = true;
1420 }
1421 }
1422 }
1423 if !explicitly_aperiodic
1424 && let Some(raw) = options.get("boundary").or_else(|| options.get("bc"))
1425 {
1426 let boundary = parse_option_list(raw);
1427 if boundary.len() == dim {
1428 for (axis, value) in boundary.iter().enumerate() {
1429 if matches!(value.as_str(), "periodic" | "cyclic" | "cc") {
1430 axes[axis] = true;
1431 }
1432 }
1433 } else if dim == 1
1434 && matches!(
1435 boundary.first().map(String::as_str),
1436 Some("periodic" | "cyclic" | "cc")
1437 )
1438 {
1439 axes[0] = true;
1440 }
1441 }
1442 fold_in_declared_periods(options, dim, &mut axes, explicitly_aperiodic)?;
1443 Ok(axes)
1444}
1445
1446fn fold_in_declared_periods(
1452 options: &BTreeMap<String, String>,
1453 dim: usize,
1454 axes: &mut [bool],
1455 explicitly_aperiodic: bool,
1456) -> Result<(), String> {
1457 let declared = axes_with_declared_period(options, dim)?;
1458 if explicitly_aperiodic && declared.iter().any(|d| *d) {
1459 return Err(TermBuilderError::incompatible_config(
1460 "periodic=false denies the periodicity that the smooth's own period declaration \
1461 asserts; drop one of the two",
1462 )
1463 .to_string());
1464 }
1465 for (axis, declared_axis) in declared.into_iter().enumerate() {
1466 axes[axis] |= declared_axis;
1467 }
1468 Ok(())
1469}
1470
1471fn parse_optional_numeric_list(
1472 options: &BTreeMap<String, String>,
1473 keys: &[&str],
1474 dim: usize,
1475) -> Result<Vec<Option<f64>>, String> {
1476 let Some(raw) = keys.iter().find_map(|key| options.get(*key)) else {
1477 return Ok(vec![None; dim]);
1478 };
1479 let values = split_list_option(raw);
1480 let mut out = vec![None; dim];
1481 if values.len() == 1 && dim == 1 {
1482 if !values[0].eq_ignore_ascii_case("none") {
1483 out[0] = Some(parse_numeric_expr(&values[0])?);
1484 }
1485 return Ok(out);
1486 }
1487 if values.len() != dim {
1488 return Err(format!(
1489 "numeric option list length {} must match smooth dimension {}",
1490 values.len(),
1491 dim
1492 ));
1493 }
1494 for (i, value) in values.iter().enumerate() {
1495 if !value.eq_ignore_ascii_case("none") {
1496 out[i] = Some(parse_numeric_expr(value)?);
1497 }
1498 }
1499 Ok(out)
1500}
1501
1502fn parse_periods(
1503 options: &BTreeMap<String, String>,
1504 periodic_axes: &[bool],
1505) -> Result<Vec<Option<f64>>, String> {
1506 let dim = periodic_axes.len();
1507 let lone_periodic_broadcast = options
1512 .get("period")
1513 .or_else(|| options.get("periods"))
1514 .and_then(|raw| {
1515 let values = split_list_option(raw);
1516 if values.len() != 1 || dim <= 1 {
1517 return None;
1518 }
1519 let mut iter = periodic_axes.iter().enumerate().filter(|(_, p)| **p);
1520 let first = iter.next()?;
1521 if iter.next().is_some() {
1522 return None;
1523 }
1524 Some((first.0, values.into_iter().next()?))
1525 });
1526 let periods = if let Some((axis, value)) = lone_periodic_broadcast {
1527 let mut out = vec![None; dim];
1528 if !value.eq_ignore_ascii_case("none") {
1529 out[axis] = Some(parse_numeric_expr(&value)?);
1530 }
1531 out
1532 } else {
1533 parse_optional_numeric_list(options, &["period", "periods"], dim)?
1534 };
1535 for (axis, (periodic, period)) in periodic_axes.iter().zip(periods.iter()).enumerate() {
1536 if *periodic
1537 && let Some(value) = period
1538 && (!value.is_finite() || *value <= 0.0)
1539 {
1540 return Err(format!(
1541 "period for periodic axis {axis} must be finite and positive, got {value}"
1542 ));
1543 }
1544 }
1545 Ok(periods)
1546}
1547
1548fn parse_period_origins(
1549 options: &BTreeMap<String, String>,
1550 periodic_axes: &[bool],
1551) -> Result<Vec<Option<f64>>, String> {
1552 parse_optional_numeric_list(
1553 options,
1554 &[
1555 "origin",
1556 "origins",
1557 "period_origin",
1558 "period-origin",
1559 "domain_origin",
1560 ],
1561 periodic_axes.len(),
1562 )
1563}
1564
1565fn parse_tensor_periodic_axes(
1576 options: &BTreeMap<String, String>,
1577 dim: usize,
1578) -> Result<Vec<bool>, String> {
1579 let mut axes = vec![false; dim];
1580 if let Some(raw) = options.get("periodic").or_else(|| options.get("cyclic")) {
1581 let lowered = raw.trim().to_ascii_lowercase();
1582 match lowered.as_str() {
1583 "true" | "yes" | "y" => {
1584 axes.fill(true);
1585 }
1586 "false" | "no" | "n" => {
1587 }
1589 _ => {
1590 let entries = parse_option_list(raw);
1591 let all_bool = !entries.is_empty()
1592 && entries.iter().all(|v| {
1593 matches!(
1594 v.as_str(),
1595 "true" | "yes" | "y" | "false" | "no" | "n" | "none"
1596 )
1597 });
1598 let all_zero_one =
1610 !entries.is_empty() && entries.iter().all(|v| v == "0" || v == "1");
1611 let has_repeat = {
1612 let mut seen = std::collections::BTreeSet::new();
1613 !entries.iter().all(|v| seen.insert(v.clone()))
1614 };
1615 let numeric_mask = all_zero_one && entries.len() == dim && has_repeat;
1616 if all_bool || numeric_mask {
1617 if entries.len() != dim {
1618 return Err(format!(
1619 "periodic list length {} must match smooth dimension {}",
1620 entries.len(),
1621 dim
1622 ));
1623 }
1624 for (i, v) in entries.iter().enumerate() {
1625 axes[i] = matches!(v.as_str(), "true" | "yes" | "y" | "1");
1626 }
1627 } else {
1628 for axis_raw in entries {
1629 let axis = axis_raw
1630 .parse::<usize>()
1631 .map_err(|err| format!("invalid periodic axis '{axis_raw}': {err}"))?;
1632 if axis >= dim {
1633 return Err(format!(
1634 "periodic axis {axis} out of range for {dim}D smooth"
1635 ));
1636 }
1637 axes[axis] = true;
1638 }
1639 }
1640 }
1641 }
1642 }
1643 if let Some(raw) = options.get("boundary").or_else(|| options.get("bc")) {
1644 let boundary = parse_option_list(raw);
1645 if boundary.len() == 1 {
1648 if matches!(boundary[0].as_str(), "periodic" | "cyclic" | "cc") {
1649 axes.fill(true);
1650 }
1651 } else if boundary.len() == dim {
1652 for (axis, value) in boundary.iter().enumerate() {
1653 if matches!(value.as_str(), "periodic" | "cyclic" | "cc") {
1654 axes[axis] = true;
1655 }
1656 }
1657 }
1658 }
1659 if let Some(raw) = options.get("bs").or_else(|| options.get("type"))
1667 && bs_selector_is_vector(raw)
1668 {
1669 let per_margin = parse_option_list(raw);
1670 if per_margin.len() == dim {
1671 for (axis, margin_bs) in per_margin.iter().enumerate() {
1672 if matches!(canonicalize_smooth_type(margin_bs), "cc" | "cp" | "cyclic") {
1673 axes[axis] = true;
1674 }
1675 }
1676 }
1677 }
1678 let explicitly_aperiodic = options
1682 .get("periodic")
1683 .or_else(|| options.get("cyclic"))
1684 .is_some_and(|raw| {
1685 matches!(
1686 raw.trim().to_ascii_lowercase().as_str(),
1687 "false" | "no" | "n"
1688 )
1689 });
1690 fold_in_declared_periods(options, dim, &mut axes, explicitly_aperiodic)?;
1691 Ok(axes)
1692}
1693
1694fn validate_tensor_boundary_tokens(
1716 options: &BTreeMap<String, String>,
1717 dim: usize,
1718) -> Result<(), String> {
1719 let Some(raw) = options.get("boundary").or_else(|| options.get("bc")) else {
1720 return Ok(());
1721 };
1722 let entries = parse_option_list(raw);
1723 if entries.len() != 1 && entries.len() != dim {
1729 return Err(TermBuilderError::invalid_option(format!(
1730 "tensor smooth bc/boundary={raw:?} has {} entries but the smooth has {dim} margins; \
1731 pass one token per margin or a single token for all of them",
1732 entries.len()
1733 ))
1734 .to_string());
1735 }
1736 for (axis, value) in entries.iter().enumerate() {
1737 let inert = matches!(
1738 value.trim().to_ascii_lowercase().as_str(),
1739 "clamped" | "open" | "natural" | "free" | "none" | "" | "periodic" | "cyclic" | "cc"
1740 );
1741 if !inert {
1742 return Err(TermBuilderError::unsupported_feature(format!(
1743 "tensor smooth margin {axis} boundary token '{value}' is not supported \
1744 (got bc/boundary={raw:?} on a {dim}-D tensor); tensor margins accept the periodic \
1745 selectors (periodic/cyclic/cc) or the non-periodic markers (clamped/open/natural/free). \
1746 Apply anchored/zero-value endpoint constraints with a 1-D s(x, bc=...) term instead."
1747 ))
1748 .to_string());
1749 }
1750 }
1751 Ok(())
1752}
1753
1754fn tensor_k_axis_option_axis(
1755 key: &str,
1756 cols: &[usize],
1757 ds: &Dataset,
1758) -> Result<Option<usize>, String> {
1759 let Some(suffix) = key.strip_prefix("k_") else {
1760 return Ok(None);
1761 };
1762 if suffix.is_empty() {
1763 return Err("tensor k axis option must be named k_<axis> or k_<variable>".to_string());
1764 }
1765 if let Ok(axis) = suffix.parse::<usize>() {
1766 return if axis < cols.len() {
1767 Ok(Some(axis))
1768 } else {
1769 Err(format!(
1770 "tensor k axis option `{key}` references axis {axis}, but the smooth has {} margins",
1771 cols.len()
1772 ))
1773 };
1774 }
1775
1776 let mut matches = cols
1777 .iter()
1778 .enumerate()
1779 .filter(|(_, col)| ds.headers.get(**col).is_some_and(|name| name == suffix))
1780 .map(|(axis, _)| axis);
1781 let first = matches.next();
1782 if matches.next().is_some() {
1783 return Err(format!(
1784 "tensor k axis option `{key}` matches more than one margin named `{suffix}`"
1785 ));
1786 }
1787 first.map(Some).ok_or_else(|| {
1788 let margin_names = cols
1789 .iter()
1790 .enumerate()
1791 .map(|(axis, col)| {
1792 let name = ds
1793 .headers
1794 .get(*col)
1795 .map(String::as_str)
1796 .unwrap_or("<unnamed>");
1797 format!("{axis}:{name}")
1798 })
1799 .collect::<Vec<_>>()
1800 .join(", ");
1801 format!(
1802 "tensor k axis option `{key}` does not match a margin index or name; tensor margins are [{margin_names}]"
1803 )
1804 })
1805}
1806
1807fn is_tensor_k_axis_option_key(key: &str) -> bool {
1808 key.strip_prefix("k_")
1809 .is_some_and(|suffix| !suffix.is_empty())
1810}
1811
1812fn parse_tensor_k_list(
1816 options: &BTreeMap<String, String>,
1817 cols: &[usize],
1818 ds: &Dataset,
1819) -> Result<(Vec<usize>, bool), String> {
1820 let mut axis_values = vec![None; cols.len()];
1821 let mut saw_axis_alias = false;
1822 for (key, value) in options {
1823 let Some(axis) = tensor_k_axis_option_axis(key, cols, ds)? else {
1824 continue;
1825 };
1826 saw_axis_alias = true;
1827 if axis_values[axis].is_some() {
1828 return Err(format!("tensor k axis {axis} is specified more than once"));
1829 }
1830 let k: usize = value
1831 .parse()
1832 .map_err(|err| format!("invalid tensor k option `{key}={value}`: {err}"))?;
1833 axis_values[axis] = Some(k);
1834 }
1835
1836 let raw = options
1837 .get("k")
1838 .or_else(|| options.get("basis_dim"))
1839 .or_else(|| options.get("basis-dim"))
1840 .or_else(|| options.get("basisdim"));
1841 if saw_axis_alias {
1842 if raw.is_some() {
1843 return Err(
1844 "tensor k axis aliases cannot be combined with k= or basis_dim=".to_string(),
1845 );
1846 }
1847 if let Some(missing_axis) = axis_values.iter().position(Option::is_none) {
1848 let margin_name = cols
1849 .get(missing_axis)
1850 .and_then(|col| ds.headers.get(*col))
1851 .map(String::as_str)
1852 .unwrap_or("<unnamed>");
1853 return Err(format!(
1854 "tensor k axis aliases must specify every margin; missing axis {missing_axis} ({margin_name})"
1855 ));
1856 }
1857 return Ok((
1858 axis_values
1859 .into_iter()
1860 .map(|k| k.expect("missing axis values rejected above"))
1861 .collect(),
1862 false,
1863 ));
1864 }
1865 let Some(raw) = raw else {
1866 let inferred = heuristic_tensor_margin_knots(cols, ds);
1867 return Ok((inferred, true));
1868 };
1869 let entries = split_list_option(raw);
1870 if entries.len() == 1 {
1871 let k: usize = entries[0]
1872 .parse()
1873 .map_err(|err| format!("invalid tensor k '{}': {err}", entries[0]))?;
1874 return Ok((vec![k; cols.len()], false));
1875 }
1876 if entries.len() != cols.len() {
1877 return Err(format!(
1878 "tensor k list length {} must match smooth dimension {}",
1879 entries.len(),
1880 cols.len()
1881 ));
1882 }
1883 let mut out = Vec::with_capacity(entries.len());
1884 for entry in entries {
1885 let k: usize = entry
1886 .parse()
1887 .map_err(|err| format!("invalid tensor k '{entry}': {err}"))?;
1888 out.push(k);
1889 }
1890 Ok((out, false))
1891}
1892
1893fn parse_tensor_identifiability(
1902 options: &BTreeMap<String, String>,
1903 kind: SmoothKind,
1904) -> Result<TensorBSplineIdentifiability, String> {
1905 let Some(raw) = options.get("identifiability").map(String::as_str) else {
1906 return Ok(match kind {
1907 SmoothKind::Ti => TensorBSplineIdentifiability::MarginalSumToZero,
1908 _ => TensorBSplineIdentifiability::default(),
1909 });
1910 };
1911 match raw.trim().to_ascii_lowercase().as_str() {
1912 "none" => Ok(TensorBSplineIdentifiability::None),
1913 "sum_tozero" | "sum-to-zero" | "center_sum_tozero" | "center-sum-to-zero" | "centered"
1914 | "sumtozero" => Ok(TensorBSplineIdentifiability::SumToZero),
1915 "marginal_sum_tozero" | "marginal-sum-to-zero" | "marginal_sumtozero"
1916 | "marginalsumtozero" | "interaction" => {
1917 Ok(TensorBSplineIdentifiability::MarginalSumToZero)
1918 }
1919 other => Err(TermBuilderError::unsupported_feature(format!(
1920 "invalid tensor identifiability '{other}'; expected one of: none, sum_tozero, marginal_sum_tozero"
1921 ))
1922 .to_string()),
1923 }
1924}
1925
1926fn parse_bspline_identifiability(
1950 options: &BTreeMap<String, String>,
1951) -> Result<Option<BSplineIdentifiability>, String> {
1952 let Some(raw) = options.get("identifiability").map(String::as_str) else {
1953 return Ok(None);
1954 };
1955 match raw.trim().to_ascii_lowercase().as_str() {
1956 "none" => Ok(Some(BSplineIdentifiability::None)),
1957 "sum_tozero" | "sum-to-zero" | "center_sum_tozero" | "center-sum-to-zero" | "centered"
1958 | "sumtozero" => Ok(Some(BSplineIdentifiability::WeightedSumToZero {
1959 weights: None,
1960 })),
1961 "linear" | "remove_linear_trend" | "remove-linear-trend" | "removelineartrend"
1962 | "center_linear_orthogonal" | "center-linear-orthogonal" => {
1963 Ok(Some(BSplineIdentifiability::RemoveLinearTrend))
1964 }
1965 "frozen" | "frozen_transform" | "orthogonal" | "orthogonal_to_design_columns" => {
1966 Err(TermBuilderError::unsupported_feature(format!(
1967 "B-spline identifiability '{}' is internal-only (it is minted by design freezing \
1968 or needs an explicit design-column block); use one of: none, sum_tozero, linear",
1969 raw.trim()
1970 ))
1971 .to_string())
1972 }
1973 other => Err(TermBuilderError::unsupported_feature(format!(
1974 "invalid B-spline identifiability '{other}'; expected one of: none, sum_tozero, linear"
1975 ))
1976 .to_string()),
1977 }
1978}
1979
1980#[derive(Debug, Clone, Copy, Default)]
1984struct BSplineIdentifiabilityContext {
1985 has_anchor: bool,
1988 periodic: bool,
1990 natural_cubic_regression: bool,
1994}
1995
1996fn resolve_bspline_identifiability(
2025 options: &BTreeMap<String, String>,
2026 structural_default: BSplineIdentifiability,
2027 context: BSplineIdentifiabilityContext,
2028) -> Result<BSplineIdentifiability, String> {
2029 let Some(explicit) = parse_bspline_identifiability(options)? else {
2030 return Ok(structural_default);
2031 };
2032 if context.has_anchor && !matches!(explicit, BSplineIdentifiability::None) {
2033 return Err(TermBuilderError::incompatible_config(
2034 "an anchored endpoint already fixes the smooth's level (the global intercept is \
2035 suppressed), so it cannot also carry a centering identifiability constraint; \
2036 drop the anchor or use identifiability='none'",
2037 )
2038 .to_string());
2039 }
2040 if matches!(explicit, BSplineIdentifiability::RemoveLinearTrend) {
2041 if context.periodic {
2042 return Err(TermBuilderError::incompatible_config(
2043 "identifiability='linear' removes the constant and linear directions using \
2044 open-knot Greville geometry, which a periodic basis does not span; use 'none' \
2045 or 'sum_tozero' on a periodic smooth",
2046 )
2047 .to_string());
2048 }
2049 if context.natural_cubic_regression {
2050 return Err(TermBuilderError::incompatible_config(
2051 "identifiability='linear' needs B-spline knot/degree geometry, which the natural \
2052 cubic regression basis (bs='cr'/'cs') does not carry; use 'none' or 'sum_tozero', \
2053 or switch to bs='ps'",
2054 )
2055 .to_string());
2056 }
2057 }
2058 Ok(explicit)
2059}
2060
2061fn bspline_boundary_declares_periodic_axis(options: &BTreeMap<String, String>) -> bool {
2062 options
2063 .get("boundary")
2064 .or_else(|| options.get("bc"))
2065 .map(|raw| {
2066 parse_option_list(raw)
2067 .into_iter()
2068 .any(|value| matches!(value.as_str(), "periodic" | "cyclic" | "cc"))
2069 })
2070 .unwrap_or(false)
2071}
2072
2073pub(crate) fn canonicalize_smooth_type(raw: &str) -> &str {
2095 match raw {
2096 "tp" => "tps",
2099 "gp" => "matern",
2104 "curv" | "constant_curvature" | "mkappa" => "curvature",
2108 "mjs" | "measure_jet" | "web" => "measurejet",
2112 other => other,
2113 }
2114}
2115
2116pub(crate) fn tensor_margin_bs_is_supported(margin_bs: &str) -> bool {
2127 matches!(
2128 canonicalize_smooth_type(margin_bs),
2129 "tps" | "ps" | "bs" | "bspline" | "cr" | "cs" | "cc" | "cp" | "cyclic"
2130 )
2131}
2132
2133pub(crate) fn smooth_options_declare_periodic(options: &BTreeMap<String, String>) -> bool {
2139 options.contains_key("periodic")
2140 || options.contains_key("cyclic")
2141 || options
2142 .get("boundary")
2143 .or_else(|| options.get("bc"))
2144 .map(|boundary| {
2145 boundary.to_ascii_lowercase().contains("periodic")
2146 || boundary.to_ascii_lowercase().contains("cyclic")
2147 })
2148 .unwrap_or(false)
2149}
2150
2151pub(crate) fn bs_selector_is_vector(raw: &str) -> bool {
2168 let trimmed = raw.trim();
2169 let bracketed = (trimmed.starts_with('[') && trimmed.ends_with(']'))
2170 || (trimmed.starts_with("c(") || trimmed.starts_with("C(")) && trimmed.ends_with(')')
2171 || (trimmed.starts_with('(') && trimmed.ends_with(')'));
2172 bracketed && !parse_option_list(trimmed).is_empty()
2173}
2174
2175pub fn resolve_smooth_type_name(
2176 kind: SmoothKind,
2177 n_cols: usize,
2178 options: &BTreeMap<String, String>,
2179) -> String {
2180 let selector = options.get("type").or_else(|| options.get("bs"));
2181 if let Some(raw) = selector
2186 && bs_selector_is_vector(raw)
2187 && matches!(kind, SmoothKind::Te | SmoothKind::Ti | SmoothKind::T2)
2188 {
2189 return "tensor".to_string();
2190 }
2191 selector
2192 .map(|s| canonicalize_smooth_type(&s.to_ascii_lowercase()).to_string())
2193 .unwrap_or_else(|| match kind {
2194 SmoothKind::Te | SmoothKind::Ti | SmoothKind::T2 => "tensor".to_string(),
2195 SmoothKind::S if n_cols == 1 => "bspline".to_string(),
2196 SmoothKind::S if smooth_options_declare_periodic(options) => "tensor".to_string(),
2200 SmoothKind::S => "tps".to_string(),
2201 })
2202}
2203
2204pub fn smooth_type_uses_spatial_center_heuristic(canonical_type: &str) -> bool {
2213 matches!(canonical_type, "tps" | "matern" | "duchon")
2214}
2215
2216pub fn build_smooth_basis(
2217 kind: SmoothKind,
2218 vars: &[String],
2219 cols: &[usize],
2220 options: &BTreeMap<String, String>,
2221 ds: &Dataset,
2222 inference_notes: &mut Vec<String>,
2223 policy: &ResourcePolicy,
2224 smooth_coordinate_count: usize,
2225) -> Result<SmoothBasisSpec, String> {
2226 let stripped_sizing_options;
2230 let (options, sizing_rows) = match options.get(DEFAULT_SIZING_ROWS_OPTION) {
2231 Some(raw) => {
2232 let rows = raw.parse::<usize>().map_err(|_| {
2233 format!("internal by-level sizing rows carrier is not a count: '{raw}'")
2234 })?;
2235 let mut cleaned = options.clone();
2236 cleaned.remove(DEFAULT_SIZING_ROWS_OPTION);
2237 stripped_sizing_options = cleaned;
2238 (&stripped_sizing_options, rows)
2239 }
2240 None => (options, ds.values.nrows()),
2241 };
2242 let coord_cols: Vec<(&String, usize)> = vars
2264 .iter()
2265 .zip(cols.iter().copied())
2266 .filter(|(_, col)| !matches!(ds.column_kinds.get(*col), Some(ColumnKindTag::Categorical)))
2267 .collect();
2268 if !coord_cols.is_empty() {
2269 let views: Vec<ArrayView1<'_, f64>> = coord_cols
2270 .iter()
2271 .map(|(_, col)| ds.values.column(*col))
2272 .collect();
2273 let n_rows = views[0].len();
2274 let mut distinct_points = std::collections::HashSet::<Vec<u64>>::new();
2275 for r in 0..n_rows {
2276 let key: Vec<u64> = views
2277 .iter()
2278 .map(|v| gam_data::canonical_level_bits(v[r]))
2279 .collect();
2280 distinct_points.insert(key);
2281 if distinct_points.len() > 1 {
2282 break;
2283 }
2284 }
2285 if distinct_points.len() <= 1 {
2286 return Err(TermBuilderError::degenerate_data(if coord_cols.len() == 1 {
2287 let var = coord_cols[0].0;
2288 format!(
2289 "smooth term over '{var}' has only one unique value in the training data \
2290 — a smooth on a constant column is degenerate and would only fit the response mean. \
2291 Remove `{var}` from the smooth, drop the term, or check the data."
2292 )
2293 } else {
2294 let names = coord_cols
2295 .iter()
2296 .map(|(v, _)| v.as_str())
2297 .collect::<Vec<_>>()
2298 .join(", ");
2299 format!(
2300 "smooth term over ({names}) has only one unique joint coordinate in the training \
2301 data — every coordinate is constant, so the smooth is degenerate and would only \
2302 fit the response mean. Drop the term or check the data."
2303 )
2304 })
2305 .to_string());
2306 }
2307
2308 if matches!(
2317 resolve_smooth_type_name(kind, cols.len(), options).as_str(),
2318 "sphere" | "s2" | "sos"
2319 ) {
2320 for (axis, (var, col)) in coord_cols.iter().enumerate() {
2321 let column = ds.values.column(*col);
2322 let mut distinct = std::collections::HashSet::<u64>::new();
2323 for &value in column.iter() {
2324 distinct.insert(gam_data::canonical_level_bits(value));
2325 if distinct.len() > 1 {
2326 break;
2327 }
2328 }
2329 if distinct.len() <= 1 {
2330 let slice = if axis == 0 {
2333 "a single parallel (constant latitude)"
2334 } else {
2335 "a single meridian (constant longitude)"
2336 };
2337 return Err(TermBuilderError::degenerate_data(format!(
2338 "sphere smooth has a constant '{var}' column — every point lies on \
2339 {slice}, so the 2-sphere term is degenerate and unidentifiable along \
2340 that axis. A spherical smooth needs genuine variation in BOTH latitude \
2341 and longitude; vary '{var}', drop the term, or fit a 1-D smooth on the \
2342 varying coordinate."
2343 ))
2344 .to_string());
2345 }
2346 }
2347 }
2348 }
2349 if let Some(by_name) = options.get("by").cloned() {
2350 let by_col = options
2351 .get("__by_col")
2352 .and_then(|raw| raw.parse::<usize>().ok())
2353 .or_else(|| vars.iter().position(|v| v == &by_name).map(|idx| cols[idx]))
2354 .ok_or_else(|| format!("unknown by= column '{by_name}'"))?;
2355 let mut inner_options = options.clone();
2356 inner_options.remove("by");
2357 inner_options.remove("__by_col");
2358 inner_options.remove("id");
2359 inject_by_level_sizing_rows(&mut inner_options, ds, by_col);
2363 let inner = build_smooth_basis(
2364 kind,
2365 vars,
2366 cols,
2367 &inner_options,
2368 ds,
2369 inference_notes,
2370 policy,
2371 smooth_coordinate_count,
2372 )?;
2373 let by_kind = match ds.column_kinds.get(by_col).copied() {
2374 Some(ColumnKindTag::Categorical) => ByVarKind::Factor {
2375 feature_col: by_col,
2376 ordered: option_bool(options, "ordered").unwrap_or(false),
2377 frozen_levels: None,
2378 },
2379 Some(ColumnKindTag::Continuous | ColumnKindTag::Binary) => ByVarKind::Numeric {
2380 feature_col: by_col,
2381 },
2382 None => {
2383 return Err(format!(
2384 "internal column-kind lookup failed for by='{by_name}'"
2385 ));
2386 }
2387 };
2388 return Ok(SmoothBasisSpec::BySmooth {
2389 smooth: Box::new(inner),
2390 by_kind,
2391 });
2392 }
2393
2394 let smooth_double_penalty = option_bool(options, "double_penalty").unwrap_or(true);
2395 let type_opt = resolve_smooth_type_name(kind, cols.len(), options);
2396
2397 if matches!(type_opt.as_str(), "fs" | "sz" | "re") {
2398 validate_known_options(type_opt.as_str(), options, SHAPE_CONSTRAINED_SMOOTH_OPTION_KEYS)?;
2399 if cols.len() != 2 {
2400 return Err(format!(
2401 "{} factor-smooth currently expects exactly two variables (one numeric, one categorical)",
2402 type_opt
2403 ));
2404 }
2405 let kinds = cols
2406 .iter()
2407 .map(|&c| ds.column_kinds.get(c).copied())
2408 .collect::<Vec<_>>();
2409 let (cont_idx, group_idx) = if type_opt == "re" {
2410 match (kinds[0], kinds[1]) {
2412 (Some(ColumnKindTag::Categorical), _) => (1usize, 0usize),
2413 (_, Some(ColumnKindTag::Categorical)) => (0usize, 1usize),
2414 _ => (1usize, 0usize),
2415 }
2416 } else {
2417 match (kinds[0], kinds[1]) {
2418 (_, Some(ColumnKindTag::Categorical)) => (0usize, 1usize),
2419 (Some(ColumnKindTag::Categorical), _) => (1usize, 0usize),
2420 _ => {
2421 return Err(format!(
2422 "{} factor-smooth requires one categorical factor variable",
2423 type_opt
2424 ));
2425 }
2426 }
2427 };
2428 let c = cols[cont_idx];
2429 let (minv, maxv) = col_minmax(ds.values.column(c))?;
2430 let degree = if type_opt == "re" {
2431 1
2432 } else {
2433 option_usize(options, "degree").unwrap_or(DEFAULT_BSPLINE_DEGREE)
2434 };
2435 let pooled_internal = heuristic_knots_for_column(ds.values.column(c));
2455 let default_internal = if type_opt == "re" {
2456 0
2469 } else {
2470 let min_group_resolution =
2471 min_per_group_unique_count(ds.values.column(c), ds.values.column(cols[group_idx]));
2472 let basis_cap = min_group_resolution.saturating_sub(2).max(degree + 2);
2480 let internal_cap = basis_cap.saturating_sub(degree + 1);
2481 let capped = pooled_internal.min(internal_cap.max(1));
2482 let fs_default_internal = FACTOR_SMOOTH_DEFAULT_BASIS_DIM
2498 .saturating_sub(degree + 1)
2499 .max(1);
2500 capped.min(fs_default_internal)
2501 };
2502 let (n_knots, _, effective_degree) =
2503 parse_ps_internal_knots(options, degree, default_internal)?;
2504 let penalty_order = option_usize(options, "penalty_order")
2505 .unwrap_or(if effective_degree > 1 { 2 } else { 1 })
2506 .min(effective_degree);
2507 let marginal_knotspec = resolve_nonperiodic_bspline_knotspec(
2541 options,
2542 ds.values.column(c),
2543 (minv, maxv),
2544 effective_degree,
2545 n_knots,
2546 )?;
2547 let marginal = BSplineBasisSpec {
2548 degree: effective_degree,
2549 penalty_order,
2550 knotspec: marginal_knotspec,
2551 double_penalty: option_bool(options, "double_penalty")
2562 .unwrap_or(type_opt.as_str() != "sz"),
2563 identifiability: BSplineIdentifiability::None,
2564 boundary_conditions: Default::default(),
2565 boundary: OneDimensionalBoundary::Open,
2566 };
2567 let flavour = match type_opt.as_str() {
2568 "fs" => FactorSmoothFlavour::Fs {
2569 m_null_penalty_orders: vec![
2570 option_usize(options, "m").unwrap_or(DEFAULT_PENALTY_ORDER),
2571 ],
2572 },
2573 "sz" => FactorSmoothFlavour::Sz,
2574 "re" => FactorSmoothFlavour::Re,
2575 other => {
2577 return Err(format!(
2578 "internal: factor-smooth flavour dispatch reached unexpected type `{}`",
2579 other
2580 ));
2581 }
2582 };
2583 return Ok(SmoothBasisSpec::FactorSmooth {
2584 spec: FactorSmoothSpec {
2585 continuous_cols: vec![c],
2586 group_col: cols[group_idx],
2587 marginal,
2588 flavour,
2589 group_frozen_levels: None,
2590 frozen_global_orthogonality: None,
2591 },
2592 });
2593 }
2594
2595 match type_opt.as_str() {
2596 "cyclic" | "cc" | "cp" | "cyclic-ps" | "periodic" => {
2604 validate_known_options("cyclic", options, CYCLIC_SMOOTH_OPTION_KEYS)?;
2605 if cols.len() != 1 {
2606 return Err(format!(
2607 "periodic smooth expects one variable, got {}",
2608 cols.len()
2609 ));
2610 }
2611 let c = cols[0];
2612 let (minv, maxv) = col_minmax(ds.values.column(c))?;
2613 let degree = option_usize(options, "degree").unwrap_or(DEFAULT_BSPLINE_DEGREE);
2614 let mut default_internal = heuristic_knots_for_column(ds.values.column(c));
2615 if ds.values.nrows() <= 32 && smooth_coordinate_count >= 5 {
2616 default_internal = default_internal.min(1);
2617 }
2618 let cyclic_default_basis_cap = CYCLIC_DEFAULT_BASIS_DIM.max(degree + 1);
2634 let default_basis = (default_internal + degree + 1).min(cyclic_default_basis_cap);
2635 let num_basis = option_usize_any(options, &["k", "basis_dim", "basis-dim", "basisdim"])
2636 .unwrap_or(default_basis);
2637 if num_basis < degree + 1 {
2638 return Err(format!(
2639 "periodic smooth: k={} too small for degree {}; expected k >= {}",
2640 num_basis,
2641 degree,
2642 degree + 1
2643 ));
2644 }
2645 let periodic_axes = [true];
2656 let periods = parse_periods(options, &periodic_axes)?;
2657 let origins = parse_period_origins(options, &periodic_axes)?;
2658 let has_endpoint_decl = ["period_start", "start", "period_end", "end"]
2673 .iter()
2674 .any(|key| options.contains_key(*key));
2675 let (domain_start, period) = if let Some(p) = periods[0] {
2676 (origins[0].unwrap_or(minv), p)
2677 } else if has_endpoint_decl {
2678 parse_periodic_domain_1d(options, minv, maxv)?
2679 } else {
2680 let span = maxv - minv;
2681 if !(span.is_finite() && span > 0.0) {
2682 return Err(format!(
2683 "cyclic smooth requires a positive observed data range to derive \
2684 its period, got [{minv}, {maxv}]"
2685 ));
2686 }
2687 (origins[0].unwrap_or(minv), span)
2688 };
2689 let identifiability = resolve_bspline_identifiability(
2696 options,
2697 BSplineIdentifiability::default(),
2698 BSplineIdentifiabilityContext {
2699 periodic: true,
2700 ..Default::default()
2701 },
2702 )?;
2703 Ok(SmoothBasisSpec::BSpline1D {
2704 feature_col: c,
2705 spec: BSplineBasisSpec {
2706 degree,
2707 penalty_order: option_usize(options, "penalty_order")
2708 .unwrap_or(DEFAULT_PENALTY_ORDER),
2709 knotspec: BSplineKnotSpec::PeriodicUniform {
2710 data_range: (domain_start, domain_start + period),
2711 num_basis,
2712 },
2713 double_penalty: smooth_double_penalty,
2714 identifiability,
2715 boundary_conditions: Default::default(),
2716 boundary: OneDimensionalBoundary::Cyclic {
2717 start: domain_start,
2718 end: domain_start + period,
2719 },
2720 },
2721 })
2722 }
2723 "bspline" | "ps" | "p-spline" | "cr" | "cs" => {
2724 let validation_name = match type_opt.as_str() {
2737 "cr" => "cr",
2738 "cs" => "cs",
2739 _ => "bspline",
2740 };
2741 validate_known_options(validation_name, options, BSPLINE_SMOOTH_OPTION_KEYS)?;
2742 if cols.len() != 1 {
2743 return Err(TermBuilderError::incompatible_config(format!(
2744 "bspline smooth expects one variable, got {}",
2745 cols.len()
2746 ))
2747 .to_string());
2748 }
2749 let c = cols[0];
2750 let (minv, maxv) = col_minmax(ds.values.column(c))?;
2751 let degree = option_usize(options, "degree").unwrap_or(DEFAULT_BSPLINE_DEGREE);
2752 let default_internal = heuristic_knots_for_column(ds.values.column(c));
2753 let (mut n_knots, inferred, effective_degree) =
2754 parse_ps_internal_knots(options, degree, default_internal)?;
2755 let periodic_axes = parse_periodic_axes(options, 1).map_err(|e| e.to_string())?;
2756 reject_unconsumable_period_declaration(validation_name, options, &periodic_axes)?;
2760 if periodic_axes[0] && effective_degree != degree {
2765 return Err(TermBuilderError::invalid_option(format!(
2766 "periodic smooth: k={} too small for degree {}; expected k >= {}",
2767 effective_degree + 1,
2768 degree,
2769 degree + 1
2770 ))
2771 .to_string());
2772 }
2773 if inferred && ds.values.nrows() <= 32 && smooth_coordinate_count >= 5 {
2774 n_knots = n_knots.min(1);
2775 }
2776 if inferred {
2777 let unique = unique_count_column(ds.values.column(c));
2778 let ceiling = ((unique as f64).cbrt() as usize).max(20);
2779 inference_notes.push(format!(
2780 "Automatically set {} internal knots for smooth '{}' from {} unique values (rule: clamp(unique/4, 4..max(20, cbrt(unique))) = clamp(unique/4, 4..{})). Override with knots=... or k=....",
2781 n_knots,
2782 vars.join(","),
2783 unique,
2784 ceiling,
2785 ));
2786 }
2787 let boundary_conditions =
2788 if periodic_axes[0] && bspline_boundary_declares_periodic_axis(options) {
2789 BSplineBoundaryConditions::default()
2790 } else {
2791 parse_bspline_boundary_conditions(options).map_err(|e| e.to_string())?
2792 };
2793 let structural_identifiability = if boundary_conditions.has_anchor() {
2803 BSplineIdentifiability::None
2804 } else {
2805 BSplineIdentifiability::default()
2806 };
2807 let identifiability = resolve_bspline_identifiability(
2813 options,
2814 structural_identifiability,
2815 BSplineIdentifiabilityContext {
2816 has_anchor: boundary_conditions.has_anchor(),
2817 periodic: periodic_axes[0],
2818 natural_cubic_regression: !periodic_axes[0]
2819 && (type_opt == "cr" || type_opt == "cs"),
2820 },
2821 )?;
2822 let periods = parse_periods(options, &periodic_axes).map_err(|e| e.to_string())?;
2823 let origins =
2824 parse_period_origins(options, &periodic_axes).map_err(|e| e.to_string())?;
2825 let (knotspec, boundary) = if periodic_axes[0] {
2826 if !boundary_conditions.is_free() {
2827 return Err(TermBuilderError::incompatible_config(
2828 "periodic B-splines cannot also declare endpoint boundary conditions",
2829 )
2830 .to_string());
2831 }
2832 {
2833 let (domain_start, p_value) = if let Some(period) = periods[0] {
2834 (origins[0].unwrap_or(minv), period)
2835 } else {
2836 parse_periodic_domain_1d(options, minv, maxv).map_err(|e| e.to_string())?
2837 };
2838 let domain_end = domain_start + p_value;
2839 (
2840 BSplineKnotSpec::PeriodicUniform {
2841 data_range: (domain_start, domain_end),
2842 num_basis: n_knots + effective_degree + 1,
2843 },
2844 OneDimensionalBoundary::Cyclic {
2845 start: domain_start,
2846 end: domain_end,
2847 },
2848 )
2849 }
2850 } else if type_opt == "cr" || type_opt == "cs" {
2851 let k_cr = (n_knots + effective_degree + 1).max(CR_MIN_KNOTS);
2868 let knotspec = match capped_cr_marginal_knotspec(
2869 ds.values.column(c),
2870 k_cr,
2871 &vars.join(","),
2872 inference_notes,
2873 )? {
2874 Some(cr_knotspec) => cr_knotspec,
2875 None => resolve_nonperiodic_bspline_knotspec(
2876 options,
2877 ds.values.column(c),
2878 (minv, maxv),
2879 effective_degree,
2880 n_knots,
2881 )?,
2882 };
2883 (knotspec, parse_cyclic_boundary(options, minv, maxv)?)
2884 } else {
2885 (
2886 resolve_nonperiodic_bspline_knotspec(
2887 options,
2888 ds.values.column(c),
2889 (minv, maxv),
2890 effective_degree,
2891 n_knots,
2892 )?,
2893 parse_cyclic_boundary(options, minv, maxv)?,
2894 )
2895 };
2896 let double_penalty = smooth_double_penalty;
2900 let penalty_order = option_usize(options, "penalty_order")
2905 .unwrap_or(DEFAULT_PENALTY_ORDER)
2906 .min(effective_degree);
2907 Ok(SmoothBasisSpec::BSpline1D {
2908 feature_col: c,
2909 spec: BSplineBasisSpec {
2910 degree: effective_degree,
2911 penalty_order,
2912 knotspec,
2913 double_penalty,
2914 identifiability,
2915 boundary,
2916 boundary_conditions,
2917 },
2918 })
2919 }
2920 "tps" | "thinplate" | "thin-plate" => {
2921 validate_known_options("thinplate", options, THINPLATE_SMOOTH_OPTION_KEYS)?;
2922 let plan = plan_spatial_basis(
2923 sizing_rows,
2924 cols.len(),
2925 CenterCountRequest::Default,
2926 DuchonNullspaceOrder::Linear,
2927 option_bool(options, "scale_dims").unwrap_or(false),
2928 policy,
2929 )
2930 .map_err(|e| e.to_string())?;
2931 let default_centers = plan.centers;
2941 let centers = parse_countwith_basis_alias(
2942 options,
2943 "centers",
2944 cap_default_spatial_centers(options, default_centers),
2945 )?;
2946 let center_strategy = if has_explicit_countwith_basis_alias(options, "centers") {
2947 spatial_center_strategy_for_dimension(centers, cols.len())
2948 } else {
2949 auto_spatial_center_strategy(centers, cols.len())
2950 };
2951 let periodic = parse_periodic_axes_option(options, cols.len())?;
2952 reject_unconsumable_radial_period_declaration(
2953 "thinplate",
2954 options,
2955 cols.len(),
2956 periodic.as_deref(),
2957 false,
2958 )?;
2959 Ok(SmoothBasisSpec::ThinPlate {
2960 feature_cols: cols.to_vec(),
2961 spec: ThinPlateBasisSpec {
2962 center_strategy,
2963 periodic,
2964 length_scale: option_f64(options, "length_scale").unwrap_or(0.0),
2972 double_penalty: smooth_double_penalty,
2973 identifiability: parse_spatial_identifiability(options)
2974 .map_err(|e| e.to_string())?,
2975 radial_reparam: None,
2976 },
2977 input_scale: None,
2978 })
2979 }
2980 "sphere" | "s2" | "sos" => {
2981 validate_known_options("sphere", options, SPHERE_SMOOTH_OPTION_KEYS)?;
2982 if cols.len() != 2 {
2983 return Err(format!(
2984 "sphere smooth expects exactly two variables (lat, lon), got {}",
2985 cols.len()
2986 ));
2987 }
2988 let radians = option_bool(options, "radians").unwrap_or_else(|| {
2989 options
2990 .get("units")
2991 .map(|u| u.eq_ignore_ascii_case("radian") || u.eq_ignore_ascii_case("radians"))
2992 .unwrap_or(false)
2993 });
2994 let degree_requested = options.contains_key("degree")
3000 || options.contains_key("l")
3001 || options.contains_key("max_degree")
3002 || options.contains_key("max-degree");
3003 let kernel = options
3004 .get("kernel")
3005 .or_else(|| options.get("method"))
3006 .map(|raw| strip_quotes(raw).trim().to_ascii_lowercase())
3007 .unwrap_or_else(|| {
3008 if degree_requested {
3009 "harmonic".to_string()
3010 } else {
3011 "sobolev".to_string()
3012 }
3013 });
3014 let (method, wahba_kernel) = match kernel.as_str() {
3015 "sobolev" | "wahba" | "wahba_sobolev" | "wahba-sobolev" => {
3016 (SphereMethod::Wahba, SphereWahbaKernel::Sobolev)
3017 }
3018 "pseudo" | "mgcv" | "sos" | "wahba_pseudo" | "wahba-pseudo" => {
3019 (SphereMethod::Wahba, SphereWahbaKernel::Pseudo)
3020 }
3021 "harmonic" | "spherical_harmonic" | "spherical-harmonic" => {
3022 (SphereMethod::Harmonic, SphereWahbaKernel::Sobolev)
3023 }
3024 other => {
3025 return Err(format!(
3026 "unsupported sphere kernel '{other}'; expected sobolev, pseudo, or harmonic"
3027 ));
3028 }
3029 };
3030 let wahba_kernel = match option_usize_any(options, &["lmax", "l_max", "l-max"]) {
3039 None => wahba_kernel,
3040 Some(_) if matches!(method, SphereMethod::Harmonic) => {
3041 return Err(
3042 "sphere smooth: lmax= states the truncation of a Wahba reproducing kernel \
3043 and does not apply to kernel=harmonic; use degree=/max_degree= to set the \
3044 harmonic degree"
3045 .to_string(),
3046 );
3047 }
3048 Some(lmax) => {
3049 if !(SPHERE_TRUNCATION_LMAX_RANGE).contains(&lmax) {
3050 return Err(format!(
3051 "sphere smooth: lmax={lmax} is out of range; the truncated Wahba \
3052 kernels support lmax in {}..={} (the device kernel bakes it in as a \
3053 compile-time bound)",
3054 SPHERE_TRUNCATION_LMAX_RANGE.start(),
3055 SPHERE_TRUNCATION_LMAX_RANGE.end()
3056 ));
3057 }
3058 let lmax = lmax as u16;
3059 match wahba_kernel {
3060 SphereWahbaKernel::Sobolev | SphereWahbaKernel::SobolevTruncated { .. } => {
3061 SphereWahbaKernel::SobolevTruncated { lmax }
3062 }
3063 SphereWahbaKernel::Pseudo | SphereWahbaKernel::PseudoTruncated { .. } => {
3064 SphereWahbaKernel::PseudoTruncated { lmax }
3065 }
3066 }
3067 }
3068 };
3069 let max_degree = if matches!(method, SphereMethod::Harmonic) {
3070 let degree =
3071 option_usize_any(options, &["degree", "l", "max_degree", "max-degree"])
3072 .or_else(|| option_usize(options, "centers"))
3073 .or_else(|| {
3074 option_usize_any(options, &["k", "basis_dim", "basis-dim", "basisdim"])
3075 .and_then(|k| (1..=128).find(|&l| l * (l + 2) >= k))
3076 })
3077 .unwrap_or_else(|| default_spherical_harmonic_degree(sizing_rows));
3078 if degree == 0 {
3079 return Err("sphere smooth requires degree/max_degree >= 1".to_string());
3080 }
3081 if degree > 32 {
3082 return Err(format!(
3083 "sphere smooth max_degree={} is too large for the dense harmonic engine (limit 32)",
3084 degree
3085 ));
3086 }
3087 Some(degree)
3088 } else {
3089 None
3090 };
3091 let penalty_order = option_usize(options, "penalty_order")
3092 .or_else(|| option_usize(options, "m"))
3093 .unwrap_or(DEFAULT_PENALTY_ORDER);
3094 let center_strategy = if matches!(method, SphereMethod::Wahba) {
3095 let mut centers = parse_countwith_basis_alias(
3096 options,
3097 "centers",
3098 default_num_centers(sizing_rows, cols.len()),
3099 )?;
3100 if penalty_order >= 4 {
3101 centers = centers.max(30);
3102 }
3103 CenterStrategy::FarthestPoint {
3104 num_centers: centers,
3105 }
3106 } else {
3107 CenterStrategy::FarthestPoint { num_centers: 0 }
3108 };
3109 Ok(SmoothBasisSpec::Sphere {
3110 feature_cols: cols.to_vec(),
3111 spec: SphericalSplineBasisSpec {
3112 center_strategy,
3113 penalty_order,
3114 double_penalty: smooth_double_penalty,
3115 radians,
3116 method,
3117 max_degree,
3118 wahba_kernel,
3119 identifiability: SphericalSplineIdentifiability::CenterSumToZero,
3120 },
3121 })
3122 }
3123 "curvature" => {
3124 validate_known_options("curvature", options, CURVATURE_SMOOTH_OPTION_KEYS)?;
3134 let kappa_opt = option_f64(options, "kappa");
3139 let kappa_fixed = kappa_opt.is_some();
3140 let kappa = kappa_opt.unwrap_or(0.0);
3141 if !kappa.is_finite() {
3142 return Err("curvature smooth requires a finite kappa".to_string());
3143 }
3144 let length_scale_opt = option_f64(options, "length_scale");
3151 let length_scale_fixed = length_scale_opt.is_some();
3152 let length_scale = length_scale_opt.unwrap_or(0.0);
3153 if !length_scale.is_finite() || length_scale < 0.0 {
3154 return Err(format!(
3155 "curvature smooth length_scale must be positive (or omitted for auto); got {length_scale}"
3156 ));
3157 }
3158 let centers = parse_countwith_basis_alias(
3159 options,
3160 "centers",
3161 default_num_centers(sizing_rows, cols.len()),
3162 )?;
3163 if centers < 2 {
3164 return Err("curvature smooth requires at least 2 centers".to_string());
3165 }
3166 let center_strategy = if has_explicit_countwith_basis_alias(options, "centers") {
3167 spatial_center_strategy_for_dimension(centers, cols.len())
3168 } else {
3169 auto_spatial_center_strategy(centers, cols.len())
3170 };
3171 Ok(SmoothBasisSpec::ConstantCurvature {
3172 feature_cols: cols.to_vec(),
3173 spec: ConstantCurvatureBasisSpec {
3174 center_strategy,
3175 kappa,
3176 kappa_fixed,
3177 length_scale,
3180 length_scale_fixed,
3181 double_penalty: option_bool(options, "double_penalty").unwrap_or(false),
3188 identifiability: ConstantCurvatureIdentifiability::CenterSumToZero,
3189 },
3190 })
3191 }
3192 "measurejet" => {
3193 validate_known_options("measurejet", options, MEASURE_JET_SMOOTH_OPTION_KEYS)?;
3199 let order_s = option_f64(options, "s").unwrap_or(0.0);
3200 if !(order_s.is_finite() && (order_s == 0.0 || (order_s > 0.0 && order_s < 2.0))) {
3203 return Err(format!(
3204 "measurejet smooth s must lie in (0, 2) (or be omitted for auto); got {order_s}"
3205 ));
3206 }
3207 let alpha =
3215 option_f64(options, "alpha").unwrap_or(MeasureJetBasisSpec::default().alpha);
3216 if !alpha.is_finite() {
3217 return Err("measurejet smooth requires a finite alpha".to_string());
3218 }
3219 let tau0 = option_f64(options, "tau").unwrap_or(1e-3);
3220 if !(tau0.is_finite() && tau0 >= 0.0) {
3221 return Err(format!(
3222 "measurejet smooth tau must be finite and nonnegative; got {tau0}"
3223 ));
3224 }
3225 let num_scales = option_usize(options, "scales").unwrap_or(0);
3226 let length_scale = option_f64(options, "length_scale").unwrap_or(0.0);
3227 if !length_scale.is_finite() || length_scale < 0.0 {
3228 return Err(format!(
3229 "measurejet smooth length_scale must be positive (or omitted for auto); got {length_scale}"
3230 ));
3231 }
3232 let centers = parse_countwith_basis_alias(
3233 options,
3234 "centers",
3235 default_num_centers(sizing_rows, cols.len()),
3236 )?;
3237 if centers < 3 {
3238 return Err("measurejet smooth requires at least 3 centers".to_string());
3239 }
3240 let center_strategy = if has_explicit_countwith_basis_alias(options, "centers") {
3241 spatial_center_strategy_for_dimension(centers, cols.len())
3242 } else {
3243 auto_spatial_center_strategy(centers, cols.len())
3244 };
3245 let multiscale = option_bool(options, "multiscale").unwrap_or(false);
3249 let learn_length_scale =
3261 option_bool(options, "learn_length_scale").unwrap_or(length_scale == 0.0);
3262 Ok(SmoothBasisSpec::MeasureJet {
3263 feature_cols: cols.to_vec(),
3264 spec: MeasureJetBasisSpec {
3265 center_strategy,
3266 order_s,
3267 alpha,
3268 tau0,
3269 num_scales,
3270 length_scale,
3273 double_penalty: smooth_double_penalty,
3274 learn_length_scale,
3275 multiscale,
3276 identifiability: MeasureJetIdentifiability::CenterSumToZero,
3277 frozen_quadrature: None,
3278 },
3279 input_scale: None,
3280 })
3281 }
3282 "matern" => {
3283 validate_known_options("matern", options, MATERN_SMOOTH_OPTION_KEYS)?;
3288 let plan = plan_spatial_basis(
3289 sizing_rows,
3290 cols.len(),
3291 CenterCountRequest::Default,
3292 DuchonNullspaceOrder::Zero,
3293 option_bool(options, "scale_dims").unwrap_or(false),
3294 policy,
3295 )
3296 .map_err(|e| e.to_string())?;
3297 let univariate_floor = if cols.len() == 1 {
3300 heuristic_knots_for_column(ds.values.column(cols[0]))
3301 .saturating_add(DEFAULT_BSPLINE_DEGREE + 1)
3302 } else {
3303 0
3304 };
3305 let centers = parse_countwith_basis_alias(
3306 options,
3307 "centers",
3308 cap_default_spatial_centers(
3309 options,
3310 default_matern_center_count(
3311 sizing_rows,
3312 cols.len(),
3313 plan.centers,
3314 univariate_floor,
3315 ),
3316 ),
3317 )?;
3318 let center_strategy = if has_explicit_countwith_basis_alias(options, "centers") {
3319 spatial_center_strategy_for_dimension(centers, cols.len())
3320 } else {
3321 auto_spatial_center_strategy(centers, cols.len())
3322 };
3323 let nu = parse_matern_nu(options.get("nu").map(String::as_str).unwrap_or("5/2"))?;
3324 if matches!(nu, MaternNu::Half) && cols.len() >= 2 {
3330 return Err(TermBuilderError::unsupported_feature(format!(
3331 "matern() with nu=1/2 is not supported for d>=2 (got {} covariates): \
3332 the exponential kernel's Laplacian is singular at center collisions, \
3333 which makes the operator-collocation penalty non-invertible. \
3334 Choose nu>=3/2 (e.g. nu=3/2 or the default nu=5/2) for multi-dimensional smooths.",
3335 cols.len()
3336 ))
3337 .to_string());
3338 }
3339 let aniso_log_scales = if option_bool(options, "scale_dims").unwrap_or(false) {
3340 Some(vec![0.0; cols.len()])
3341 } else {
3342 None
3343 };
3344 let periodic = parse_periodic_axes_option(options, cols.len())?;
3345 reject_unconsumable_radial_period_declaration(
3346 "matern",
3347 options,
3348 cols.len(),
3349 periodic.as_deref(),
3350 false,
3351 )?;
3352 Ok(SmoothBasisSpec::Matern {
3353 feature_cols: cols.to_vec(),
3354 spec: MaternBasisSpec {
3355 center_strategy,
3356 periodic,
3357 length_scale: option_f64(options, "length_scale")
3373 .map(MaternLengthScale::fixed)
3374 .unwrap_or_else(MaternLengthScale::auto),
3375 nu,
3376 include_intercept: option_bool(options, "include_intercept").unwrap_or(false),
3377 double_penalty: smooth_double_penalty,
3378 identifiability: parse_matern_identifiability(options)
3379 .map_err(|e| e.to_string())?,
3380 aniso_log_scales,
3381 },
3386 input_scale: None,
3387 })
3388 }
3389 "duchon" | "ds" => {
3390 validate_known_options("duchon", options, DUCHON_SMOOTH_OPTION_KEYS)?;
3391 if options.contains_key("double_penalty") {
3392 return Err(TermBuilderError::incompatible_config(format!(
3393 "Duchon smooth '{}' does not support double_penalty; the Duchon smoother already ships its native reproducing-norm penalty plus a null-space shrinkage ridge.",
3394 vars.join(", ")
3395 ))
3396 .to_string());
3397 }
3398 let requested_nullspace_order = parse_duchon_order_opt(options)?;
3399 let length_scale = option_f64_strict(options, "length_scale")?;
3400 let (nullspace_order, power) = match parse_duchon_power_policy(options)? {
3413 DuchonPowerPolicy::Explicit(req_power) => {
3414 if length_scale.is_some() && req_power.fract() != 0.0 {
3415 return Err(TermBuilderError::incompatible_config(format!(
3416 "hybrid Duchon-Matern smooth '{}' (length_scale=...) requires an integer power, got power={}; \
3417 drop length_scale to use the scale-free structural kernel with a fractional power.",
3418 vars.join(", "),
3419 req_power,
3420 ))
3421 .to_string());
3422 }
3423 (
3424 requested_nullspace_order.unwrap_or(DuchonNullspaceOrder::Linear),
3425 req_power,
3426 )
3427 }
3428 DuchonPowerPolicy::CubicStructuralDefault => {
3429 match length_scale {
3444 None => {
3445 let (default_order, s) =
3446 crate::basis::duchon_cubic_default(cols.len());
3447 (requested_nullspace_order.unwrap_or(default_order), s)
3448 }
3449 Some(_) => {
3450 let (default_order, s_frac) =
3476 crate::basis::duchon_cubic_default(cols.len());
3477 (
3478 requested_nullspace_order.unwrap_or(default_order),
3479 s_frac.floor(),
3480 )
3481 }
3482 }
3483 }
3484 };
3485 let plan = plan_spatial_basis(
3486 sizing_rows,
3487 cols.len(),
3488 CenterCountRequest::Default,
3489 nullspace_order,
3490 option_bool(options, "scale_dims").unwrap_or(false),
3491 policy,
3492 )
3493 .map_err(|e| e.to_string())?;
3494 let centers_explicit = has_explicit_countwith_basis_alias(options, "centers");
3495 let polynomial_cols = match nullspace_order {
3496 DuchonNullspaceOrder::Zero => 1,
3497 DuchonNullspaceOrder::Linear => cols.len() + 1,
3498 DuchonNullspaceOrder::Degree(degree) => {
3499 crate::basis::duchon_nullspace_dimension(cols.len(), degree)
3500 }
3501 };
3502 let univariate_floor = if cols.len() == 1 {
3505 heuristic_knots_for_column(ds.values.column(cols[0]))
3506 .saturating_add(DEFAULT_BSPLINE_DEGREE + 1)
3507 } else {
3508 0
3509 };
3510 let default_centers = default_duchon_center_count(
3511 sizing_rows,
3512 cols.len(),
3513 plan.centers,
3514 polynomial_cols,
3515 univariate_floor,
3516 );
3517 let spectral_rank = option_usize(options, "rank");
3518 let center_default = if spectral_rank.is_some() {
3519 count_unique_coordinate_rows(ds.values.view(), &cols).min(2000)
3531 } else {
3532 cap_default_spatial_centers(options, default_centers)
3533 };
3534 let requested_centers =
3535 parse_countwith_basis_alias(options, "centers", center_default)?;
3536 if requested_centers > ds.values.nrows() {
3537 return Err(TermBuilderError::incompatible_config(format!(
3538 "Duchon smooth '{}' requested {requested_centers} centers but only {} rows are available",
3539 vars.join(", "),
3540 ds.values.nrows(),
3541 ))
3542 .to_string());
3543 }
3544 if requested_centers <= polynomial_cols {
3545 return Err(TermBuilderError::incompatible_config(format!(
3546 "Duchon smooth '{}' requested basis dimension {} but order={:?} in {}D needs {} polynomial null-space columns; choose centers/k > {}",
3547 vars.join(", "),
3548 requested_centers,
3549 nullspace_order,
3550 cols.len(),
3551 polynomial_cols,
3552 polynomial_cols,
3553 ))
3554 .to_string());
3555 }
3556 if let Some(rank) = spectral_rank
3557 && (rank <= polynomial_cols || rank > requested_centers)
3558 {
3559 return Err(TermBuilderError::incompatible_config(format!(
3560 "Duchon smooth '{}' spectral rank must satisfy {} < rank <= centers (got rank={rank}, centers={requested_centers})",
3561 vars.join(", "),
3562 polynomial_cols,
3563 ))
3564 .to_string());
3565 }
3566 let mut centers = requested_centers;
3567 if !centers_explicit && ds.values.nrows() <= 32 && smooth_coordinate_count >= 5 {
3568 centers = centers.max(polynomial_cols + 4);
3569 }
3570 let aniso_log_scales = if option_bool(options, "scale_dims").unwrap_or(false) {
3571 Some(vec![0.0; cols.len()])
3572 } else {
3573 None
3574 };
3575 let operator_penalties = DuchonOperatorPenaltySpec::all_disabled();
3584 let mut periodic = parse_periodic_axes_option(options, cols.len())?;
3592 if cols.len() == 1
3593 && let Some(axes) = periodic.as_mut()
3594 && axes.len() == 1
3595 && axes[0].is_none()
3596 {
3597 let (minv, maxv) = col_minmax(ds.values.column(cols[0]))?;
3598 if maxv > minv {
3599 axes[0] = Some(maxv - minv);
3600 }
3601 }
3602 let boundary = if cols.len() == 1 {
3603 let c = cols[0];
3604 let (minv, maxv) = col_minmax(ds.values.column(c))?;
3605 parse_cyclic_boundary(options, minv, maxv)?
3606 } else {
3607 OneDimensionalBoundary::Open
3608 };
3609 let is_periodic = periodic
3610 .as_ref()
3611 .is_some_and(|axes| axes.iter().any(Option::is_some))
3612 || matches!(boundary, OneDimensionalBoundary::Cyclic { .. });
3613 reject_unconsumable_radial_period_declaration(
3614 "duchon",
3615 options,
3616 cols.len(),
3617 periodic.as_deref(),
3618 matches!(boundary, OneDimensionalBoundary::Cyclic { .. }),
3619 )?;
3620 if spectral_rank.is_some() && is_periodic {
3621 return Err(TermBuilderError::incompatible_config(
3622 "Duchon spectral rank is defined for the scale-free open-domain kernel, \
3623 not a periodic image expansion"
3624 .to_string(),
3625 )
3626 .to_string());
3627 }
3628 let center_strategy = if spectral_rank.is_some() {
3629 let mut coordinates = Array2::<f64>::zeros((ds.values.nrows(), cols.len()));
3637 for (axis, &column) in cols.iter().enumerate() {
3638 coordinates
3639 .column_mut(axis)
3640 .assign(&ds.values.column(column));
3641 }
3642 let sampled = select_r_uniform_subsample_centers(coordinates.view(), centers, 1)
3643 .map_err(|error| error.to_string())?;
3644 CenterStrategy::UserProvided(sampled)
3645 } else if is_periodic {
3646 if centers_explicit {
3647 spatial_center_strategy_for_dimension(centers, cols.len())
3648 } else {
3649 auto_spatial_center_strategy(centers, cols.len())
3650 }
3651 } else {
3652 duchon_center_strategy(centers, cols.len(), !centers_explicit)
3653 };
3654 let center_strategy = match spectral_rank {
3655 Some(rank) => CenterStrategy::DuchonSpectral {
3656 knots: Box::new(center_strategy),
3657 basis: DuchonSpectralBasis::Fresh { rank },
3658 },
3659 None => center_strategy,
3660 };
3661 Ok(SmoothBasisSpec::Duchon {
3662 feature_cols: cols.to_vec(),
3663 spec: DuchonBasisSpec {
3664 center_strategy,
3665 periodic,
3666 length_scale,
3667 power,
3668 nullspace_order,
3669 identifiability: parse_spatial_identifiability(options)
3670 .map_err(|e| e.to_string())?,
3671 aniso_log_scales,
3672 operator_penalties,
3673 boundary,
3674 radial_reparam: None,
3675 },
3676 input_scale: None,
3677 })
3678 }
3679 "tensor" | "te" | "ti" | "t2" => {
3680 validate_known_options("tensor", options, TENSOR_SMOOTH_OPTION_KEYS)?;
3681 if cols.len() < 2 {
3682 return Err(TermBuilderError::incompatible_config(format!(
3683 "tensor smooth expects at least 2 variables, got {}",
3684 cols.len()
3685 ))
3686 .to_string());
3687 }
3688 let dim = cols.len();
3689
3690 if let Some(raw) = options.get("bs").or_else(|| options.get("type"))
3713 && bs_selector_is_vector(raw)
3714 {
3715 let per_margin = parse_option_list(raw);
3716 if per_margin.len() != dim {
3717 return Err(TermBuilderError::invalid_option(format!(
3718 "tensor smooth per-margin bs vector has {} entries but the smooth has {} margins",
3719 per_margin.len(),
3720 dim
3721 ))
3722 .to_string());
3723 }
3724 for (axis, margin_bs) in per_margin.iter().enumerate() {
3725 if !tensor_margin_bs_is_supported(margin_bs) {
3726 return Err(TermBuilderError::unsupported_feature(format!(
3727 "tensor smooth margin {axis} basis '{margin_bs}' is not a supported penalized-spline margin; \
3728 tensor margins accept tp/tps/ps/bs/cr/cc"
3729 ))
3730 .to_string());
3731 }
3732 }
3733 }
3734 validate_tensor_boundary_tokens(options, dim)?;
3738 let periodic_axes = parse_tensor_periodic_axes(options, dim)?;
3739 reject_unconsumable_period_declaration("tensor", options, &periodic_axes)?;
3740 if let Some(key) = PERIOD_ENDPOINT_OPTION_KEYS
3744 .iter()
3745 .find(|key| options.contains_key(**key))
3746 {
3747 return Err(TermBuilderError::invalid_option(format!(
3748 "tensor(): `{key}=` declares one axis's periodic domain and has no per-margin \
3749 form; on a tensor smooth give periods=[...] (with origins=[...] for the \
3750 domain start), which name their margin"
3751 ))
3752 .to_string());
3753 }
3754 let periods_opt = parse_periods(options, &periodic_axes)?;
3755 let origins_opt = parse_period_origins(options, &periodic_axes)?;
3756 let requested_degrees = parse_tensor_per_axis_usize(options, "degree", dim)?;
3761 let requested_penalty_orders =
3762 parse_tensor_per_axis_usize(options, "penalty_order", dim)?;
3763 let axis_degree = |axis: usize| -> usize {
3764 requested_degrees[axis].unwrap_or(DEFAULT_BSPLINE_DEGREE)
3765 };
3766 let axis_penalty_order = |axis: usize| -> usize {
3767 requested_penalty_orders[axis]
3768 .unwrap_or(if axis_degree(axis) > 1 { 2 } else { 1 })
3769 };
3770 let (mut k_list, k_inferred) = parse_tensor_k_list(options, cols, ds)?;
3771 if ds.values.nrows() <= 32 && smooth_coordinate_count >= 5 {
3772 for (axis, k) in k_list.iter_mut().enumerate() {
3773 *k = (*k).min(axis_degree(axis) + 2);
3774 }
3775 }
3776 if k_inferred {
3777 inference_notes.push(format!(
3778 "Automatically set per-margin basis sizes {:?} for tensor smooth '{}' \
3779 (dimension-aware tensor budget: total ∏k kept near the mgcv-te default \
3780 and within the data support, distributed geometrically across margins and \
3781 capped per margin by each column's resolution). \
3782 Override with k=<int> or k=[k0,k1,...].",
3783 k_list,
3784 vars.join(",")
3785 ));
3786 }
3787 let per_axis_bs: Vec<Option<String>> =
3800 match options.get("bs").or_else(|| options.get("type")) {
3801 Some(raw) if bs_selector_is_vector(raw) => {
3802 let list = parse_option_list(raw);
3803 (0..dim).map(|a| list.get(a).cloned()).collect()
3804 }
3805 Some(raw) => {
3806 let scalar = raw
3807 .trim()
3808 .trim_matches('"')
3809 .trim_matches('\'')
3810 .to_ascii_lowercase();
3811 vec![Some(scalar); dim]
3812 }
3813 None => vec![None; dim],
3814 };
3815 let margin_wants_cr = |bs: &Option<String>| -> bool {
3821 matches!(
3822 bs.as_deref(),
3823 None | Some("cr") | Some("cs") | Some("tp") | Some("tps")
3824 )
3825 };
3826 let requested_knot_placement = explicit_knot_placement(options)?;
3827 let mut margins: Vec<BSplineBasisSpec> = Vec::with_capacity(dim);
3828 let mut emitted_periods: Vec<Option<f64>> = Vec::with_capacity(dim);
3829 for axis in 0..dim {
3830 let c = cols[axis];
3831 let (data_min, data_max) = col_minmax(ds.values.column(c))?;
3832 let k_requested = k_list[axis];
3848 let n_distinct_axis = unique_count_column(ds.values.column(c));
3849 let k_axis = k_requested.min(n_distinct_axis).max(2);
3850 if k_axis < k_requested {
3851 log::info!(
3852 "tensor smooth: margin axis {axis} requested k={k_requested}, but the \
3853 covariate has only {n_distinct_axis} distinct value(s); reducing this \
3854 margin to k={k_axis} (mgcv-style data-support cap on the per-axis basis)."
3855 );
3856 }
3857 if k_axis < 2 {
3873 return Err(TermBuilderError::invalid_option(format!(
3874 "tensor smooth: k[{axis}]={k_axis} too small; tensor margins require k >= 2"
3875 ))
3876 .to_string());
3877 }
3878 let degree = axis_degree(axis);
3879 let penalty_order = axis_penalty_order(axis);
3880 let effective_degree = degree.min(k_axis - 1).max(1);
3881 let effective_penalty_order = penalty_order.min(effective_degree);
3882 let margin_is_cc = matches!(
3889 canonicalize_smooth_type(per_axis_bs[axis].as_deref().unwrap_or("")),
3890 "cc" | "cp" | "cyclic"
3891 );
3892 let (knotspec, boundary, axis_period) = if periodic_axes[axis] {
3893 let (domain_start, period_value) = match periods_opt[axis] {
3904 Some(period_value) => {
3905 if !period_value.is_finite() || period_value <= 0.0 {
3906 return Err(format!(
3907 "tensor smooth axis {axis}: period must be a positive finite value, got {period_value}"
3908 ));
3909 }
3910 (origins_opt[axis].unwrap_or(data_min), period_value)
3911 }
3912 None if margin_is_cc => {
3913 let span = data_max - data_min;
3914 if !span.is_finite() || span <= 0.0 {
3915 return Err(format!(
3916 "tensor smooth axis {axis}: cyclic margin requires a positive \
3917 observed data range to derive its period, got [{data_min}, {data_max}]"
3918 ));
3919 }
3920 (origins_opt[axis].unwrap_or(data_min), span)
3921 }
3922 None => {
3923 return Err(format!(
3924 "tensor smooth axis {axis} is periodic but requires an explicit \
3925 period: pass period=<value> (scalar) or period=[..., <value>, ...]. \
3926 Deriving the period from the observed data range is sample-dependent \
3927 (off-by-ε seam), so it is not inferred."
3928 ));
3929 }
3930 };
3931 let domain_end = domain_start + period_value;
3932 (
3933 BSplineKnotSpec::PeriodicUniform {
3934 data_range: (domain_start, domain_end),
3935 num_basis: k_axis,
3936 },
3937 OneDimensionalBoundary::Cyclic {
3938 start: domain_start,
3939 end: domain_end,
3940 },
3941 Some(period_value),
3942 )
3943 } else if margin_wants_cr(&per_axis_bs[axis])
3944 && requested_knot_placement.is_none()
3945 && requested_degrees[axis].is_none_or(|d| d == CR_MARGIN_DEGREE)
3946 && requested_penalty_orders[axis]
3947 .is_none_or(|m| m == CR_MARGIN_PENALTY_ORDER)
3948 && k_axis >= 3
3949 {
3950 let cr_knots = crate::basis::select_cr_knots(ds.values.column(c), k_axis)
3964 .map_err(|e| e.to_string())?;
3965 (
3966 BSplineKnotSpec::NaturalCubicRegression { knots: cr_knots },
3967 OneDimensionalBoundary::Open,
3968 None,
3969 )
3970 } else {
3971 let num_internal_knots = if effective_degree < degree {
3978 k_axis.saturating_sub(effective_degree + 1)
3979 } else {
3980 k_axis.saturating_sub(degree + 1).max(1)
3981 };
3982 let knotspec = match requested_knot_placement
3983 .unwrap_or(crate::basis::BSplineKnotPlacement::Uniform)
3984 {
3985 crate::basis::BSplineKnotPlacement::Uniform => BSplineKnotSpec::Generate {
3986 data_range: (data_min, data_max),
3987 num_internal_knots,
3988 },
3989 crate::basis::BSplineKnotPlacement::Quantile => {
3990 crate::basis::auto_knot_vector_1d_quantile(
3991 ds.values.column(c),
3992 num_internal_knots,
3993 effective_degree,
3994 )
3995 .map_err(|e| e.to_string())?;
3996 BSplineKnotSpec::Automatic {
3997 num_internal_knots: Some(num_internal_knots),
3998 placement: crate::basis::BSplineKnotPlacement::Quantile,
3999 }
4000 }
4001 };
4002 (knotspec, OneDimensionalBoundary::Open, None)
4003 };
4004 margins.push(BSplineBasisSpec {
4009 degree: effective_degree,
4010 penalty_order: effective_penalty_order,
4011 knotspec,
4012 double_penalty: false,
4013 identifiability: BSplineIdentifiability::None,
4014 boundary,
4015 boundary_conditions: BSplineBoundaryConditions::default(),
4016 });
4017 emitted_periods.push(axis_period);
4018 }
4019 let canon_cols: Vec<usize> = {
4040 let mut perm: Vec<usize> = (0..dim).collect();
4041 perm.sort_by_key(|&a| cols[a]);
4042 if perm.iter().enumerate().any(|(i, &a)| i != a) {
4043 margins = perm.iter().map(|&a| margins[a].clone()).collect();
4044 emitted_periods = perm.iter().map(|&a| emitted_periods[a]).collect();
4045 }
4046 perm.iter().map(|&a| cols[a]).collect()
4047 };
4048 let any_periodic = emitted_periods.iter().any(|p| p.is_some());
4049 let periods_vec = if any_periodic {
4050 emitted_periods
4051 } else {
4052 Vec::new()
4053 };
4054 let tensor_double_penalty = smooth_double_penalty;
4058 Ok(SmoothBasisSpec::TensorBSpline {
4059 feature_cols: canon_cols,
4060 spec: TensorBSplineSpec {
4061 marginalspecs: margins,
4062 periods: periods_vec,
4063 double_penalty: tensor_double_penalty,
4064 identifiability: parse_tensor_identifiability(options, kind)?,
4065 penalty_decomposition: if matches!(kind, SmoothKind::T2)
4075 || type_opt.as_str() == "t2"
4076 {
4077 TensorBSplinePenaltyDecomposition::Separable
4078 } else {
4079 TensorBSplinePenaltyDecomposition::MarginalKroneckerSum
4080 },
4081 },
4082 })
4083 }
4084 "pca" => {
4085 validate_known_options("pca", options, PCA_SMOOTH_OPTION_KEYS)?;
4086 let path = options
4087 .get("lazy_path")
4088 .or_else(|| options.get("pca_basis_path"))
4089 .or_else(|| options.get("path"))
4090 .map(|raw| PathBuf::from(strip_quotes(raw)));
4091 let Some(path) = path else {
4092 return Err(TermBuilderError::incompatible_config(
4093 "pca smooth requires lazy_path=... on the formula path",
4094 )
4095 .to_string());
4096 };
4097 let k = option_usize_any(options, &["k", "basis_dim", "basis-dim", "basisdim"])
4098 .unwrap_or(0);
4099 let chunk_size = option_usize(options, "chunk_size").unwrap_or(DEFAULT_PCA_CHUNK_SIZE);
4100 Ok(SmoothBasisSpec::Pca {
4101 feature_cols: cols.to_vec(),
4102 basis_matrix: Array2::<f64>::zeros((cols.len(), k)),
4103 centered: option_bool(options, "centered").unwrap_or(true),
4104 smooth_penalty: option_f64(options, "smooth_penalty").unwrap_or(1.0),
4105 center_mean: None,
4106 pca_basis_path: Some(path),
4107 chunk_size,
4108 })
4109 }
4110 other => Err(TermBuilderError::unsupported_feature(format!(
4111 "unsupported smooth type '{other}'"
4112 ))
4113 .to_string()),
4114 }
4115}
4116
4117pub fn enable_scale_dimensions(spec: &mut TermCollectionSpec) {
4119 for smooth in spec.smooth_terms.iter_mut() {
4120 promote_thin_plate_for_scale_dimensions(&mut smooth.basis);
4127 match &mut smooth.basis {
4128 SmoothBasisSpec::Matern {
4129 feature_cols,
4130 spec: matern,
4131 ..
4132 } => {
4133 if matern.aniso_log_scales.is_none() {
4134 let d = feature_cols.len();
4135 matern.aniso_log_scales = Some(vec![0.0; d]);
4136 }
4137 }
4138 SmoothBasisSpec::Duchon {
4139 feature_cols,
4140 spec: duchon,
4141 ..
4142 } => {
4143 if duchon.aniso_log_scales.is_none() {
4144 let d = feature_cols.len();
4145 duchon.aniso_log_scales = Some(vec![0.0; d]);
4146 }
4147 }
4148 SmoothBasisSpec::ByVariable { .. }
4153 | SmoothBasisSpec::FactorSumToZero { .. }
4154 | SmoothBasisSpec::BSpline1D { .. }
4155 | SmoothBasisSpec::BySmooth { .. }
4156 | SmoothBasisSpec::FactorSmooth { .. }
4157 | SmoothBasisSpec::ThinPlate { .. }
4158 | SmoothBasisSpec::Sphere { .. }
4159 | SmoothBasisSpec::ConstantCurvature { .. }
4160 | SmoothBasisSpec::MeasureJet { .. }
4161 | SmoothBasisSpec::Pca { .. }
4162 | SmoothBasisSpec::TensorBSpline { .. } => {}
4163 }
4164 }
4165}
4166
4167fn promote_thin_plate_for_scale_dimensions(basis: &mut SmoothBasisSpec) {
4202 let SmoothBasisSpec::ThinPlate {
4203 feature_cols,
4204 spec,
4205 input_scale,
4206 } = &*basis
4207 else {
4208 return;
4209 };
4210 let d = feature_cols.len();
4211 if d <= 1 {
4212 return;
4213 }
4214 let m = thin_plate_penalty_order(d);
4219 let nullspace_order = match m {
4220 0 | 1 => DuchonNullspaceOrder::Zero,
4221 2 => DuchonNullspaceOrder::Linear,
4222 _ => DuchonNullspaceOrder::Degree(m - 1),
4223 };
4224 let duchon_spec = DuchonBasisSpec {
4225 center_strategy: spec.center_strategy.clone(),
4226 periodic: spec.periodic.clone(),
4227 length_scale: None,
4232 power: 0.0,
4234 nullspace_order,
4235 identifiability: spec.identifiability.clone(),
4236 aniso_log_scales: Some(vec![0.0; d]),
4240 operator_penalties: DuchonOperatorPenaltySpec::default(),
4241 boundary: OneDimensionalBoundary::Open,
4242 radial_reparam: None,
4243 };
4244 let feature_cols = feature_cols.clone();
4245 let input_scale = *input_scale;
4246 *basis = SmoothBasisSpec::Duchon {
4249 feature_cols,
4250 spec: duchon_spec,
4251 input_scale,
4252 };
4253}
4254
4255pub fn spatial_center_strategy_for_dimension(num_centers: usize, d: usize) -> CenterStrategy {
4260 if d <= 3 {
4261 CenterStrategy::FarthestPoint { num_centers }
4268 } else {
4269 default_spatial_center_strategy(num_centers, d)
4270 }
4271}
4272
4273fn duchon_center_strategy(num_centers: usize, d: usize, automatic: bool) -> CenterStrategy {
4289 let realized = if d == 1 {
4290 CenterStrategy::UniformGrid {
4291 points_per_dim: num_centers,
4292 }
4293 } else {
4294 spatial_center_strategy_for_dimension(num_centers, d)
4295 };
4296 if automatic {
4297 CenterStrategy::Auto(Box::new(realized))
4298 } else {
4299 realized
4300 }
4301}
4302
4303pub fn col_minmax(col: ArrayView1<'_, f64>) -> Result<(f64, f64), String> {
4304 let min = col.iter().fold(f64::INFINITY, |a, &b| a.min(b));
4305 let max = col.iter().fold(f64::NEG_INFINITY, |a, &b| a.max(b));
4306 if !min.is_finite() || !max.is_finite() {
4307 return Err(TermBuilderError::degenerate_data(
4308 "non-finite data encountered while inferring knot range",
4309 )
4310 .to_string());
4311 }
4312 if (max - min).abs() < 1e-12 {
4313 Ok((min, min + 1e-6))
4314 } else {
4315 Ok((min, max))
4316 }
4317}
4318
4319pub fn unique_count_column(col: ArrayView1<'_, f64>) -> usize {
4320 use std::collections::HashSet;
4321 let mut set = HashSet::<u64>::with_capacity(col.len());
4322 for &v in col {
4323 set.insert(gam_data::canonical_level_bits(v));
4324 }
4325 set.len().max(1)
4326}
4327
4328pub(crate) const CR_MIN_KNOTS: usize = 3;
4334
4335fn capped_cr_marginal_knotspec(
4362 col: ArrayView1<'_, f64>,
4363 k_cr_requested: usize,
4364 label: &str,
4365 inference_notes: &mut Vec<String>,
4366) -> Result<Option<BSplineKnotSpec>, String> {
4367 let n_distinct = unique_count_column(col);
4368 let k_cr = k_cr_requested.min(n_distinct);
4369 if k_cr < CR_MIN_KNOTS {
4370 inference_notes.push(format!(
4371 "Smooth '{label}': cubic-regression ('cr'/'cs'/'sz') basis requested k={k_cr_requested}, \
4372 but the covariate has only {n_distinct} distinct value(s) — too few to support a cubic \
4373 regression spline (needs >= {CR_MIN_KNOTS} distinct values). Degraded to the linear \
4374 B-spline marginal the default basis builds on the same data."
4375 ));
4376 return Ok(None);
4377 }
4378 if k_cr < k_cr_requested {
4379 inference_notes.push(format!(
4380 "Smooth '{label}': cubic-regression ('cr'/'cs'/'sz') basis reduced from k={k_cr_requested} \
4381 to k={k_cr} to match the covariate's {n_distinct} distinct value(s) (mgcv-style \
4382 data-support cap; a cr basis cannot place more value-knots than the data has)."
4383 ));
4384 }
4385 let cr_knots = crate::basis::select_cr_knots(col, k_cr).map_err(|e| e.to_string())?;
4386 Ok(Some(BSplineKnotSpec::NaturalCubicRegression {
4387 knots: cr_knots,
4388 }))
4389}
4390
4391fn min_per_group_unique_count(
4398 feature_col: ArrayView1<'_, f64>,
4399 group_col: ArrayView1<'_, f64>,
4400) -> usize {
4401 use std::collections::{HashMap, HashSet};
4402 let mut per_group: HashMap<u64, HashSet<u64>> = HashMap::new();
4403 for (xi, gi) in feature_col.iter().zip(group_col.iter()) {
4404 per_group
4405 .entry(gam_data::canonical_level_bits(*gi))
4406 .or_default()
4407 .insert(gam_data::canonical_level_bits(*xi));
4408 }
4409 per_group
4410 .values()
4411 .map(|s| s.len())
4412 .min()
4413 .unwrap_or(1)
4414 .max(1)
4415}
4416
4417pub fn heuristic_knots_for_column(col: ArrayView1<'_, f64>) -> usize {
4445 const MAX_DEFAULT_INTERNAL_KNOTS: usize = 8;
4448 let unique = unique_count_column(col);
4449 (unique / 4).clamp(4, MAX_DEFAULT_INTERNAL_KNOTS)
4450}
4451
4452fn heuristic_tensor_margin_knots(cols: &[usize], ds: &Dataset) -> Vec<usize> {
4473 let d = cols.len().max(1);
4474 let degree = DEFAULT_BSPLINE_DEGREE;
4475 let min_k = degree + 2; let n = ds.values.nrows();
4477
4478 let per_margin_cap: Vec<usize> = cols
4482 .iter()
4483 .map(|&c| heuristic_knots_for_column(ds.values.column(c)).max(min_k))
4484 .collect();
4485
4486 let mgcv_like_per_margin = match d {
4493 2 => 7usize,
4494 3 => 5usize,
4495 _ => 4usize,
4496 };
4497 let mgcv_like_total = (mgcv_like_per_margin as f64).powi(d as i32);
4498 let data_budget = (n as f64) * 0.8;
4499 let p_target = mgcv_like_total
4500 .max(min_k.pow(d as u32) as f64)
4501 .min(data_budget);
4502
4503 let geo_per_margin = p_target.powf(1.0 / d as f64).round() as usize;
4506 let unclamped: Vec<usize> = per_margin_cap
4507 .iter()
4508 .map(|&cap| geo_per_margin.clamp(min_k, cap))
4509 .collect();
4510
4511 let mut k_list = unclamped;
4516 loop {
4517 let product: f64 = k_list.iter().map(|&k| k as f64).product();
4518 if product >= p_target {
4519 break;
4520 }
4521 let Some(idx) = k_list
4524 .iter()
4525 .zip(per_margin_cap.iter())
4526 .enumerate()
4527 .filter(|&(_, (k, cap))| k < cap)
4528 .max_by_key(|&(_, (k, cap))| (cap - k, *cap))
4529 .map(|(i, _)| i)
4530 else {
4531 break;
4532 };
4533 k_list[idx] += 1;
4534 }
4535 k_list
4536}
4537
4538pub fn heuristic_centers(n: usize, d: usize) -> usize {
4539 default_num_centers(n, d)
4540}
4541
4542fn parse_endpoint_side(
4547 value: &str,
4548 context: &str,
4549) -> Result<BSplineEndpointBoundaryCondition, String> {
4550 match value.trim().to_ascii_lowercase().as_str() {
4551 "" | "none" | "open" | "unconstrained" | "free" => {
4552 Ok(BSplineEndpointBoundaryCondition::Free)
4553 }
4554 "clamped" | "clamp" | "zero_derivative" | "zero-derivative" => {
4555 Ok(BSplineEndpointBoundaryCondition::Clamped)
4556 }
4557 "anchored" | "anchor" | "zero" | "zero_value" | "zero-value" => {
4558 Ok(BSplineEndpointBoundaryCondition::Anchored { value: 0.0 })
4559 }
4560 other => Err(format!(
4561 "unsupported {context} boundary condition '{other}'; expected free, clamped, or anchored"
4562 )),
4563 }
4564}
4565
4566fn boundary_anchor_value(
4567 options: &BTreeMap<String, String>,
4568 side: &str,
4569 fallback: Option<f64>,
4570) -> Option<f64> {
4571 [
4572 format!("anchor_{side}"),
4573 format!("{side}_anchor"),
4574 format!("anchor-value-{side}"),
4575 ]
4576 .iter()
4577 .find_map(|key| option_f64(options, key))
4578 .or(fallback)
4579}
4580
4581fn apply_anchor_value(
4582 cond: BSplineEndpointBoundaryCondition,
4583 value: Option<f64>,
4584) -> BSplineEndpointBoundaryCondition {
4585 match cond {
4586 BSplineEndpointBoundaryCondition::Anchored { .. } => {
4587 BSplineEndpointBoundaryCondition::Anchored {
4588 value: value.unwrap_or(0.0),
4589 }
4590 }
4591 other => other,
4592 }
4593}
4594
4595fn parse_bspline_boundary_conditions(
4596 options: &BTreeMap<String, String>,
4597) -> Result<BSplineBoundaryConditions, String> {
4598 let fallback_anchor = option_f64(options, "anchor")
4599 .or_else(|| option_f64(options, "anchor_value"))
4600 .or_else(|| option_f64(options, "value"));
4601 let global_boundary_conditions = options
4608 .get("boundary_conditions")
4609 .or_else(|| options.get("bc"))
4610 .or_else(|| options.get("boundary"));
4611 let mut boundary_conditions = BSplineBoundaryConditions::default();
4612
4613 if let Some(raw_boundary_conditions) = global_boundary_conditions {
4614 let cond = parse_endpoint_side(raw_boundary_conditions, "boundary_conditions")?;
4615 let side = options
4616 .get("side")
4617 .map(|s| s.trim().to_ascii_lowercase())
4618 .unwrap_or_else(|| "both".to_string());
4619 match side.as_str() {
4620 "both" | "all" | "endpoints" => {
4621 boundary_conditions.left = cond;
4622 boundary_conditions.right = cond;
4623 }
4624 "left" | "start" | "lower" => boundary_conditions.left = cond,
4625 "right" | "end" | "upper" => boundary_conditions.right = cond,
4626 other => {
4627 return Err(format!(
4628 "unsupported B-spline boundary side '{other}'; expected left, right, or both"
4629 ));
4630 }
4631 }
4632 }
4633
4634 if let Some(raw) = options
4635 .get("bc_left")
4636 .or_else(|| options.get("left_bc"))
4637 .or_else(|| options.get("bc_start"))
4638 .or_else(|| options.get("start_bc"))
4639 {
4640 boundary_conditions.left = parse_endpoint_side(raw, "left endpoint")?;
4641 }
4642 if let Some(raw) = options
4643 .get("bc_right")
4644 .or_else(|| options.get("right_bc"))
4645 .or_else(|| options.get("bc_end"))
4646 .or_else(|| options.get("end_bc"))
4647 {
4648 boundary_conditions.right = parse_endpoint_side(raw, "right endpoint")?;
4649 }
4650
4651 boundary_conditions.left = apply_anchor_value(
4652 boundary_conditions.left,
4653 boundary_anchor_value(options, "left", fallback_anchor),
4654 );
4655 boundary_conditions.right = apply_anchor_value(
4656 boundary_conditions.right,
4657 boundary_anchor_value(options, "right", fallback_anchor),
4658 );
4659
4660 if options.contains_key("side") && global_boundary_conditions.is_none() {
4666 return Err(TermBuilderError::invalid_option(
4667 "`side=` selects which endpoint a boundary condition applies to, but this smooth declares none; add bc=<condition> or drop it",
4668 )
4669 .to_string());
4670 }
4671 if !boundary_conditions.has_anchor()
4672 && let Some(key) = ANCHOR_VALUE_OPTION_KEYS
4673 .iter()
4674 .find(|key| options.contains_key(**key))
4675 {
4676 return Err(TermBuilderError::invalid_option(format!(
4677 "`{key}=` sets the value an ANCHORED endpoint is pinned to, but no endpoint of this smooth is anchored; add bc=anchored (or bc_left=/bc_right=anchored) or drop it"
4678 ))
4679 .to_string());
4680 }
4681
4682 Ok(boundary_conditions)
4683}
4684
4685const ANCHOR_VALUE_OPTION_KEYS: [&str; 8] = [
4688 "anchor",
4689 "anchor_value",
4690 "value",
4691 "anchor_left",
4692 "left_anchor",
4693 "anchor_right",
4694 "right_anchor",
4695 "anchor-value-left",
4696];
4697
4698fn parse_ps_internal_knots(
4712 options: &BTreeMap<String, String>,
4713 degree: usize,
4714 default_internal_knots: usize,
4715) -> Result<(usize, bool, usize), String> {
4716 const MIN_EXPRESSIVE_INTERNAL_KNOTS: usize = 2;
4717 let knots_internal = if knots_option_is_list(options) {
4727 None
4728 } else {
4729 option_usize_strict(options, "knots")?
4730 };
4731 let basis_dim = option_usize_any_strict(options, &["k", "basis_dim", "basis-dim", "basisdim"])?;
4732 if knots_internal.is_some() && basis_dim.is_some() {
4733 return Err(TermBuilderError::incompatible_config(
4734 "ps/bspline smooth: specify either knots=<internal_knots> or k=<basis_dim> (not both)",
4735 )
4736 .to_string());
4737 }
4738 if let Some(k) = basis_dim {
4739 if k < 2 {
4740 return Err(TermBuilderError::invalid_option(format!(
4741 "ps/bspline smooth: k={} too small; B-spline basis requires k >= 2",
4742 k
4743 ))
4744 .to_string());
4745 }
4746 let effective_degree = degree.min(k - 1).max(1);
4752 let num_internal_knots = if effective_degree < degree {
4753 k.saturating_sub(effective_degree + 1)
4756 } else {
4757 (k - degree - 1).max(MIN_EXPRESSIVE_INTERNAL_KNOTS)
4758 };
4759 Ok((num_internal_knots, false, effective_degree))
4760 } else {
4761 Ok((
4762 knots_internal.unwrap_or(default_internal_knots),
4763 knots_internal.is_none(),
4764 degree,
4765 ))
4766 }
4767}
4768
4769fn knots_option_is_list(options: &BTreeMap<String, String>) -> bool {
4775 options
4776 .get("knots")
4777 .map(|raw| {
4778 let t = raw.trim();
4779 t.starts_with('[') || t.starts_with("c(") || t.starts_with("C(") || t.starts_with('(')
4780 })
4781 .unwrap_or(false)
4782}
4783
4784fn parse_explicit_internal_knots(
4789 options: &BTreeMap<String, String>,
4790) -> Result<Option<Vec<f64>>, String> {
4791 if !knots_option_is_list(options) {
4792 return Ok(None);
4793 }
4794 let raw = options
4795 .get("knots")
4796 .expect("knots_option_is_list implies the key is present");
4797 let tokens = split_list_option(raw);
4798 if tokens.is_empty() {
4799 return Err(TermBuilderError::invalid_option(format!(
4800 "knots={raw} is an empty list; supply at least one internal knot position \
4801 (e.g. knots=[0.2, 0.5, 0.8]) or a scalar count (e.g. knots=8)"
4802 ))
4803 .to_string());
4804 }
4805 let mut positions = Vec::with_capacity(tokens.len());
4806 for tok in &tokens {
4807 let value = parse_numeric_expr(tok).map_err(|err| {
4808 TermBuilderError::invalid_option(format!(
4809 "knots list entry '{tok}' is not a numeric position: {err}"
4810 ))
4811 .to_string()
4812 })?;
4813 positions.push(value);
4814 }
4815 Ok(Some(positions))
4816}
4817
4818fn parse_tensor_per_axis_usize(
4839 options: &BTreeMap<String, String>,
4840 key: &str,
4841 dim: usize,
4842) -> Result<Vec<Option<usize>>, String> {
4843 let Some(raw) = options.get(key) else {
4844 return Ok(vec![None; dim]);
4845 };
4846 let values = split_list_option(raw);
4847 let parse_one = |value: &str| -> Result<Option<usize>, String> {
4848 let trimmed = value.trim().trim_matches('"').trim_matches('\'').trim();
4849 if trimmed.is_empty() || trimmed.eq_ignore_ascii_case("none") {
4850 return Ok(None);
4851 }
4852 trimmed.parse::<usize>().map(Some).map_err(|err| {
4853 TermBuilderError::invalid_option(format!(
4854 "tensor smooth `{key}={raw}`: '{trimmed}' is not a non-negative integer ({err})"
4855 ))
4856 .to_string()
4857 })
4858 };
4859 if values.len() == 1 {
4860 let shared = parse_one(&values[0])?;
4861 return Ok(vec![shared; dim]);
4862 }
4863 if values.len() != dim {
4864 return Err(TermBuilderError::invalid_option(format!(
4865 "tensor smooth `{key}={raw}` has {} entries but the smooth has {dim} margins; pass one \
4866 value per margin or a single value for all of them",
4867 values.len()
4868 ))
4869 .to_string());
4870 }
4871 values.iter().map(|value| parse_one(value)).collect()
4872}
4873
4874const CR_MARGIN_DEGREE: usize = 3;
4878
4879const CR_MARGIN_PENALTY_ORDER: usize = 2;
4884
4885fn parse_knot_placement(
4886 options: &BTreeMap<String, String>,
4887) -> Result<crate::basis::BSplineKnotPlacement, String> {
4888 use crate::basis::BSplineKnotPlacement;
4889 match options
4890 .get("knot_placement")
4891 .or_else(|| options.get("knot-placement"))
4892 .or_else(|| options.get("knotplacement"))
4893 {
4894 None => Ok(BSplineKnotPlacement::Uniform),
4895 Some(raw) => match raw
4896 .trim()
4897 .trim_matches('"')
4898 .trim_matches('\'')
4899 .to_ascii_lowercase()
4900 .as_str()
4901 {
4902 "uniform" | "even" | "equal" => Ok(BSplineKnotPlacement::Uniform),
4903 "quantile" | "quantiles" | "data" | "empirical" => Ok(BSplineKnotPlacement::Quantile),
4904 other => Err(TermBuilderError::invalid_option(format!(
4905 "knot_placement={other} is not recognised; expected \"uniform\" or \"quantile\""
4906 ))
4907 .to_string()),
4908 },
4909 }
4910}
4911
4912fn explicit_knot_placement(
4922 options: &BTreeMap<String, String>,
4923) -> Result<Option<crate::basis::BSplineKnotPlacement>, String> {
4924 let declared = ["knot_placement", "knot-placement", "knotplacement"]
4925 .iter()
4926 .any(|key| options.contains_key(*key));
4927 if !declared {
4928 return Ok(None);
4929 }
4930 parse_knot_placement(options).map(Some)
4931}
4932
4933fn resolve_nonperiodic_bspline_knotspec(
4944 options: &BTreeMap<String, String>,
4945 data: ArrayView1<'_, f64>,
4946 data_range: (f64, f64),
4947 degree: usize,
4948 n_knots: usize,
4949) -> Result<BSplineKnotSpec, String> {
4950 use crate::basis::{BSplineKnotPlacement, clamped_knot_vector_from_internal_positions};
4951 if let Some(positions) = parse_explicit_internal_knots(options)? {
4952 if option_usize_any_strict(options, &["k", "basis_dim", "basis-dim", "basisdim"])?.is_some()
4953 {
4954 return Err(TermBuilderError::incompatible_config(
4955 "ps/bspline smooth: specify either explicit knots=[...] positions or \
4956 k=<basis_dim> (not both); the basis size is fixed by the knot vector",
4957 )
4958 .to_string());
4959 }
4960 let knots = clamped_knot_vector_from_internal_positions(data_range, &positions, degree)
4961 .map_err(|e| e.to_string())?;
4962 return Ok(BSplineKnotSpec::Provided(knots));
4963 }
4964 match parse_knot_placement(options)? {
4965 BSplineKnotPlacement::Uniform => Ok(BSplineKnotSpec::Generate {
4966 data_range,
4967 num_internal_knots: n_knots,
4968 }),
4969 BSplineKnotPlacement::Quantile => {
4970 crate::basis::auto_knot_vector_1d_quantile(data, n_knots, degree)
4974 .map_err(|e| e.to_string())?;
4975 Ok(BSplineKnotSpec::Automatic {
4976 num_internal_knots: Some(n_knots),
4977 placement: BSplineKnotPlacement::Quantile,
4978 })
4979 }
4980 }
4981}
4982
4983pub(crate) const SHAPE_CONSTRAINED_SMOOTH_OPTION_KEYS: &[&str] = &[
4999 "type",
5000 "bs",
5001 "k",
5002 "basis_dim",
5003 "basis-dim",
5004 "basisdim",
5005 "knots",
5006 "knot_placement",
5007 "knot-placement",
5008 "knotplacement",
5009 "degree",
5010 "penalty_order",
5011 "m",
5012 "double_penalty",
5013 "ordered",
5014];
5015
5016pub(crate) const CYCLIC_SMOOTH_OPTION_KEYS: &[&str] = &[
5017 "type",
5018 "bs",
5019 "by",
5020 "k",
5021 "basis_dim",
5022 "basis-dim",
5023 "basisdim",
5024 "degree",
5025 "penalty_order",
5026 "period",
5027 "periods",
5028 "period_start",
5029 "period_end",
5030 "start",
5031 "end",
5032 "origin",
5033 "origins",
5034 "period_origin",
5035 "period-origin",
5036 "domain_origin",
5037 "double_penalty",
5038 "id",
5039 "__by_col",
5040 "identifiability",
5041];
5042
5043pub(crate) const BSPLINE_SMOOTH_OPTION_KEYS: &[&str] = &[
5044 "type",
5045 "bs",
5046 "by",
5047 "k",
5048 "basis_dim",
5049 "basis-dim",
5050 "basisdim",
5051 "knots",
5052 "knot_placement",
5053 "knot-placement",
5054 "knotplacement",
5055 "degree",
5056 "penalty_order",
5057 "boundary",
5058 "bc",
5059 "boundary_conditions",
5060 "bc_left",
5061 "bc_right",
5062 "left_bc",
5063 "right_bc",
5064 "start_bc",
5065 "end_bc",
5066 "side",
5067 "anchor",
5068 "anchor_value",
5069 "value",
5070 "anchor_left",
5071 "left_anchor",
5072 "anchor_right",
5073 "right_anchor",
5074 "periodic",
5075 "period",
5076 "periods",
5077 "period_start",
5078 "period_end",
5079 "origin",
5080 "double_penalty",
5081 "id",
5082 "__by_col",
5083 "identifiability",
5084];
5085
5086pub(crate) const THINPLATE_SMOOTH_OPTION_KEYS: &[&str] = &[
5087 "type",
5088 "bs",
5089 "by",
5090 "length_scale",
5091 "centers",
5092 "k",
5093 "basis_dim",
5094 "basis-dim",
5095 "basisdim",
5096 "knots",
5097 "include_intercept",
5098 "double_penalty",
5099 "id",
5100 "__by_col",
5101 "identifiability",
5102 "periodic",
5103 "cyclic",
5104 "period",
5105 "period_start",
5106 "period_end",
5107 "scale_dims",
5108];
5109
5110pub(crate) const SPHERE_SMOOTH_OPTION_KEYS: &[&str] = &[
5111 "type",
5112 "bs",
5113 "by",
5114 "centers",
5115 "k",
5116 "basis_dim",
5117 "basis-dim",
5118 "basisdim",
5119 "knots",
5120 "penalty_order",
5121 "m",
5122 "double_penalty",
5123 "id",
5124 "__by_col",
5125 "kernel",
5126 "method",
5127 "radians",
5128 "units",
5129 "degree",
5130 "l",
5131 "max_degree",
5132 "max-degree",
5133 "lmax",
5134 "l_max",
5135 "l-max",
5136];
5137
5138pub(crate) const CURVATURE_SMOOTH_OPTION_KEYS: &[&str] = &[
5139 "type",
5140 "bs",
5141 "by",
5142 "centers",
5143 "k",
5144 "basis_dim",
5145 "basis-dim",
5146 "basisdim",
5147 "knots",
5148 "kappa",
5149 "length_scale",
5150 "double_penalty",
5151 "id",
5152 "__by_col",
5153];
5154
5155pub(crate) const MEASURE_JET_SMOOTH_OPTION_KEYS: &[&str] = &[
5156 "type",
5157 "bs",
5158 "by",
5159 "centers",
5160 "k",
5161 "basis_dim",
5162 "basis-dim",
5163 "basisdim",
5164 "knots",
5165 "s",
5166 "alpha",
5167 "tau",
5168 "scales",
5169 "length_scale",
5170 "double_penalty",
5171 "multiscale",
5172 "learn_length_scale",
5173 "id",
5174 "__by_col",
5175];
5176
5177pub(crate) const MATERN_SMOOTH_OPTION_KEYS: &[&str] = &[
5178 "type",
5179 "bs",
5180 "by",
5181 "nu",
5182 "length_scale",
5183 "centers",
5184 "k",
5185 "basis_dim",
5186 "basis-dim",
5187 "basisdim",
5188 "knots",
5189 "include_intercept",
5190 "double_penalty",
5191 "id",
5192 "__by_col",
5193 "identifiability",
5194 "periodic",
5195 "cyclic",
5196 "period",
5197 "period_start",
5198 "period_end",
5199 "scale_dims",
5200];
5201
5202pub(crate) const DUCHON_SMOOTH_OPTION_KEYS: &[&str] = &[
5203 "type",
5204 "bs",
5205 "by",
5206 "length_scale",
5207 "centers",
5208 "k",
5209 "basis_dim",
5210 "basis-dim",
5211 "basisdim",
5212 "knots",
5213 "rank",
5214 "power",
5215 "p",
5216 "nullspace_order",
5217 "order",
5218 "identifiability",
5219 "periodic",
5220 "cyclic",
5221 "period",
5222 "period_start",
5223 "period_end",
5224 "scale_dims",
5225 "double_penalty",
5226 "id",
5227 "__by_col",
5228];
5229
5230pub(crate) const TENSOR_SMOOTH_OPTION_KEYS: &[&str] = &[
5231 "type",
5232 "bs",
5233 "by",
5234 "k",
5235 "basis_dim",
5236 "basis-dim",
5237 "basisdim",
5238 "knot_placement",
5239 "knot-placement",
5240 "knotplacement",
5241 "degree",
5242 "penalty_order",
5243 "double_penalty",
5244 "periodic",
5245 "cyclic",
5246 "period",
5247 "periods",
5248 "period_start",
5249 "period_end",
5250 "origin",
5251 "origins",
5252 "period_origin",
5253 "period-origin",
5254 "domain_origin",
5255 "boundary",
5256 "bc",
5257 "identifiability",
5258 "id",
5259 "__by_col",
5260];
5261
5262pub(crate) const PCA_SMOOTH_OPTION_KEYS: &[&str] = &[
5263 "type",
5264 "bs",
5265 "by",
5266 "k",
5267 "basis_dim",
5268 "basis-dim",
5269 "basisdim",
5270 "lazy_path",
5271 "path",
5272 "pca_basis_path",
5273 "chunk_size",
5274 "smooth_penalty",
5275 "centered",
5276 "double_penalty",
5277 "id",
5278 "__by_col",
5279];
5280
5281pub fn validate_known_options(
5282 term_name: &str,
5283 options: &BTreeMap<String, String>,
5284 known: &[&str],
5285) -> Result<(), String> {
5286 let known_set: std::collections::BTreeSet<&&str> = known.iter().collect();
5287 for key in options.keys() {
5288 if !known_set.contains(&key.as_str()) {
5289 if term_name == "tensor" && is_tensor_k_axis_option_key(key) {
5290 continue;
5291 }
5292 let key_l = key.to_ascii_lowercase();
5294 let mut suggestions: Vec<&str> = known
5295 .iter()
5296 .filter(|k| {
5297 let kl = k.to_ascii_lowercase();
5298 kl.contains(&key_l) || key_l.contains(&kl) || {
5299 let n = kl
5300 .chars()
5301 .zip(key_l.chars())
5302 .take_while(|(a, b)| a == b)
5303 .count();
5304 n >= 3
5305 }
5306 })
5307 .copied()
5308 .collect();
5309 suggestions.sort_unstable();
5310 suggestions.dedup();
5311 let hint = if suggestions.is_empty() {
5312 String::new()
5313 } else {
5314 format!(" — did you mean one of [{}]?", suggestions.join(", "))
5315 };
5316 return Err(TermBuilderError::invalid_option(format!(
5317 "{term_name}() does not accept option `{key}`{hint}. Valid options: [{}]",
5318 {
5319 let mut sorted = known.to_vec();
5320 sorted.sort_unstable();
5321 sorted.join(", ")
5322 }
5323 ))
5324 .to_string());
5325 }
5326 }
5327 Ok(())
5328}
5329
5330pub const SECONDARY_CENTER_CAP_OPTION: &str = "__secondary_center_cap";
5340
5341pub(crate) fn cap_default_spatial_centers(
5346 options: &BTreeMap<String, String>,
5347 default_count: usize,
5348) -> usize {
5349 match option_usize(options, SECONDARY_CENTER_CAP_OPTION) {
5350 Some(cap) => default_count.min(cap),
5351 None => default_count,
5352 }
5353}
5354
5355fn default_matern_center_count(
5356 n: usize,
5357 d: usize,
5358 planned_count: usize,
5359 univariate_floor: usize,
5360) -> usize {
5361 let low_n_floor = (d + 4).min(n);
5368 planned_count
5379 .max(low_n_floor)
5380 .max(univariate_floor.min(n))
5381 .max(1)
5382}
5383
5384fn default_duchon_center_count(
5385 n: usize,
5386 d: usize,
5387 planned_count: usize,
5388 polynomial_cols: usize,
5389 univariate_floor: usize,
5390) -> usize {
5391 let mgcv_default = 10usize.saturating_mul(3usize.saturating_pow(d.saturating_sub(1) as u32));
5402 let low_n_floor = (polynomial_cols + 1).min(n).max(1);
5403 planned_count
5411 .min(mgcv_default)
5412 .max(low_n_floor)
5413 .max(univariate_floor.min(n))
5414}
5415
5416pub fn parse_countwith_basis_alias(
5417 options: &BTreeMap<String, String>,
5418 primarykey: &str,
5419 default_count: usize,
5420) -> Result<usize, String> {
5421 let primary = option_usize_strict(options, primarykey)?;
5426 let basis_dim = option_usize_any_strict(
5427 options,
5428 &["k", "basis_dim", "basis-dim", "basisdim", "knots"],
5429 )?;
5430 if primary.is_some() && basis_dim.is_some() {
5431 return Err(TermBuilderError::incompatible_config(format!(
5432 "specify either {}=<count> or k=<basis_dim> (not both)",
5433 primarykey
5434 ))
5435 .to_string());
5436 }
5437 Ok(primary.or(basis_dim).unwrap_or(default_count))
5438}
5439
5440pub fn has_explicit_countwith_basis_alias(
5441 options: &BTreeMap<String, String>,
5442 primarykey: &str,
5443) -> bool {
5444 options.contains_key(primarykey)
5445 || ["k", "basis_dim", "basis-dim", "basisdim", "knots"]
5446 .iter()
5447 .any(|alias| options.contains_key(*alias))
5448}
5449
5450pub fn parse_cyclic_boundary(
5451 options: &BTreeMap<String, String>,
5452 minv: f64,
5453 maxv: f64,
5454) -> Result<OneDimensionalBoundary, String> {
5455 let cyclic = option_bool(options, "cyclic")
5456 .or_else(|| option_bool(options, "periodic"))
5457 .unwrap_or(false);
5458 if !cyclic {
5459 return Ok(OneDimensionalBoundary::Open);
5460 }
5461 let start = match option_numeric_expr(options, "period_start")? {
5462 Some(v) => v,
5463 None => option_numeric_expr(options, "start")?.unwrap_or(minv),
5464 };
5465 let end = match option_numeric_expr(options, "period_end")? {
5466 Some(v) => v,
5467 None => option_numeric_expr(options, "end")?.unwrap_or(maxv),
5468 };
5469 if end <= start {
5470 return Err(format!(
5471 "cyclic smooth requires period_end/end ({end}) > period_start/start ({start})"
5472 ));
5473 }
5474 Ok(OneDimensionalBoundary::Cyclic { start, end })
5475}
5476
5477pub fn parse_periodic_domain_1d(
5484 options: &BTreeMap<String, String>,
5485 minv: f64,
5486 maxv: f64,
5487) -> Result<(f64, f64), String> {
5488 let start_opt = match option_numeric_expr(options, "period_start")? {
5489 Some(v) => Some(v),
5490 None => option_numeric_expr(options, "start")?,
5491 };
5492 let end_opt = match option_numeric_expr(options, "period_end")? {
5493 Some(v) => Some(v),
5494 None => option_numeric_expr(options, "end")?,
5495 };
5496 if end_opt.is_none() && start_opt.is_none() {
5507 return Err(
5508 "periodic B-spline smooth requires an explicit period: pass period=<value> \
5509 (e.g. period=2*pi) or period_start=/period_end=. Deriving the period from the \
5510 observed data range is sample-dependent and produces an off-by-ε seam, so it is \
5511 not inferred."
5512 .to_string(),
5513 );
5514 }
5515 let start = start_opt.unwrap_or(minv);
5516 let end = end_opt.unwrap_or(maxv);
5517 if !(start.is_finite() && end.is_finite()) {
5518 return Err(format!(
5519 "periodic smooth domain requires finite endpoints, got ({start}, {end})"
5520 ));
5521 }
5522 if end <= start {
5523 return Err(format!(
5524 "periodic smooth requires period_end/end ({end}) > period_start/start ({start})"
5525 ));
5526 }
5527 Ok((start, end - start))
5528}
5529
5530fn parse_matern_nu(raw: &str) -> Result<MaternNu, String> {
5531 let trimmed = raw.trim();
5532 let lowered = trimmed.to_ascii_lowercase();
5533 let named = match lowered.as_str() {
5536 "1/2" | "0.5" | "half" => Some(MaternNu::Half),
5537 "3/2" | "1.5" => Some(MaternNu::ThreeHalves),
5538 "5/2" | "2.5" => Some(MaternNu::FiveHalves),
5539 "7/2" | "3.5" => Some(MaternNu::SevenHalves),
5540 "9/2" | "4.5" => Some(MaternNu::NineHalves),
5541 _ => None,
5542 };
5543 if let Some(nu) = named {
5544 return Ok(nu);
5545 }
5546
5547 let value = if let Some((num, den)) = trimmed.split_once('/') {
5548 let num = num
5549 .trim()
5550 .parse::<f64>()
5551 .map_err(|err| format!("{}: {err}", unsupported_matern_nu_message(raw)))?;
5552 let den = den
5553 .trim()
5554 .parse::<f64>()
5555 .map_err(|err| format!("{}: {err}", unsupported_matern_nu_message(raw)))?;
5556 if den == 0.0 || !num.is_finite() || !den.is_finite() {
5557 return Err(unsupported_matern_nu_message(raw));
5558 }
5559 num / den
5560 } else {
5561 trimmed
5562 .parse::<f64>()
5563 .map_err(|err| format!("{}: {err}", unsupported_matern_nu_message(raw)))?
5564 };
5565
5566 const TOL: f64 = 1e-12;
5567 if (value - 0.5).abs() <= TOL {
5568 Ok(MaternNu::Half)
5569 } else if (value - 1.5).abs() <= TOL {
5570 Ok(MaternNu::ThreeHalves)
5571 } else if (value - 2.5).abs() <= TOL {
5572 Ok(MaternNu::FiveHalves)
5573 } else if (value - 3.5).abs() <= TOL {
5574 Ok(MaternNu::SevenHalves)
5575 } else if (value - 4.5).abs() <= TOL {
5576 Ok(MaternNu::NineHalves)
5577 } else {
5578 Err(unsupported_matern_nu_message(raw))
5579 }
5580}
5581
5582fn unsupported_matern_nu_message(raw: &str) -> String {
5583 TermBuilderError::unsupported_feature(format!(
5584 "unsupported Matern nu '{raw}'; supported half-integer values are 1/2, 3/2, 5/2, 7/2, and 9/2"
5585 ))
5586 .to_string()
5587}
5588
5589#[derive(Clone, Debug, serde::Serialize, serde::Deserialize)]
5590pub enum DuchonPowerPolicy {
5591 Explicit(f64),
5592 CubicStructuralDefault,
5596}
5597
5598pub fn parse_duchon_power_policy(
5599 options: &BTreeMap<String, String>,
5600) -> Result<DuchonPowerPolicy, String> {
5601 if let Some(raw_nu) = options.get("nu") {
5602 return Err(TermBuilderError::incompatible_config(format!(
5603 "Duchon smooths use power=<number>, not nu='{}'. Use power=1.5, power=2, etc.",
5604 raw_nu
5605 ))
5606 .to_string());
5607 }
5608 match options.get("power").or_else(|| options.get("p")) {
5612 Some(raw) => {
5613 let value = raw.parse::<f64>().map_err(|err| {
5614 TermBuilderError::invalid_option(format!(
5615 "invalid Duchon power '{}'; expected a non-negative number such as power=1.5 or power=2: {}",
5616 raw, err
5617 ))
5618 .to_string()
5619 })?;
5620 if !value.is_finite() || value < 0.0 {
5621 return Err(TermBuilderError::invalid_option(format!(
5622 "invalid Duchon power '{}'; expected a finite non-negative number such as power=1.5 or power=2",
5623 raw
5624 ))
5625 .to_string());
5626 }
5627 Ok(DuchonPowerPolicy::Explicit(value))
5628 }
5629 None => Ok(DuchonPowerPolicy::CubicStructuralDefault),
5630 }
5631}
5632
5633pub fn parse_duchon_power(options: &BTreeMap<String, String>) -> Result<f64, String> {
5634 match parse_duchon_power_policy(options)? {
5635 DuchonPowerPolicy::Explicit(power) => Ok(power),
5636 DuchonPowerPolicy::CubicStructuralDefault => Ok(1.5),
5642 }
5643}
5644
5645pub fn parse_duchon_order_opt(
5654 options: &BTreeMap<String, String>,
5655) -> Result<Option<DuchonNullspaceOrder>, String> {
5656 if !options.contains_key("order") && !options.contains_key("nullspace_order") {
5657 return Ok(None);
5658 }
5659 parse_duchon_order(options).map(Some)
5660}
5661
5662pub fn parse_duchon_order(
5663 options: &BTreeMap<String, String>,
5664) -> Result<DuchonNullspaceOrder, String> {
5665 match options.get("order").or_else(|| options.get("nullspace_order")) {
5668 None => Ok(DuchonNullspaceOrder::Linear),
5672 Some(raw) => match raw.parse::<usize>() {
5673 Ok(0) => Ok(DuchonNullspaceOrder::Zero),
5674 Ok(1) => Ok(DuchonNullspaceOrder::Linear),
5675 Ok(other) => Ok(DuchonNullspaceOrder::Degree(other)),
5676 Err(_) => Err(TermBuilderError::invalid_option(format!(
5677 "invalid Duchon order '{}'; expected a non-negative integer such as order=0, order=1, or order=2",
5678 raw
5679 ))
5680 .to_string()),
5681 },
5682 }
5683}
5684
5685fn parse_matern_identifiability(
5686 options: &BTreeMap<String, String>,
5687) -> Result<MaternIdentifiability, TermBuilderError> {
5688 let Some(raw) = options.get("identifiability").map(String::as_str) else {
5689 return Ok(MaternIdentifiability::default());
5690 };
5691 match raw.trim().to_ascii_lowercase().as_str() {
5692 "none" => Ok(MaternIdentifiability::None),
5693 "sum_tozero" | "sum-to-zero" | "center_sum_tozero" | "center-sum-to-zero" | "centered" => {
5694 Ok(MaternIdentifiability::CenterSumToZero)
5695 }
5696 "linear" | "center_linear_orthogonal" | "center-linear-orthogonal" => {
5697 Ok(MaternIdentifiability::CenterLinearOrthogonal)
5698 }
5699 other => Err(TermBuilderError::unsupported_feature(format!(
5700 "invalid Matérn identifiability '{other}'; expected one of: none, sum_tozero, linear"
5701 ))),
5702 }
5703}
5704
5705fn parse_spatial_identifiability(
5706 options: &BTreeMap<String, String>,
5707) -> Result<SpatialIdentifiability, TermBuilderError> {
5708 let Some(raw) = options.get("identifiability").map(String::as_str) else {
5709 return Ok(SpatialIdentifiability::default());
5710 };
5711 match raw.trim().to_ascii_lowercase().as_str() {
5712 "none" => Ok(SpatialIdentifiability::None),
5713 "orthogonal"
5714 | "orthogonal_to_parametric"
5715 | "orthogonal-to-parametric"
5716 | "parametric_orthogonal" => Ok(SpatialIdentifiability::OrthogonalToParametric),
5717 "frozen" => Err(TermBuilderError::unsupported_feature(
5718 "spatial identifiability 'frozen' is internal-only; use none or orthogonal_to_parametric",
5719 )),
5720 other => Err(TermBuilderError::unsupported_feature(format!(
5721 "invalid spatial identifiability '{other}'; expected one of: none, orthogonal_to_parametric"
5722 ))),
5723 }
5724}
5725
5726#[cfg(test)]
5727mod tests {
5728 use super::*;
5729 use crate::basis::{OperatorPenaltySpec, PenaltySource};
5730 use crate::inference::formula_dsl::parse_formula;
5731 use gam_data::{DataSchema, SchemaColumn};
5732 use ndarray::{Array1, Array2};
5733 use std::collections::BTreeMap;
5734
5735 #[test]
5742 fn unique_count_column_uses_canonical_level_bits() {
5743 let signed_zero = Array1::from(vec![0.0, -0.0, 0.0]);
5745 assert_eq!(
5746 unique_count_column(signed_zero.view()),
5747 1,
5748 "+0.0 and -0.0 must collapse to a single level"
5749 );
5750
5751 let nan_a = f64::from_bits(0x7ff8_0000_0000_0001);
5752 let nan_b = f64::from_bits(0xfff8_0000_0000_dead);
5753 assert!(nan_a.is_nan() && nan_b.is_nan() && nan_a.to_bits() != nan_b.to_bits());
5754 let nans = Array1::from(vec![nan_a, nan_b]);
5755 assert_eq!(
5756 unique_count_column(nans.view()),
5757 1,
5758 "distinct NaN payloads must collapse to a single level"
5759 );
5760
5761 let finite = Array1::from(vec![1.0, 2.0, 2.0, 3.0]);
5763 assert_eq!(unique_count_column(finite.view()), 3);
5764 }
5765
5766 #[test]
5775 fn radial_1d_default_not_starved_below_univariate_spline_resolution_1867() {
5776 let n = 30usize;
5777 let d = 1usize;
5778 let planned = default_num_centers(n, d);
5780 assert!(
5781 planned < 11,
5782 "precondition: conditioning cap starves the raw radial default (got {planned})"
5783 );
5784 let col: Array1<f64> = Array1::from_iter((0..n).map(|i| i as f64 / (n as f64 - 1.0)));
5787 let univariate_floor =
5788 heuristic_knots_for_column(col.view()).saturating_add(DEFAULT_BSPLINE_DEGREE + 1);
5789 assert_eq!(univariate_floor, 11, "univariate spline resolution at n=30");
5790
5791 assert_eq!(default_matern_center_count(n, d, planned, 0), planned);
5793 assert!(default_duchon_center_count(n, d, planned, 2, 0) <= planned);
5794
5795 assert!(
5799 default_matern_center_count(n, d, planned, univariate_floor) >= univariate_floor,
5800 "matern 1-D default must not be starved below the spline resolution"
5801 );
5802 assert!(
5803 default_duchon_center_count(n, d, planned, 2, univariate_floor) >= univariate_floor,
5804 "duchon 1-D default must not be starved below the spline resolution"
5805 );
5806
5807 assert_eq!(default_matern_center_count(200, 2, 40, 0), 40);
5810 }
5811
5812 #[test]
5821 fn duchon_2d_default_is_low_rank_not_generic_spatial_width_1757() {
5822 let n = 500usize;
5823 let d = 2usize;
5824 let polynomial_cols = d + 1;
5825 let generic_plan = default_num_centers(n, d);
5826 let duchon_default = default_duchon_center_count(n, d, generic_plan, polynomial_cols, 0);
5827 let spline_rank = 10usize.saturating_mul(3usize.saturating_pow((d - 1) as u32));
5828
5829 assert!(
5830 generic_plan > spline_rank,
5831 "precondition: generic spatial plan should be wider than the Duchon low-rank spline rank"
5832 );
5833 assert_eq!(
5834 duchon_default, spline_rank,
5835 "2-D Duchon default must use the low-rank spline representer size, not the generic spatial width"
5836 );
5837 assert!(
5838 duchon_default > polynomial_cols,
5839 "the capped default must still contain the affine polynomial null space"
5840 );
5841 }
5842
5843 #[test]
5856 fn measure_jet_reml_selects_the_representer_range_by_default_2761() {
5857 let ds = continuous_dataset(
5858 &["y", "x1", "x2"],
5859 (0..40)
5860 .map(|i| {
5861 let t = i as f64 / 39.0;
5862 vec![(6.0 * t).sin(), t, 0.5 + 0.5 * (6.0 * t).cos()]
5863 })
5864 .collect(),
5865 );
5866 let col_map = ds.column_map();
5867 let learns = |body: &str| -> bool {
5868 let parsed = parse_formula(&format!("y ~ {body}")).expect("parse mjs formula");
5869 let terms = build_termspec(
5870 &parsed.terms,
5871 &ds,
5872 &col_map,
5873 &mut Vec::new(),
5874 &gam_runtime::resource::ResourcePolicy::default_library(),
5875 )
5876 .expect("build mjs term");
5877 let SmoothBasisSpec::MeasureJet { spec, .. } = &terms.smooth_terms[0].basis else {
5878 panic!("expected a measure-jet smooth for '{body}'");
5879 };
5880 let learns = crate::smooth::measure_jet_learns_length_scale(spec);
5883 assert_eq!(
5884 spec.learn_length_scale, learns,
5885 "'{body}': the ψ accessor and the spec field must not disagree"
5886 );
5887 assert_eq!(
5888 crate::smooth::measure_jet_psi_dim(spec),
5889 usize::from(learns),
5890 "'{body}': single-scale ψ dimension is exactly the ℓ coordinate"
5891 );
5892 assert_eq!(
5893 crate::smooth::measure_jet_enrolls_psi(spec),
5894 learns,
5895 "'{body}': single-scale enrollment is exactly the ℓ coordinate"
5896 );
5897 learns
5898 };
5899
5900 assert!(
5901 learns("mjs(x1, x2, centers=8)"),
5902 "a plain measure-jet smooth must REML-select its representer range: λ shrinks \
5903 inside a span and cannot move one, so a frozen ℓ is an error no smoothing \
5904 parameter can repair (#2761 measured 13.4x held-out RMSE, with the design's \
5905 own least-squares span floor sitting AT the fitted value)"
5906 );
5907 assert!(
5908 !learns("mjs(x1, x2, centers=8, length_scale=0.3)"),
5909 "a typed length_scale= is a request, not a seed, and must pin ℓ — the same \
5910 short-circuit an explicitly-scaled Matérn gets"
5911 );
5912 assert!(
5913 !learns("mjs(x1, x2, centers=8, learn_length_scale=false)"),
5914 "an explicit opt-out must be honored"
5915 );
5916 assert!(
5917 learns("mjs(x1, x2, centers=8, length_scale=0.3, learn_length_scale=true)"),
5918 "an explicit opt-in must beat the length_scale= pin, so a caller can seed the \
5919 search at a range of their choosing"
5920 );
5921 }
5922
5923 fn continuous_dataset(headers: &[&str], rows: Vec<Vec<f64>>) -> Dataset {
5924 let nrows = rows.len();
5925 let ncols = headers.len();
5926 let values = Array2::from_shape_vec(
5927 (nrows, ncols),
5928 rows.into_iter().flat_map(|row| row.into_iter()).collect(),
5929 )
5930 .expect("rectangular test data");
5931 Dataset {
5932 headers: headers.iter().map(|name| name.to_string()).collect(),
5933 values,
5934 schema: DataSchema {
5935 columns: headers
5936 .iter()
5937 .map(|name| SchemaColumn {
5938 name: name.to_string(),
5939 kind: ColumnKindTag::Continuous,
5940 levels: vec![],
5941 })
5942 .collect(),
5943 },
5944 column_kinds: vec![ColumnKindTag::Continuous; ncols],
5945 }
5946 }
5947
5948 fn factor_dataset() -> Dataset {
5949 let rows = (0..24)
5950 .map(|i| {
5951 let x = i as f64 / 23.0;
5952 let g = (i % 2) as f64;
5953 vec![x + g, x, g]
5954 })
5955 .collect::<Vec<_>>();
5956 Dataset {
5957 headers: vec!["y".into(), "x".into(), "g".into()],
5958 values: Array2::from_shape_vec(
5959 (rows.len(), 3),
5960 rows.into_iter().flat_map(|row| row.into_iter()).collect(),
5961 )
5962 .expect("rectangular factor test data"),
5963 schema: DataSchema {
5964 columns: vec![
5965 SchemaColumn {
5966 name: "y".into(),
5967 kind: ColumnKindTag::Continuous,
5968 levels: vec![],
5969 },
5970 SchemaColumn {
5971 name: "x".into(),
5972 kind: ColumnKindTag::Continuous,
5973 levels: vec![],
5974 },
5975 SchemaColumn {
5976 name: "g".into(),
5977 kind: ColumnKindTag::Categorical,
5978 levels: vec!["a".into(), "b".into()],
5979 },
5980 ],
5981 },
5982 column_kinds: vec![
5983 ColumnKindTag::Continuous,
5984 ColumnKindTag::Continuous,
5985 ColumnKindTag::Categorical,
5986 ],
5987 }
5988 }
5989
5990 fn build_two_dimensional_spatial_basis(
5991 ds: &Dataset,
5992 selector: &str,
5993 count_option: Option<&str>,
5994 ) -> SmoothBasisSpec {
5995 let mut options = BTreeMap::new();
5996 options.insert("bs".to_string(), selector.to_string());
5997 if let Some(option) = count_option {
5998 options.insert(option.to_string(), "7".to_string());
5999 }
6000 let mut notes = Vec::new();
6001 build_smooth_basis(
6002 SmoothKind::S,
6003 &["x".to_string(), "z".to_string()],
6004 &[1, 2],
6005 &options,
6006 ds,
6007 &mut notes,
6008 &ResourcePolicy::default_library(),
6009 1,
6010 )
6011 .unwrap_or_else(|error| {
6012 panic!("failed to build {selector} with count option {count_option:?}: {error}")
6013 })
6014 }
6015
6016 fn curvature_or_measurejet_center_strategy(basis: &SmoothBasisSpec) -> &CenterStrategy {
6017 match basis {
6018 SmoothBasisSpec::ConstantCurvature { spec, .. } => &spec.center_strategy,
6019 SmoothBasisSpec::MeasureJet { spec, .. } => &spec.center_strategy,
6020 other => panic!("expected curvature or measure-jet basis, got {other:?}"),
6021 }
6022 }
6023
6024 fn build_sphere_over_lat_lon(ds: &Dataset) -> Result<SmoothBasisSpec, String> {
6026 let mut options = BTreeMap::new();
6027 options.insert("bs".to_string(), "sphere".to_string());
6028 options.insert("k".to_string(), "10".to_string());
6029 options.insert("kernel".to_string(), "sobolev".to_string());
6030 let mut notes = Vec::new();
6031 build_smooth_basis(
6032 SmoothKind::S,
6033 &["lat".to_string(), "lon".to_string()],
6034 &[1, 2],
6035 &options,
6036 ds,
6037 &mut notes,
6038 &ResourcePolicy::default_library(),
6039 1,
6040 )
6041 }
6042
6043 #[test]
6049 fn sphere_rejects_constant_longitude_but_accepts_varying() {
6050 let rows_const_lon: Vec<Vec<f64>> = (0..60)
6052 .map(|i| {
6053 let lat = -70.0 + 140.0 * (i as f64) / 59.0;
6054 vec![0.0, lat, 0.0] })
6056 .collect();
6057 let ds_const = continuous_dataset(&["y", "lat", "lon"], rows_const_lon);
6058 let err = build_sphere_over_lat_lon(&ds_const)
6059 .expect_err("a constant-longitude sphere smooth must be rejected as degenerate");
6060 let lower = err.to_lowercase();
6061 assert!(
6062 (lower.contains("constant")
6063 || lower.contains("degenerate")
6064 || lower.contains("unique"))
6065 && lower.contains("lon"),
6066 "rejection must flag degeneracy and name the constant longitude coordinate: {err}"
6067 );
6068
6069 let rows_ok: Vec<Vec<f64>> = (0..60)
6071 .map(|i| {
6072 let lat = -70.0 + 140.0 * (i as f64) / 59.0;
6073 let lon = -170.0 + 340.0 * ((i * 17 % 60) as f64) / 59.0;
6076 vec![0.0, lat, lon]
6077 })
6078 .collect();
6079 let ds_ok = continuous_dataset(&["y", "lat", "lon"], rows_ok);
6080 build_sphere_over_lat_lon(&ds_ok)
6081 .expect("a sphere smooth over varying latitude and longitude must build");
6082 }
6083
6084 #[test]
6085 fn curvature_and_measurejet_omitted_counts_retain_auto_provenance() {
6086 let ds = continuous_dataset(
6087 &["y", "x", "z"],
6088 (0..64)
6089 .map(|i| {
6090 let x = i as f64 / 63.0;
6091 let z = ((i * 17) % 64) as f64 / 63.0;
6092 vec![x.sin() + z.cos(), x, z]
6093 })
6094 .collect(),
6095 );
6096 let expected = default_num_centers(ds.values.nrows(), 2);
6097
6098 for selector in ["curv", "mjs"] {
6099 let basis = build_two_dimensional_spatial_basis(&ds, selector, None);
6100 let strategy = curvature_or_measurejet_center_strategy(&basis);
6101 assert!(
6102 matches!(strategy, CenterStrategy::Auto(_)),
6103 "an omitted count on {selector} must retain Auto provenance, got {strategy:?}",
6104 );
6105 assert_eq!(
6106 strategy.planned_num_centers(2),
6107 expected,
6108 "Auto provenance must preserve {selector}'s resolved default count",
6109 );
6110 }
6111 }
6112
6113 #[test]
6114 fn curvature_and_measurejet_explicit_count_aliases_remain_pinned() {
6115 let ds = continuous_dataset(
6116 &["y", "x", "z"],
6117 (0..32)
6118 .map(|i| {
6119 let x = i as f64 / 31.0;
6120 let z = ((i * 11) % 32) as f64 / 31.0;
6121 vec![x - z, x, z]
6122 })
6123 .collect(),
6124 );
6125
6126 for selector in ["curv", "mjs"] {
6127 for alias in [
6128 "centers",
6129 "k",
6130 "basis_dim",
6131 "basis-dim",
6132 "basisdim",
6133 "knots",
6134 ] {
6135 let basis = build_two_dimensional_spatial_basis(&ds, selector, Some(alias));
6136 let strategy = curvature_or_measurejet_center_strategy(&basis);
6137 assert!(
6138 !matches!(strategy, CenterStrategy::Auto(_)),
6139 "explicit {alias}= on {selector} must remain pinned, got {strategy:?}",
6140 );
6141 assert_eq!(
6142 strategy.planned_num_centers(2),
6143 7,
6144 "explicit {alias}= must remain the exact {selector} center count",
6145 );
6146 }
6147 }
6148 }
6149
6150 #[test]
6158 fn default_univariate_thinplate_basis_dim_is_modest() {
6159 let n = 300usize;
6162 let rows: Vec<Vec<f64>> = (0..n)
6163 .map(|i| {
6164 let x = -3.0 + 6.0 * (i as f64) / ((n - 1) as f64);
6165 vec![x.sin(), x]
6166 })
6167 .collect();
6168 let ds = continuous_dataset(&["y", "x"], rows);
6169
6170 let mut options = BTreeMap::new();
6171 options.insert("bs".to_string(), "tp".to_string());
6172
6173 let mut notes = Vec::new();
6174 let basis = build_smooth_basis(
6175 SmoothKind::S,
6176 &["x".to_string()],
6177 &[1],
6178 &options,
6179 &ds,
6180 &mut notes,
6181 &ResourcePolicy::default_library(),
6182 1,
6183 )
6184 .expect("build default univariate tp smooth");
6185
6186 let centers = match &basis {
6187 SmoothBasisSpec::ThinPlate { spec, .. } => match &spec.center_strategy {
6188 CenterStrategy::Auto(inner) => match inner.as_ref() {
6189 CenterStrategy::FarthestPoint { num_centers }
6190 | CenterStrategy::EqualMass { num_centers }
6191 | CenterStrategy::EqualMassCovarRepresentative { num_centers }
6192 | CenterStrategy::KMeans { num_centers, .. } => *num_centers,
6193 other => panic!("unexpected auto inner center strategy: {other:?}"),
6194 },
6195 CenterStrategy::FarthestPoint { num_centers }
6196 | CenterStrategy::EqualMass { num_centers }
6197 | CenterStrategy::EqualMassCovarRepresentative { num_centers }
6198 | CenterStrategy::KMeans { num_centers, .. } => *num_centers,
6199 other => panic!("unexpected center strategy: {other:?}"),
6200 },
6201 other => panic!("expected ThinPlate basis, got {other:?}"),
6202 };
6203
6204 assert!(
6208 centers >= 1,
6209 "default univariate tp must still build a usable basis (centers={centers})",
6210 );
6211 }
6212
6213 #[test]
6220 fn default_matern_2d_seeds_resolving_length_scale_not_overscaled_diameter() {
6221 let side = 24usize; let mut rows: Vec<Vec<f64>> = Vec::with_capacity(side * side);
6226 for i in 0..side {
6227 for j in 0..side {
6228 let x1 = i as f64 / (side - 1) as f64; let x2 = j as f64 / (side - 1) as f64; let y = (6.0 * x1).sin() * (6.0 * x2).cos();
6231 rows.push(vec![y, x1, x2]);
6232 }
6233 }
6234 let n = rows.len();
6235 let ds = continuous_dataset(&["y", "x1", "x2"], rows);
6236
6237 let mut options = BTreeMap::new();
6238 options.insert("bs".to_string(), "gp".to_string()); let mut notes = Vec::new();
6240 let mut basis = build_smooth_basis(
6241 SmoothKind::S,
6242 &["x1".to_string(), "x2".to_string()],
6243 &[1, 2],
6244 &options,
6245 &ds,
6246 &mut notes,
6247 &ResourcePolicy::default_library(),
6248 1,
6249 )
6250 .expect("build default 2-D matern smooth");
6251
6252 let (feature_cols, seeded_length_scale) = match &basis {
6255 SmoothBasisSpec::Matern {
6256 feature_cols, spec, ..
6257 } => (feature_cols.clone(), spec.length_scale),
6258 other => panic!("expected Matern basis, got {other:?}"),
6259 };
6260 assert_eq!(seeded_length_scale, MaternLengthScale::auto());
6261
6262 crate::smooth::auto_init_length_scale_in_basis(ds.values.view(), &mut basis);
6271 let (realized, requested_centers) = match &basis {
6272 SmoothBasisSpec::Matern { spec, .. } => (
6273 spec.length_scale
6274 .resolved()
6275 .expect("auto-init must resolve Matérn length scale"),
6276 match &spec.center_strategy {
6277 CenterStrategy::FarthestPoint { num_centers }
6278 | CenterStrategy::EqualMass { num_centers }
6279 | CenterStrategy::EqualMassCovarRepresentative { num_centers }
6280 | CenterStrategy::KMeans { num_centers, .. } => *num_centers,
6281 CenterStrategy::Auto(inner) => match inner.as_ref() {
6282 CenterStrategy::FarthestPoint { num_centers }
6283 | CenterStrategy::EqualMass { num_centers }
6284 | CenterStrategy::EqualMassCovarRepresentative { num_centers }
6285 | CenterStrategy::KMeans { num_centers, .. } => *num_centers,
6286 other => panic!("unexpected inner center strategy: {other:?}"),
6287 },
6288 other => panic!("unexpected center strategy: {other:?}"),
6289 },
6290 ),
6291 other => panic!("expected Matern basis after auto-init, got {other:?}"),
6292 };
6293 let expected = crate::smooth::auto_initial_length_scale_for_centers(
6294 ds.values.view(),
6295 &feature_cols,
6296 requested_centers,
6297 );
6298 assert!(
6299 (realized - expected).abs() <= 1e-12,
6300 "auto-init must seed the density-adaptive rotation-invariant \
6301 wiggly-side length scale (expected {expected}, got {realized})",
6302 );
6303
6304 let max_range = 1.0_f64; assert!(
6309 realized < max_range / 4.0,
6310 "matern seed length_scale {realized} must be in the resolving regime, \
6311 not the over-smoothed diameter corner (n={n}, max_range≈{max_range})",
6312 );
6313 }
6314
6315 #[test]
6319 fn matern_length_scale_provenance_drives_prebuild_kappa_locking() {
6320 let ds = continuous_dataset(
6321 &["y", "x1", "x2"],
6322 vec![
6323 vec![0.0, -1.0, -0.5],
6324 vec![1.0, -0.2, 0.7],
6325 vec![0.0, 0.6, -0.8],
6326 vec![1.0, 1.1, 0.4],
6327 ],
6328 );
6329 let build = |length_scale: Option<&str>| {
6330 let mut options = BTreeMap::new();
6331 options.insert("bs".to_string(), "gp".to_string());
6332 if let Some(value) = length_scale {
6333 options.insert("length_scale".to_string(), value.to_string());
6334 }
6335 let mut notes = Vec::new();
6336 build_smooth_basis(
6337 SmoothKind::S,
6338 &["x1".to_string(), "x2".to_string()],
6339 &[1, 2],
6340 &options,
6341 &ds,
6342 &mut notes,
6343 &ResourcePolicy::default_library(),
6344 1,
6345 )
6346 .expect("build Matérn provenance fixture")
6347 };
6348 let collection = |basis| TermCollectionSpec {
6349 linear_terms: Vec::new(),
6350 random_effect_terms: Vec::new(),
6351 smooth_terms: vec![SmoothTermSpec {
6352 frozen_parametric_residualization: None,
6353 name: "spatial".to_string(),
6354 basis,
6355 shape: ShapeConstraint::None,
6356 joint_null_rotation: None,
6357 }],
6358 };
6359
6360 let mut auto = collection(build(None));
6361 assert!(matches!(
6362 &auto.smooth_terms[0].basis,
6363 SmoothBasisSpec::Matern {
6364 spec: MaternBasisSpec {
6365 length_scale: MaternLengthScale::Auto { resolved: None },
6366 ..
6367 },
6368 ..
6369 }
6370 ));
6371 assert!(
6372 !crate::smooth::all_spatial_terms_kappa_fixed(&auto),
6373 "BMS pre-design query must enroll omitted Matérn κ"
6374 );
6375 crate::smooth::auto_init_length_scale_in_place(ds.values.view(), &mut auto.smooth_terms[0]);
6376 assert!(matches!(
6377 &auto.smooth_terms[0].basis,
6378 SmoothBasisSpec::Matern {
6379 spec: MaternBasisSpec {
6380 length_scale: MaternLengthScale::Auto {
6381 resolved: Some(value)
6382 },
6383 ..
6384 },
6385 ..
6386 } if value.is_finite() && *value > 0.0
6387 ));
6388 assert!(
6389 !crate::smooth::all_spatial_terms_kappa_fixed(&auto),
6390 "resolved Auto Matérn κ must remain optimizer-owned"
6391 );
6392
6393 for explicit in ["0.75", "0.0"] {
6394 let fixed = collection(build(Some(explicit)));
6395 assert!(matches!(
6396 &fixed.smooth_terms[0].basis,
6397 SmoothBasisSpec::Matern {
6398 spec: MaternBasisSpec {
6399 length_scale: MaternLengthScale::Fixed(value),
6400 ..
6401 },
6402 ..
6403 } if *value == explicit.parse::<f64>().unwrap()
6404 ));
6405 assert!(
6406 crate::smooth::all_spatial_terms_kappa_fixed(&fixed),
6407 "explicit Matérn length_scale={explicit} must lock κ before design build"
6408 );
6409 }
6410 }
6411
6412 #[test]
6420 fn matern_and_thinplate_accept_periodic_option() {
6421 let n = 200usize;
6422 let rows: Vec<Vec<f64>> = (0..n)
6423 .map(|i| {
6424 let x = -3.0 + 6.0 * (i as f64) / ((n - 1) as f64);
6425 vec![x.sin(), x]
6426 })
6427 .collect();
6428 let ds = continuous_dataset(&["y", "x"], rows);
6429
6430 let mut matern_opts = BTreeMap::new();
6432 matern_opts.insert("bs".to_string(), "gp".to_string()); matern_opts.insert("periodic".to_string(), "true".to_string());
6434 let mut notes = Vec::new();
6435 let matern_basis = build_smooth_basis(
6436 SmoothKind::S,
6437 &["x".to_string()],
6438 &[1],
6439 &matern_opts,
6440 &ds,
6441 &mut notes,
6442 &ResourcePolicy::default_library(),
6443 1,
6444 )
6445 .expect("matern(x, periodic=true) must be accepted");
6446 match &matern_basis {
6447 SmoothBasisSpec::Matern { spec, .. } => assert!(
6448 spec.periodic.is_some(),
6449 "periodic=true must thread a Some(periodic) into the matern spec",
6450 ),
6451 other => panic!("expected Matern basis, got {other:?}"),
6452 }
6453
6454 let mut tps_opts = BTreeMap::new();
6456 tps_opts.insert("bs".to_string(), "tp".to_string());
6457 tps_opts.insert("periodic".to_string(), "true".to_string());
6458 let mut notes = Vec::new();
6459 let tps_basis = build_smooth_basis(
6460 SmoothKind::S,
6461 &["x".to_string()],
6462 &[1],
6463 &tps_opts,
6464 &ds,
6465 &mut notes,
6466 &ResourcePolicy::default_library(),
6467 1,
6468 )
6469 .expect("thinplate(x, periodic=true) must be accepted");
6470 match &tps_basis {
6471 SmoothBasisSpec::ThinPlate { spec, .. } => assert!(
6472 spec.periodic.is_some(),
6473 "periodic=true must thread a Some(periodic) into the thinplate spec",
6474 ),
6475 other => panic!("expected ThinPlate basis, got {other:?}"),
6476 }
6477 }
6478
6479 #[test]
6488 fn scalar_periodic_false_builds_non_periodic_radial_smooth() {
6489 let n = 200usize;
6490 let rows: Vec<Vec<f64>> = (0..n)
6491 .map(|i| {
6492 let x = -3.0 + 6.0 * (i as f64) / ((n - 1) as f64);
6493 vec![x.sin(), x]
6494 })
6495 .collect();
6496 let ds = continuous_dataset(&["y", "x"], rows);
6497
6498 let build = |bs: &str| -> SmoothBasisSpec {
6499 let mut opts = BTreeMap::new();
6500 opts.insert("bs".to_string(), bs.to_string());
6501 opts.insert("periodic".to_string(), "false".to_string());
6502 let mut notes = Vec::new();
6503 build_smooth_basis(
6504 SmoothKind::S,
6505 &["x".to_string()],
6506 &[1],
6507 &opts,
6508 &ds,
6509 &mut notes,
6510 &ResourcePolicy::default_library(),
6511 1,
6512 )
6513 .unwrap_or_else(|e| panic!("s(x, bs={bs}, periodic=false) must be accepted: {e}"))
6514 };
6515
6516 match &build("gp") {
6517 SmoothBasisSpec::Matern { spec, .. } => assert!(
6518 spec.periodic.is_none(),
6519 "periodic=false must leave the matern spec non-periodic, got {:?}",
6520 spec.periodic
6521 ),
6522 other => panic!("expected Matern basis, got {other:?}"),
6523 }
6524 match &build("tp") {
6525 SmoothBasisSpec::ThinPlate { spec, .. } => assert!(
6526 spec.periodic.is_none(),
6527 "periodic=false must leave the thinplate spec non-periodic, got {:?}",
6528 spec.periodic
6529 ),
6530 other => panic!("expected ThinPlate basis, got {other:?}"),
6531 }
6532 match &build("duchon") {
6533 SmoothBasisSpec::Duchon { spec, .. } => assert!(
6534 spec.periodic.is_none(),
6535 "periodic=false must leave the duchon spec non-periodic (no data-range \
6536 back-fill), got {:?}",
6537 spec.periodic
6538 ),
6539 other => panic!("expected Duchon basis, got {other:?}"),
6540 }
6541 }
6542
6543 fn inferred_tensor_basis_product(ds: &Dataset) -> usize {
6544 let parsed = parse_formula("y ~ te(theta, h)").expect("parse tensor formula");
6545 let col_map = ds.column_map();
6546 let mut notes = Vec::new();
6547 let terms = build_termspec(
6548 &parsed.terms,
6549 ds,
6550 &col_map,
6551 &mut notes,
6552 &ResourcePolicy::default_library(),
6553 )
6554 .expect("build tensor termspec");
6555 let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
6556 panic!("expected tensor smooth");
6557 };
6558 spec.marginalspecs
6559 .iter()
6560 .map(|marginal| match marginal.knotspec {
6561 BSplineKnotSpec::Generate {
6562 num_internal_knots, ..
6563 } => num_internal_knots + marginal.degree + 1,
6564 BSplineKnotSpec::PeriodicUniform { num_basis, .. } => num_basis,
6565 BSplineKnotSpec::Automatic {
6566 num_internal_knots: Some(num_internal_knots),
6567 ..
6568 } => num_internal_knots + marginal.degree + 1,
6569 BSplineKnotSpec::Automatic {
6570 num_internal_knots: None,
6571 ..
6572 } => panic!("test helper cannot infer automatic knot count"),
6573 BSplineKnotSpec::Provided(ref knots) => {
6574 knots.len().saturating_sub(marginal.degree + 1)
6575 }
6576 BSplineKnotSpec::NaturalCubicRegression { ref knots } => knots.len(),
6578 })
6579 .product()
6580 }
6581
6582 fn tensor_margin_basis_sizes(ds: &Dataset, formula: &str) -> Vec<usize> {
6583 let parsed = parse_formula(formula).expect("parse tensor formula");
6584 let col_map = ds.column_map();
6585 let mut notes = Vec::new();
6586 let terms = build_termspec(
6587 &parsed.terms,
6588 ds,
6589 &col_map,
6590 &mut notes,
6591 &ResourcePolicy::default_library(),
6592 )
6593 .expect("build tensor termspec");
6594 let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
6595 panic!("expected tensor smooth");
6596 };
6597 spec.marginalspecs
6598 .iter()
6599 .map(|marginal| match marginal.knotspec {
6600 BSplineKnotSpec::Generate {
6601 num_internal_knots, ..
6602 } => num_internal_knots + marginal.degree + 1,
6603 BSplineKnotSpec::PeriodicUniform { num_basis, .. } => num_basis,
6604 BSplineKnotSpec::Automatic {
6605 num_internal_knots: Some(num_internal_knots),
6606 ..
6607 } => num_internal_knots + marginal.degree + 1,
6608 BSplineKnotSpec::Automatic {
6609 num_internal_knots: None,
6610 ..
6611 } => panic!("test helper cannot infer automatic knot count"),
6612 BSplineKnotSpec::Provided(ref knots) => {
6613 knots.len().saturating_sub(marginal.degree + 1)
6614 }
6615 BSplineKnotSpec::NaturalCubicRegression { ref knots } => knots.len(),
6617 })
6618 .collect()
6619 }
6620
6621 #[test]
6622 fn validate_known_options_lists_valid_option_names_for_unknown_parameter() {
6623 let mut options = BTreeMap::new();
6624 options.insert("lengt_scale".to_string(), "0.25".to_string());
6625 let err = validate_known_options(
6626 "matern",
6627 &options,
6628 &["type", "bs", "length_scale", "centers", "k", "nu"],
6629 )
6630 .expect_err("unknown smooth option should be rejected");
6631 assert!(
6632 err.contains("matern() does not accept option `lengt_scale`"),
6633 "error should name the invalid option, got: {err}"
6634 );
6635 assert!(
6636 err.contains("did you mean one of [length_scale]"),
6637 "error should suggest the closest valid option, got: {err}"
6638 );
6639 assert!(
6640 err.contains("Valid options: ["),
6641 "error should list valid option names, got: {err}"
6642 );
6643 }
6644
6645 #[test]
6646 fn tensor_k_accepts_square_bracket_per_margin_list() {
6647 let ds = continuous_dataset(
6648 &["y", "x", "z"],
6649 (0..40)
6650 .map(|i| {
6651 let x = i as f64 / 39.0;
6652 let z = ((i * 7) % 40) as f64 / 39.0;
6653 vec![x.sin() + z.cos(), x, z]
6654 })
6655 .collect(),
6656 );
6657
6658 assert_eq!(
6659 tensor_margin_basis_sizes(&ds, "y ~ te(x, z, k=[5, 6])"),
6660 vec![5, 6],
6661 "square-bracket k lists should materialize the requested per-margin values"
6662 );
6663 }
6664
6665 #[test]
6676 fn bare_doubly_cyclic_tensor_derives_period_from_data_range_1776() {
6677 let ds = continuous_dataset(
6678 &["y", "x", "z"],
6679 (0..40)
6680 .map(|i| {
6681 let x = i as f64 / 39.0;
6682 let z = ((i * 7) % 40) as f64 / 39.0;
6683 vec![x.sin() + z.cos(), x, z]
6684 })
6685 .collect(),
6686 );
6687
6688 let parsed = parse_formula("y ~ te(x, z, bs=c('cc','cc'))")
6689 .expect("parse doubly-cyclic tensor formula");
6690 let col_map = ds.column_map();
6691 let mut notes = Vec::new();
6692 let terms = build_termspec(
6695 &parsed.terms,
6696 &ds,
6697 &col_map,
6698 &mut notes,
6699 &ResourcePolicy::default_library(),
6700 )
6701 .expect(
6702 "bare cc-cc tensor must build via the data-range period fallback (#1776/#1752), \
6703 not hard-error on a missing explicit period",
6704 );
6705 let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
6706 panic!("expected tensor smooth");
6707 };
6708 assert_eq!(
6709 spec.marginalspecs.len(),
6710 2,
6711 "te(x, z) builds exactly two tensor margins"
6712 );
6713 for (axis, marginal) in spec.marginalspecs.iter().enumerate() {
6714 assert!(
6715 matches!(marginal.knotspec, BSplineKnotSpec::PeriodicUniform { .. }),
6716 "cyclic margin {axis} must build a periodic (wrapped) knotspec from the \
6717 data range, got {:?}",
6718 marginal.knotspec
6719 );
6720 }
6721 }
6722
6723 #[test]
6724 fn parse_cylinder_periodic_options_match_requested_forms() {
6725 let mut opts = BTreeMap::new();
6726 opts.insert("periodic".to_string(), "[0]".to_string());
6727 opts.insert("period".to_string(), "[2*pi, None]".to_string());
6728 let axes = parse_periodic_axes(&opts, 2).expect("axes");
6729 let periods = parse_periods(&opts, &axes).expect("periods");
6730 assert_eq!(axes, vec![true, false]);
6731 assert!((periods[0].unwrap() - 2.0 * std::f64::consts::PI).abs() < 1e-12);
6732 assert_eq!(periods[1], None);
6733
6734 let mut boundary_opts = BTreeMap::new();
6735 boundary_opts.insert(
6736 "boundary".to_string(),
6737 "['periodic', 'natural']".to_string(),
6738 );
6739 boundary_opts.insert("period".to_string(), "[2*pi, None]".to_string());
6740 let boundary_axes = parse_periodic_axes(&boundary_opts, 2).expect("boundary axes");
6741 let boundary_periods =
6742 parse_periods(&boundary_opts, &boundary_axes).expect("boundary periods");
6743 assert_eq!(boundary_axes, vec![true, false]);
6744 assert!((boundary_periods[0].unwrap() - 2.0 * std::f64::consts::PI).abs() < 1e-12);
6745 assert_eq!(boundary_periods[1], None);
6746
6747 let mut unicode_opts = BTreeMap::new();
6748 unicode_opts.insert("periodic".to_string(), "[0,1]".to_string());
6749 unicode_opts.insert("period".to_string(), "[2π, τ]".to_string());
6750 let unicode_axes = parse_periodic_axes(&unicode_opts, 2).expect("unicode axes");
6751 let unicode_periods = parse_periods(&unicode_opts, &unicode_axes).expect("unicode periods");
6752 assert_eq!(unicode_axes, vec![true, true]);
6753 assert!((unicode_periods[0].unwrap() - 2.0 * std::f64::consts::PI).abs() < 1e-12);
6754 assert!((unicode_periods[1].unwrap() - std::f64::consts::TAU).abs() < 1e-12);
6755 }
6756
6757 #[test]
6766 fn tensor_boundary_tokens_accept_clamped_open_reject_anchored() {
6767 fn boundary(raw: &str, dim: usize) -> Result<(), String> {
6768 let mut opts = BTreeMap::new();
6769 opts.insert("boundary".to_string(), raw.to_string());
6770 validate_tensor_boundary_tokens(&opts, dim)
6771 }
6772
6773 for raw in [
6776 "['periodic', 'clamped']",
6777 "['periodic', 'open']",
6778 "['cc', 'clamped']",
6779 "['clamped', 'natural']",
6780 "[Periodic, CLAMPED]",
6781 "c('cc', 'clamped')", ] {
6783 assert!(
6784 boundary(raw, 2).is_ok(),
6785 "boundary={raw:?} must be accepted (clamped/open/inert non-periodic markers)"
6786 );
6787 }
6788
6789 let mut bc_opts = BTreeMap::new();
6791 bc_opts.insert("bc".to_string(), "['periodic', 'clamped']".to_string());
6792 assert!(validate_tensor_boundary_tokens(&bc_opts, 2).is_ok());
6793
6794 let err = boundary("['periodic', 'anchored']", 2)
6798 .expect_err("anchored endpoint constraint must be rejected on a tensor margin");
6799 assert!(
6800 err.contains("anchored") && err.contains("not supported"),
6801 "rejection must name the offending token and be an unsupported-feature error: {err}"
6802 );
6803
6804 assert!(validate_tensor_boundary_tokens(&BTreeMap::new(), 2).is_ok());
6806 }
6807
6808 #[test]
6809 fn parse_single_axis_periodic_zero_as_axis_not_false() {
6810 let mut opts = BTreeMap::new();
6811 opts.insert("periodic".to_string(), "[0]".to_string());
6812 opts.insert("period".to_string(), "2*pi".to_string());
6813 opts.insert("origin".to_string(), "0".to_string());
6814 let axes = parse_periodic_axes(&opts, 1).expect("axes");
6815 let periods = parse_periods(&opts, &axes).expect("periods");
6816 let origins = parse_period_origins(&opts, &axes).expect("origins");
6817 assert_eq!(axes, vec![true]);
6818 assert!((periods[0].unwrap() - 2.0 * std::f64::consts::PI).abs() < 1e-12);
6819 assert_eq!(origins[0], Some(0.0));
6820 }
6821
6822 #[test]
6823 fn one_dimensional_bspline_accepts_boundary_periodic() {
6824 let ds = continuous_dataset(
6825 &["y", "theta"],
6826 (0..16)
6827 .map(|i| {
6828 let theta = std::f64::consts::TAU * i as f64 / 16.0;
6829 vec![theta.sin(), theta]
6830 })
6831 .collect(),
6832 );
6833 let parsed = parse_formula("y ~ s(theta, boundary=periodic, period=2*pi, origin=0, k=8)")
6834 .expect("parse");
6835 let col_map = ds.column_map();
6836 let mut notes = Vec::new();
6837 let terms = build_termspec(
6838 &parsed.terms,
6839 &ds,
6840 &col_map,
6841 &mut notes,
6842 &gam_runtime::resource::ResourcePolicy::default_library(),
6843 )
6844 .expect("periodic boundary should build");
6845 let SmoothBasisSpec::BSpline1D { spec, .. } = &terms.smooth_terms[0].basis else {
6846 panic!("expected 1D B-spline");
6847 };
6848 assert!(matches!(
6849 &spec.knotspec,
6850 BSplineKnotSpec::PeriodicUniform {
6851 data_range,
6852 num_basis: 8
6853 } if *data_range == (0.0, std::f64::consts::TAU)
6854 ));
6855 }
6856
6857 #[test]
6858 fn univariate_smooth_accepts_mgcv_cubic_regression_aliases() {
6859 let ds = continuous_dataset(
6860 &["y", "x"],
6861 (0..32)
6862 .map(|i| {
6863 let x = i as f64 / 31.0;
6864 vec![x * x, x]
6865 })
6866 .collect(),
6867 );
6868 let col_map = ds.column_map();
6869
6870 for selector in ["cr", "cs"] {
6871 let formula = format!("y ~ s(x, bs='{selector}')");
6872 let parsed = parse_formula(&formula).expect("parse cr/cs smooth");
6873 let mut notes = Vec::new();
6874 let terms = build_termspec(
6875 &parsed.terms,
6876 &ds,
6877 &col_map,
6878 &mut notes,
6879 &gam_runtime::resource::ResourcePolicy::default_library(),
6880 )
6881 .unwrap_or_else(|err| panic!("bs='{selector}' must build a 1-D smooth, got: {err:?}"));
6882 let SmoothBasisSpec::BSpline1D { spec, .. } = &terms.smooth_terms[0].basis else {
6883 panic!(
6884 "bs='{selector}' must lower to a BSpline1D; got {:?}",
6885 terms.smooth_terms[0].basis
6886 );
6887 };
6888 assert!(
6889 spec.double_penalty,
6890 "bs='{selector}' must recover its null space by default"
6891 );
6892
6893 let opt_out = format!("y ~ s(x, bs='{selector}', double_penalty=false)");
6894 let parsed = parse_formula(&opt_out).expect("parse explicit null-shrinkage opt-out");
6895 let mut notes = Vec::new();
6896 let terms = build_termspec(
6897 &parsed.terms,
6898 &ds,
6899 &col_map,
6900 &mut notes,
6901 &gam_runtime::resource::ResourcePolicy::default_library(),
6902 )
6903 .expect("explicit cr/cs opt-out should build");
6904 let SmoothBasisSpec::BSpline1D { spec, .. } = &terms.smooth_terms[0].basis else {
6905 panic!("bs='{selector}' must lower to a BSpline1D");
6906 };
6907 assert!(!spec.double_penalty, "explicit opt-out must be preserved");
6908 }
6909 }
6910
6911 #[test]
6912 fn non_intercept_linear_effects_default_to_mle_with_explicit_null_recovery() {
6913 let ds = continuous_dataset(
6914 &["y", "x", "z"],
6915 (0..24)
6916 .map(|i| {
6917 let x = i as f64 / 23.0;
6918 let z = 1.0 - x;
6919 vec![x - z, x, z]
6920 })
6921 .collect(),
6922 );
6923 let parsed = parse_formula("y ~ x + z + x:z").expect("parse linear defaults");
6924 let mut notes = Vec::new();
6925 let terms = build_termspec(
6926 &parsed.terms,
6927 &ds,
6928 &ds.column_map(),
6929 &mut notes,
6930 &gam_runtime::resource::ResourcePolicy::default_library(),
6931 )
6932 .expect("build linear defaults");
6933 assert!(!terms.linear_terms.is_empty());
6934 assert!(
6935 terms.linear_terms.iter().all(|term| !term.double_penalty),
6936 "ordinary parametric effects must be unpenalized by default: {:?}",
6937 terms
6938 .linear_terms
6939 .iter()
6940 .map(|term| (&term.name, term.double_penalty))
6941 .collect::<Vec<_>>()
6942 );
6943
6944 let bounded_parsed =
6948 parse_formula("y ~ bounded(z, min=-2, max=2)").expect("parse bounded defaults");
6949 let mut bounded_notes = Vec::new();
6950 let bounded_terms = build_termspec(
6951 &bounded_parsed.terms,
6952 &ds,
6953 &ds.column_map(),
6954 &mut bounded_notes,
6955 &gam_runtime::resource::ResourcePolicy::default_library(),
6956 )
6957 .expect("build bounded defaults");
6958 assert_eq!(bounded_terms.linear_terms.len(), 1);
6959 assert!(
6960 !bounded_terms.linear_terms[0].double_penalty,
6961 "bounded() must default double_penalty=false since it cannot combine with the interval transform"
6962 );
6963
6964 for formula in [
6965 "y ~ linear(x, double_penalty=true)",
6966 "y ~ linear(x:z, double_penalty=true)",
6967 ] {
6968 let parsed = parse_formula(formula).expect("parse explicit linear shrinkage");
6969 let mut notes = Vec::new();
6970 let terms = build_termspec(
6971 &parsed.terms,
6972 &ds,
6973 &ds.column_map(),
6974 &mut notes,
6975 &gam_runtime::resource::ResourcePolicy::default_library(),
6976 )
6977 .unwrap_or_else(|error| panic!("{formula} must build: {error}"));
6978 assert_eq!(terms.linear_terms.len(), 1, "{formula}");
6979 assert!(
6980 terms.linear_terms[0].double_penalty,
6981 "{formula} must preserve the explicit shrinkage opt-in"
6982 );
6983 }
6984
6985 assert!(
6986 parse_formula("y ~ linear(x, double_penalty=ture)").is_err(),
6987 "a misspelled opt-in must be rejected instead of silently using the default"
6988 );
6989 }
6990
6991 #[test]
6992 fn tensor_smooths_default_to_joint_null_recovery_with_explicit_opt_out() {
6993 let ds = continuous_dataset(
6994 &["y", "x", "z"],
6995 (0..36)
6996 .map(|i| {
6997 let x = i as f64 / 35.0;
6998 let z = ((i * 11) % 36) as f64 / 35.0;
6999 vec![x * z, x, z]
7000 })
7001 .collect(),
7002 );
7003 let col_map = ds.column_map();
7004 for constructor in ["te", "ti", "t2"] {
7005 for (option, expected) in [("", true), (", double_penalty=false", false)] {
7006 let formula = format!("y ~ {constructor}(x, z{option})");
7007 let parsed = parse_formula(&formula).expect("parse tensor default");
7008 let mut notes = Vec::new();
7009 let terms = build_termspec(
7010 &parsed.terms,
7011 &ds,
7012 &col_map,
7013 &mut notes,
7014 &gam_runtime::resource::ResourcePolicy::default_library(),
7015 )
7016 .unwrap_or_else(|error| panic!("{formula} must build: {error}"));
7017 let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis
7018 else {
7019 panic!("{formula} must lower to TensorBSpline");
7020 };
7021 assert_eq!(spec.double_penalty, expected, "{formula}");
7022 }
7023 }
7024 }
7025
7026 #[test]
7027 fn univariate_ps_small_k_degree_reduces_through_build() {
7028 let ds = continuous_dataset(
7037 &["y", "x"],
7038 (0..32)
7039 .map(|i| {
7040 let x = i as f64 / 31.0;
7041 vec![x * x, x]
7042 })
7043 .collect(),
7044 );
7045 let col_map = ds.column_map();
7046
7047 for formula in ["y ~ s(x, bs='ps', k=3)", "y ~ s(x, k=3)"] {
7048 let parsed = parse_formula(formula).expect("parse small-k ps/cr smooth");
7049 let mut notes = Vec::new();
7050 let terms = build_termspec(
7051 &parsed.terms,
7052 &ds,
7053 &col_map,
7054 &mut notes,
7055 &gam_runtime::resource::ResourcePolicy::default_library(),
7056 )
7057 .unwrap_or_else(|err| {
7058 panic!("`{formula}` must degree-reduce, not error; got: {err:?}")
7059 });
7060 let SmoothBasisSpec::BSpline1D { spec, .. } = &terms.smooth_terms[0].basis else {
7061 panic!(
7062 "`{formula}` must lower to a BSpline1D; got {:?}",
7063 terms.smooth_terms[0].basis
7064 );
7065 };
7066 assert_eq!(
7067 spec.degree, 2,
7068 "`{formula}` must drop the cubic default to a quadratic basis"
7069 );
7070 let num_internal = match &spec.knotspec {
7071 BSplineKnotSpec::Generate {
7072 num_internal_knots, ..
7073 } => *num_internal_knots,
7074 BSplineKnotSpec::Automatic {
7075 num_internal_knots: Some(n),
7076 ..
7077 } => *n,
7078 other => panic!("`{formula}` unexpected knotspec: {other:?}"),
7079 };
7080 assert_eq!(
7081 num_internal, 0,
7082 "`{formula}` must have zero internal knots (num_basis = k = 3)"
7083 );
7084 assert!(
7086 spec.penalty_order >= 1 && spec.penalty_order <= spec.degree,
7087 "`{formula}` penalty_order {} must satisfy 1 <= order <= degree={}",
7088 spec.penalty_order,
7089 spec.degree
7090 );
7091 }
7092 }
7093
7094 #[test]
7095 fn formula_shape_constraint_round_trips_and_rejects_bogus() {
7096 let ds = continuous_dataset(
7097 &["y", "x"],
7098 (0..32)
7099 .map(|i| {
7100 let x = i as f64 / 31.0;
7101 vec![x * x, x]
7102 })
7103 .collect(),
7104 );
7105 let col_map = ds.column_map();
7106
7107 let parsed =
7108 parse_formula("y ~ s(x, shape=monotone_increasing)").expect("parse monotone smooth");
7109 let mut notes = Vec::new();
7110 let terms = build_termspec(
7111 &parsed.terms,
7112 &ds,
7113 &col_map,
7114 &mut notes,
7115 &gam_runtime::resource::ResourcePolicy::default_library(),
7116 )
7117 .expect("monotone smooth should build");
7118 assert_eq!(
7119 terms.smooth_terms[0].shape,
7120 ShapeConstraint::MonotoneIncreasing
7121 );
7122
7123 let parsed_bad = parse_formula("y ~ s(x, shape=bogus)").expect("parse bogus shape");
7124 let mut notes_bad = Vec::new();
7125 let err = build_termspec(
7126 &parsed_bad.terms,
7127 &ds,
7128 &col_map,
7129 &mut notes_bad,
7130 &gam_runtime::resource::ResourcePolicy::default_library(),
7131 )
7132 .expect_err("bogus shape must error");
7133 assert!(
7134 format!("{err:?}").contains("unknown shape constraint"),
7135 "got: {err:?}"
7136 );
7137 }
7138
7139 #[test]
7140 fn default_sphere_smooth_uses_spherical_farthest_point_centers() {
7141 let ds = continuous_dataset(
7142 &["y", "lat", "lon"],
7143 (0..24)
7144 .map(|i| {
7145 let t = i as f64 / 24.0;
7146 let lat = -60.0 + 120.0 * t;
7147 let lon = -180.0 + 360.0 * ((7 * i) % 24) as f64 / 24.0;
7148 vec![lat.to_radians().sin(), lat, lon]
7149 })
7150 .collect(),
7151 );
7152 let parsed = parse_formula("y ~ sphere(lat, lon)").expect("parse");
7153 let col_map = ds.column_map();
7154 let mut notes = Vec::new();
7155 let terms = build_termspec(
7156 &parsed.terms,
7157 &ds,
7158 &col_map,
7159 &mut notes,
7160 &gam_runtime::resource::ResourcePolicy::default_library(),
7161 )
7162 .expect("build sphere termspec");
7163 let SmoothBasisSpec::Sphere { spec, .. } = &terms.smooth_terms[0].basis else {
7164 panic!("expected sphere term");
7165 };
7166 assert!(matches!(
7167 spec.center_strategy,
7168 CenterStrategy::FarthestPoint { .. }
7169 ));
7170 }
7171
7172 #[test]
7173 fn one_dimensional_duchon_defaults_to_scale_free_length_scale() {
7174 let ds = continuous_dataset(
7175 &["y", "x"],
7176 (0..32)
7177 .map(|i| {
7178 let x = i as f64 / 31.0;
7179 vec![(std::f64::consts::TAU * x).sin(), x]
7180 })
7181 .collect(),
7182 );
7183 let parsed = parse_formula("y ~ duchon(x)").expect("parse");
7184 let col_map = ds.column_map();
7185 let mut notes = Vec::new();
7186 let terms = build_termspec(
7187 &parsed.terms,
7188 &ds,
7189 &col_map,
7190 &mut notes,
7191 &gam_runtime::resource::ResourcePolicy::default_library(),
7192 )
7193 .expect("build default duchon termspec");
7194 let SmoothBasisSpec::Duchon { spec, .. } = &terms.smooth_terms[0].basis else {
7195 panic!("expected Duchon term");
7196 };
7197 assert_eq!(spec.length_scale, None);
7198 assert!(matches!(
7199 spec.center_strategy,
7200 CenterStrategy::Auto(ref inner)
7201 if matches!(
7202 inner.as_ref(),
7203 CenterStrategy::UniformGrid { .. }
7204 )
7205 ));
7206 }
7207
7208 #[test]
7209 fn formula_duchon_default_does_not_enable_collocation_operators() {
7210 let ds = continuous_dataset(
7211 &["y", "x", "z"],
7212 (0..40)
7213 .map(|i| {
7214 let x = (i as f64 / 39.0).fract();
7215 let z = ((7 * i) as f64 / 39.0).fract();
7216 vec![x + z, x, z]
7217 })
7218 .collect(),
7219 );
7220 let parsed = parse_formula("y ~ duchon(x, z)").expect("parse");
7221 let col_map = ds.column_map();
7222 let mut notes = Vec::new();
7223 let terms = build_termspec(
7224 &parsed.terms,
7225 &ds,
7226 &col_map,
7227 &mut notes,
7228 &gam_runtime::resource::ResourcePolicy::default_library(),
7229 )
7230 .expect("build default 2D duchon termspec");
7231 let SmoothBasisSpec::Duchon { spec, .. } = &terms.smooth_terms[0].basis else {
7232 panic!("expected Duchon term");
7233 };
7234 assert!(matches!(
7235 spec.operator_penalties.mass,
7236 OperatorPenaltySpec::Disabled
7237 ));
7238 assert!(matches!(
7239 spec.operator_penalties.tension,
7240 OperatorPenaltySpec::Disabled
7241 ));
7242 assert!(matches!(
7243 spec.operator_penalties.stiffness,
7244 OperatorPenaltySpec::Disabled
7245 ));
7246 }
7247
7248 #[test]
7249 fn one_dimensional_duchon_length_scale_opts_into_hybrid_mode() {
7250 let ds = continuous_dataset(
7251 &["y", "x"],
7252 (0..32)
7253 .map(|i| {
7254 let x = i as f64 / 31.0;
7255 vec![(std::f64::consts::TAU * x).sin(), x]
7256 })
7257 .collect(),
7258 );
7259 let parsed = parse_formula("y ~ duchon(x, length_scale=0.25)").expect("parse");
7260 let col_map = ds.column_map();
7261 let mut notes = Vec::new();
7262 let terms = build_termspec(
7263 &parsed.terms,
7264 &ds,
7265 &col_map,
7266 &mut notes,
7267 &gam_runtime::resource::ResourcePolicy::default_library(),
7268 )
7269 .expect("build hybrid duchon termspec");
7270 let SmoothBasisSpec::Duchon { spec, .. } = &terms.smooth_terms[0].basis else {
7271 panic!("expected Duchon term");
7272 };
7273 assert_eq!(spec.length_scale, Some(0.25));
7274 }
7275
7276 #[test]
7277 fn multidimensional_duchon_default_uses_low_rank_mgcv_sized_basis() {
7278 let ds = continuous_dataset(
7279 &["y", "x1", "x2"],
7280 (0..500)
7281 .map(|i| {
7282 let x1 = 2.0 * (i as f64 / 499.0) - 1.0;
7283 let x2 = (((37 * i) % 500) as f64 / 499.0) * 2.0 - 1.0;
7284 vec![(2.0 * x1).sin() + (1.5 * x2).cos(), x1, x2]
7285 })
7286 .collect(),
7287 );
7288 let parsed = parse_formula("y ~ duchon(x1, x2)").expect("parse");
7289 let col_map = ds.column_map();
7290 let mut notes = Vec::new();
7291 let terms = build_termspec(
7292 &parsed.terms,
7293 &ds,
7294 &col_map,
7295 &mut notes,
7296 &gam_runtime::resource::ResourcePolicy::default_library(),
7297 )
7298 .expect("build default 2D duchon termspec");
7299 let SmoothBasisSpec::Duchon { spec, .. } = &terms.smooth_terms[0].basis else {
7300 panic!("expected Duchon term");
7301 };
7302 let CenterStrategy::Auto(inner) = &spec.center_strategy else {
7303 panic!("expected auto center strategy");
7304 };
7305 assert!(matches!(
7306 inner.as_ref(),
7307 CenterStrategy::FarthestPoint { num_centers: 30 }
7308 ));
7309 }
7310
7311 #[test]
7312 fn spectral_duchon_reproduces_fixed_seed_uniform_landmarks() {
7313 let ds = continuous_dataset(
7314 &["y", "x1", "x2", "x3", "x4"],
7315 (0..64)
7316 .map(|i| {
7317 let x = i as f64 / 63.0;
7318 vec![
7319 x.sin(),
7320 x,
7321 (3.0 * x).sin(),
7322 (5.0 * x).cos(),
7323 (7.0 * x).sin(),
7324 ]
7325 })
7326 .collect(),
7327 );
7328 let parsed = parse_formula("y ~ duchon(x1, x2, x3, x4, rank=6, order=0)").expect("parse");
7329 let col_map = ds.column_map();
7330 let mut notes = Vec::new();
7331 let terms = build_termspec(
7332 &parsed.terms,
7333 &ds,
7334 &col_map,
7335 &mut notes,
7336 &gam_runtime::resource::ResourcePolicy::default_library(),
7337 )
7338 .expect("build spectral Duchon termspec");
7339 let SmoothBasisSpec::Duchon { spec, .. } = &terms.smooth_terms[0].basis else {
7340 panic!("expected Duchon term");
7341 };
7342 let CenterStrategy::DuchonSpectral { knots, basis } = &spec.center_strategy else {
7343 panic!("expected spectral center strategy");
7344 };
7345 assert_eq!(basis.rank(), 6);
7346 let CenterStrategy::UserProvided(centers) = knots.as_ref() else {
7347 panic!("expected frozen sampled centers");
7348 };
7349 assert_eq!(centers.dim(), (64, 4));
7350 }
7351
7352 #[test]
7353 fn parse_matern_nu_accepts_equivalent_half_integer_forms() {
7354 let cases = [
7355 ("1/2", MaternNu::Half),
7356 (" 1 / 2 ", MaternNu::Half),
7357 (".5", MaternNu::Half),
7358 ("0.50", MaternNu::Half),
7359 ("half", MaternNu::Half),
7360 ("3 / 2", MaternNu::ThreeHalves),
7361 ("1.50", MaternNu::ThreeHalves),
7362 ("5 / 2", MaternNu::FiveHalves),
7363 ("2.500000000000", MaternNu::FiveHalves),
7364 ("7 / 2", MaternNu::SevenHalves),
7365 ("3.50", MaternNu::SevenHalves),
7366 ("9 / 2", MaternNu::NineHalves),
7367 ("4.50", MaternNu::NineHalves),
7368 ];
7369 for (raw, expected) in cases {
7370 let parsed = parse_matern_nu(raw).expect(raw);
7371 assert!(
7372 matches!(
7373 (parsed, expected),
7374 (MaternNu::Half, MaternNu::Half)
7375 | (MaternNu::ThreeHalves, MaternNu::ThreeHalves)
7376 | (MaternNu::FiveHalves, MaternNu::FiveHalves)
7377 | (MaternNu::SevenHalves, MaternNu::SevenHalves)
7378 | (MaternNu::NineHalves, MaternNu::NineHalves)
7379 ),
7380 "parsed {raw:?} as {parsed:?}, expected {expected:?}"
7381 );
7382 }
7383 }
7384
7385 #[test]
7386 fn parse_matern_nu_rejects_unsupported_or_invalid_values() {
7387 for raw in ["1", "2", "11/2", "1/0", "nan", "fast"] {
7388 let err = parse_matern_nu(raw).expect_err(raw);
7389 assert!(
7390 err.contains("supported half-integer values"),
7391 "unexpected error for {raw:?}: {err}"
7392 );
7393 }
7394 }
7395
7396 #[test]
7397 fn parse_ps_k_promotes_underexpressive_cubic_basis() {
7398 let mut opts = BTreeMap::new();
7399 opts.insert("k".to_string(), "4".to_string());
7400 let (internal, inferred, eff_degree) = parse_ps_internal_knots(&opts, 3, 20).expect("k=4");
7401 assert_eq!(internal, 2);
7402 assert_eq!(eff_degree, 3);
7403 assert!(!inferred);
7404
7405 opts.insert("k".to_string(), "6".to_string());
7406 let (internal, inferred, eff_degree) = parse_ps_internal_knots(&opts, 3, 20).expect("k=6");
7407 assert_eq!(internal, 2);
7408 assert_eq!(eff_degree, 3);
7409 assert!(!inferred);
7410
7411 opts.insert("k".to_string(), "10".to_string());
7412 let (internal, inferred, eff_degree) = parse_ps_internal_knots(&opts, 3, 20).expect("k=10");
7413 assert_eq!(internal, 6);
7414 assert_eq!(eff_degree, 3);
7415 assert!(!inferred);
7416 }
7417
7418 #[test]
7419 fn parse_ps_internal_knots_drops_degree_for_small_k() {
7420 let mut opts = BTreeMap::new();
7425 opts.insert("k".to_string(), "3".to_string());
7426 let (internal, inferred, eff_degree) = parse_ps_internal_knots(&opts, 3, 20).expect("k=3");
7427 assert_eq!(eff_degree, 2);
7428 assert_eq!(internal, 0);
7429 assert!(!inferred);
7430
7431 opts.insert("k".to_string(), "2".to_string());
7434 let (internal, inferred, eff_degree) = parse_ps_internal_knots(&opts, 3, 20).expect("k=2");
7435 assert_eq!(eff_degree, 1);
7436 assert_eq!(internal, 0);
7437 assert!(!inferred);
7438
7439 opts.insert("k".to_string(), "1".to_string());
7443 let err = parse_ps_internal_knots(&opts, 3, 20)
7444 .expect_err("k=1 is below the irreducible spline floor");
7445 assert!(err.contains("requires k >= 2"), "unexpected error: {err}");
7446
7447 opts.insert("k".to_string(), "4".to_string());
7450 let (internal, inferred, eff_degree) = parse_ps_internal_knots(&opts, 3, 20).expect("k=4");
7451 assert_eq!(eff_degree, 3);
7452 assert_eq!(internal, 2);
7453 assert!(!inferred);
7454 }
7455
7456 #[test]
7457 fn factor_smooth_marginal_degree_reduces_for_small_k() {
7458 let ds = factor_dataset();
7459 let col_map = ds.column_map();
7460
7461 for (k, expected_degree) in [(3usize, 2usize), (2usize, 1usize)] {
7462 let parsed =
7463 parse_formula(&format!("y ~ s(x, g, bs=fs, k={k})")).expect("parse factor smooth");
7464 let mut notes = Vec::new();
7465 let terms = build_termspec(
7466 &parsed.terms,
7467 &ds,
7468 &col_map,
7469 &mut notes,
7470 &gam_runtime::resource::ResourcePolicy::default_library(),
7471 )
7472 .unwrap_or_else(|err| panic!("fs k={k} should degree-reduce, got: {err:?}"));
7473 let SmoothBasisSpec::FactorSmooth { spec } = &terms.smooth_terms[0].basis else {
7474 panic!(
7475 "expected factor smooth, got {:?}",
7476 terms.smooth_terms[0].basis
7477 );
7478 };
7479 assert_eq!(spec.marginal.degree, expected_degree);
7480 assert!(
7481 spec.marginal.penalty_order <= spec.marginal.degree,
7482 "penalty_order {} must be clamped to degree {}",
7483 spec.marginal.penalty_order,
7484 spec.marginal.degree
7485 );
7486 let basis_size = match spec.marginal.knotspec {
7487 BSplineKnotSpec::Generate {
7488 num_internal_knots, ..
7489 } => num_internal_knots + spec.marginal.degree + 1,
7490 BSplineKnotSpec::Automatic {
7491 num_internal_knots: Some(num_internal_knots),
7492 ..
7493 } => num_internal_knots + spec.marginal.degree + 1,
7494 ref other => panic!("unexpected factor-smooth knotspec: {other:?}"),
7495 };
7496 assert_eq!(basis_size, k);
7497 }
7498 }
7499
7500 fn ternary_factor_dataset() -> Dataset {
7503 let rows = (0..120)
7504 .map(|i| {
7505 let x = (i % 3) as f64;
7506 let g = (i % 2) as f64;
7507 vec![x + g, x, g]
7508 })
7509 .collect::<Vec<_>>();
7510 Dataset {
7511 headers: vec!["y".into(), "x".into(), "g".into()],
7512 values: Array2::from_shape_vec(
7513 (rows.len(), 3),
7514 rows.into_iter().flat_map(|row| row.into_iter()).collect(),
7515 )
7516 .expect("rectangular ternary factor test data"),
7517 schema: DataSchema {
7518 columns: vec![
7519 SchemaColumn {
7520 name: "y".into(),
7521 kind: ColumnKindTag::Continuous,
7522 levels: vec![],
7523 },
7524 SchemaColumn {
7525 name: "x".into(),
7526 kind: ColumnKindTag::Continuous,
7527 levels: vec![],
7528 },
7529 SchemaColumn {
7530 name: "g".into(),
7531 kind: ColumnKindTag::Categorical,
7532 levels: vec!["a".into(), "b".into()],
7533 },
7534 ],
7535 },
7536 column_kinds: vec![
7537 ColumnKindTag::Continuous,
7538 ColumnKindTag::Continuous,
7539 ColumnKindTag::Categorical,
7540 ],
7541 }
7542 }
7543
7544 #[test]
7545 fn univariate_cr_smooth_caps_knots_to_data_support() {
7546 let ds = continuous_dataset(
7552 &["y", "x"],
7553 (0..90)
7554 .map(|i| vec![(i % 3) as f64, (i % 3) as f64])
7555 .collect(),
7556 );
7557 let col_map = ds.column_map();
7558 let parsed = parse_formula("y ~ s(x, bs=cr, k=10)").expect("parse cr smooth");
7559 let mut notes = Vec::new();
7560 let terms = build_termspec(
7561 &parsed.terms,
7562 &ds,
7563 &col_map,
7564 &mut notes,
7565 &gam_runtime::resource::ResourcePolicy::default_library(),
7566 )
7567 .expect("cr k=10 must cap to data support instead of erroring");
7568 let SmoothBasisSpec::BSpline1D { spec, .. } = &terms.smooth_terms[0].basis else {
7569 panic!("expected BSpline1D for s(x, bs=cr)");
7570 };
7571 let BSplineKnotSpec::NaturalCubicRegression { knots } = &spec.knotspec else {
7572 panic!("expected cr knotspec, got {:?}", spec.knotspec);
7573 };
7574 assert_eq!(knots.len(), 3, "cr basis not capped to 3 distinct values");
7576 assert_eq!(knots.as_slice().unwrap(), &[0.0, 1.0, 2.0]);
7577 assert!(
7579 notes.iter().any(|n| n.contains("data-support cap")),
7580 "cap not reported in inference notes: {notes:?}"
7581 );
7582 }
7583
7584 #[test]
7585 fn univariate_cr_smooth_binary_covariate_degrades_to_bspline() {
7586 let ds = continuous_dataset(
7590 &["y", "x"],
7591 (0..80)
7592 .map(|i| vec![(i % 2) as f64, (i % 2) as f64])
7593 .collect(),
7594 );
7595 let col_map = ds.column_map();
7596 let parsed = parse_formula("y ~ s(x, bs=cr, k=10)").expect("parse cr smooth");
7597 let mut notes = Vec::new();
7598 let terms = build_termspec(
7599 &parsed.terms,
7600 &ds,
7601 &col_map,
7602 &mut notes,
7603 &gam_runtime::resource::ResourcePolicy::default_library(),
7604 )
7605 .expect("binary cr must degrade to B-spline instead of erroring");
7606 let SmoothBasisSpec::BSpline1D { spec, .. } = &terms.smooth_terms[0].basis else {
7607 panic!("expected BSpline1D for s(x, bs=cr)");
7608 };
7609 assert!(
7610 !matches!(
7611 spec.knotspec,
7612 BSplineKnotSpec::NaturalCubicRegression { .. }
7613 ),
7614 "binary covariate must NOT build a cr basis, got {:?}",
7615 spec.knotspec
7616 );
7617 assert!(
7618 notes
7619 .iter()
7620 .any(|n| n.contains("Degraded to the linear B-spline")),
7621 "degradation not reported in inference notes: {notes:?}"
7622 );
7623 }
7624
7625 #[test]
7630 fn one_dimensional_identifiability_option_is_parsed_and_validated() {
7631 let mut options = BTreeMap::new();
7632
7633 assert!(matches!(
7635 parse_bspline_identifiability(&options).expect("absent option parses"),
7636 None
7637 ));
7638 assert!(matches!(
7639 resolve_bspline_identifiability(
7640 &options,
7641 BSplineIdentifiability::None,
7642 BSplineIdentifiabilityContext::default(),
7643 )
7644 .expect("absent option keeps the structural default"),
7645 BSplineIdentifiability::None
7646 ));
7647
7648 for token in ["none", "None", " NONE "] {
7649 options.insert("identifiability".to_string(), token.to_string());
7650 assert!(
7651 matches!(
7652 parse_bspline_identifiability(&options).expect("none parses"),
7653 Some(BSplineIdentifiability::None)
7654 ),
7655 "token {token:?} should select the unconstrained policy"
7656 );
7657 }
7658 for token in [
7659 "sum_tozero",
7660 "sum-to-zero",
7661 "sumtozero",
7662 "centered",
7663 "center_sum_tozero",
7664 "center-sum-to-zero",
7665 ] {
7666 options.insert("identifiability".to_string(), token.to_string());
7667 assert!(
7668 matches!(
7669 parse_bspline_identifiability(&options).expect("sum-to-zero parses"),
7670 Some(BSplineIdentifiability::WeightedSumToZero { weights: None })
7671 ),
7672 "token {token:?} should select sum-to-zero centering"
7673 );
7674 }
7675 for token in [
7676 "linear",
7677 "remove_linear_trend",
7678 "remove-linear-trend",
7679 "center_linear_orthogonal",
7680 ] {
7681 options.insert("identifiability".to_string(), token.to_string());
7682 assert!(
7683 matches!(
7684 parse_bspline_identifiability(&options).expect("linear parses"),
7685 Some(BSplineIdentifiability::RemoveLinearTrend)
7686 ),
7687 "token {token:?} should select the constant+linear removal"
7688 );
7689 }
7690
7691 options.insert("identifiability".to_string(), "none".to_string());
7693 assert!(matches!(
7694 resolve_bspline_identifiability(
7695 &options,
7696 BSplineIdentifiability::default(),
7697 BSplineIdentifiabilityContext::default(),
7698 )
7699 .expect("explicit none overrides the centering default"),
7700 BSplineIdentifiability::None
7701 ));
7702 options.insert("identifiability".to_string(), "sum_tozero".to_string());
7703 assert!(matches!(
7704 resolve_bspline_identifiability(
7705 &options,
7706 BSplineIdentifiability::None,
7707 BSplineIdentifiabilityContext::default(),
7708 )
7709 .expect("explicit sum_tozero overrides an unconstrained default"),
7710 BSplineIdentifiability::WeightedSumToZero { weights: None }
7711 ));
7712
7713 for token in ["frozen", "orthogonal"] {
7715 options.insert("identifiability".to_string(), token.to_string());
7716 let err = parse_bspline_identifiability(&options)
7717 .expect_err("internal-only policy must be refused");
7718 assert!(
7719 err.contains("internal-only"),
7720 "token {token:?} should be refused as internal-only, got: {err}"
7721 );
7722 }
7723
7724 options.insert("identifiability".to_string(), "totally_bogus".to_string());
7727 let err = parse_bspline_identifiability(&options)
7728 .expect_err("an unknown identifiability token must be refused");
7729 assert!(
7730 err.contains("totally_bogus") && err.contains("none, sum_tozero, linear"),
7731 "unknown-token error should name the token and the vocabulary, got: {err}"
7732 );
7733
7734 options.insert("identifiability".to_string(), "sum_tozero".to_string());
7736 let err = resolve_bspline_identifiability(
7737 &options,
7738 BSplineIdentifiability::None,
7739 BSplineIdentifiabilityContext {
7740 has_anchor: true,
7741 ..Default::default()
7742 },
7743 )
7744 .expect_err("anchor + centering is over-constrained");
7745 assert!(
7746 err.contains("anchored endpoint"),
7747 "anchor conflict should explain itself, got: {err}"
7748 );
7749 options.insert("identifiability".to_string(), "none".to_string());
7751 assert!(matches!(
7752 resolve_bspline_identifiability(
7753 &options,
7754 BSplineIdentifiability::None,
7755 BSplineIdentifiabilityContext {
7756 has_anchor: true,
7757 ..Default::default()
7758 },
7759 )
7760 .expect("anchor + none agrees with the structural default"),
7761 BSplineIdentifiability::None
7762 ));
7763
7764 options.insert("identifiability".to_string(), "linear".to_string());
7766 let err = resolve_bspline_identifiability(
7767 &options,
7768 BSplineIdentifiability::default(),
7769 BSplineIdentifiabilityContext {
7770 periodic: true,
7771 ..Default::default()
7772 },
7773 )
7774 .expect_err("a linear trend is not in the span of a cyclic basis");
7775 assert!(err.contains("periodic"), "got: {err}");
7776 let err = resolve_bspline_identifiability(
7777 &options,
7778 BSplineIdentifiability::default(),
7779 BSplineIdentifiabilityContext {
7780 natural_cubic_regression: true,
7781 ..Default::default()
7782 },
7783 )
7784 .expect_err("cr carries no Greville chart");
7785 assert!(err.contains("cr"), "got: {err}");
7786 }
7787
7788 #[test]
7792 fn one_dimensional_identifiability_option_reaches_the_built_spec() {
7793 let ds = continuous_dataset(
7794 &["y", "x"],
7795 (0..120)
7796 .map(|i| {
7797 let x = i as f64 / 119.0;
7798 vec![x.sin(), x]
7799 })
7800 .collect(),
7801 );
7802 let col_map = ds.column_map();
7803 let policy = gam_runtime::resource::ResourcePolicy::default_library();
7804
7805 let built = |formula: &str| -> BSplineIdentifiability {
7806 let parsed = parse_formula(formula).expect("parse");
7807 let mut notes = Vec::new();
7808 let terms = build_termspec(&parsed.terms, &ds, &col_map, &mut notes, &policy)
7809 .unwrap_or_else(|e| panic!("{formula} should build: {e}"));
7810 let SmoothBasisSpec::BSpline1D { spec, .. } = &terms.smooth_terms[0].basis else {
7811 panic!("expected BSpline1D for {formula}");
7812 };
7813 spec.identifiability.clone()
7814 };
7815
7816 assert!(matches!(
7817 built("y ~ s(x, k=8)"),
7818 BSplineIdentifiability::WeightedSumToZero { .. }
7819 ));
7820 assert!(matches!(
7821 built("y ~ s(x, k=8, identifiability='none')"),
7822 BSplineIdentifiability::None
7823 ));
7824 assert!(matches!(
7825 built("y ~ s(x, k=8, identifiability='linear')"),
7826 BSplineIdentifiability::RemoveLinearTrend
7827 ));
7828 assert!(matches!(
7829 built("y ~ cyclic(x, k=8, period=1)"),
7830 BSplineIdentifiability::WeightedSumToZero { .. }
7831 ));
7832 assert!(matches!(
7833 built("y ~ cyclic(x, k=8, period=1, identifiability='none')"),
7834 BSplineIdentifiability::None
7835 ));
7836
7837 for formula in [
7838 "y ~ s(x, k=8, identifiability='totally_bogus')",
7839 "y ~ cyclic(x, k=8, period=1, identifiability='totally_bogus')",
7840 "y ~ cyclic(x, k=8, period=1, identifiability='linear')",
7841 "y ~ s(x, k=8, bc_left=anchored, anchor_left=0, identifiability='sum_tozero')",
7842 ] {
7843 let parsed = parse_formula(formula).expect("parse");
7844 let mut notes = Vec::new();
7845 build_termspec(&parsed.terms, &ds, &col_map, &mut notes, &policy)
7846 .expect_err(&format!("{formula} must be refused, not silently accepted"));
7847 }
7848 }
7849
7850 #[test]
7853 fn a_declared_period_makes_its_axis_periodic() {
7854 let opts = |pairs: &[(&str, &str)]| -> BTreeMap<String, String> {
7855 pairs
7856 .iter()
7857 .map(|(k, v)| (k.to_string(), v.to_string()))
7858 .collect()
7859 };
7860
7861 assert_eq!(
7864 parse_periodic_axes(&opts(&[("period", "24")]), 1).expect("period=24"),
7865 vec![true]
7866 );
7867 assert_eq!(
7868 parse_periodic_axes(&opts(&[("periods", "24")]), 1).expect("periods=24"),
7869 vec![true]
7870 );
7871 assert_eq!(
7872 parse_periodic_axes(&opts(&[("period_start", "0"), ("period_end", "24")]), 1)
7873 .expect("endpoint form"),
7874 vec![true]
7875 );
7876 assert_eq!(
7878 parse_periodic_axes(&opts(&[("k", "8")]), 1).expect("no declaration"),
7879 vec![false]
7880 );
7881
7882 assert_eq!(
7884 parse_tensor_periodic_axes(&opts(&[("periods", "[2*pi, None]")]), 2)
7885 .expect("per-margin periods"),
7886 vec![true, false]
7887 );
7888 assert_eq!(
7889 parse_tensor_periodic_axes(&opts(&[("period", "[None, 24]")]), 2)
7890 .expect("per-margin period"),
7891 vec![false, true]
7892 );
7893 assert_eq!(
7896 parse_tensor_periodic_axes(&opts(&[("period", "24")]), 2).expect("scalar on 2-D"),
7897 vec![false, false]
7898 );
7899 assert_eq!(
7901 parse_tensor_periodic_axes(&opts(&[("bc", "periodic")]), 2).expect("scalar bc"),
7902 vec![true, true]
7903 );
7904
7905 let err = parse_periodic_axes(&opts(&[("periodic", "false"), ("period", "24")]), 1)
7908 .expect_err("periodic=false + period= is a contradiction");
7909 assert!(err.contains("denies the periodicity"), "got: {err}");
7910
7911 let err = reject_unconsumable_period_declaration(
7913 "tensor",
7914 &opts(&[("period", "24")]),
7915 &[false, false],
7916 )
7917 .expect_err("a scalar period on a 2-margin tensor names no margin");
7918 assert!(err.contains("does not say which"), "got: {err}");
7919 let err = reject_unconsumable_period_declaration(
7920 "bspline",
7921 &opts(&[("origin", "0")]),
7922 &[false],
7923 )
7924 .expect_err("an origin with no period is unconsumable");
7925 assert!(err.contains("declares no period"), "got: {err}");
7926 reject_unconsumable_period_declaration(
7928 "bspline",
7929 &opts(&[("period", "24"), ("origin", "0")]),
7930 &[true],
7931 )
7932 .expect("a periodic axis consumes its own declaration");
7933 }
7934
7935 #[test]
7939 fn tensor_per_axis_integer_options_parse_scalar_and_list_forms() {
7940 let mut options = BTreeMap::new();
7941 assert_eq!(
7942 parse_tensor_per_axis_usize(&options, "degree", 2).expect("absent"),
7943 vec![None, None]
7944 );
7945
7946 options.insert("degree".to_string(), "2".to_string());
7947 assert_eq!(
7948 parse_tensor_per_axis_usize(&options, "degree", 3).expect("scalar broadcasts"),
7949 vec![Some(2), Some(2), Some(2)]
7950 );
7951
7952 for spelling in ["[1, 3]", "c(1, 3)", "(1,3)"] {
7953 options.insert("degree".to_string(), spelling.to_string());
7954 assert_eq!(
7955 parse_tensor_per_axis_usize(&options, "degree", 2)
7956 .unwrap_or_else(|e| panic!("{spelling}: {e}")),
7957 vec![Some(1), Some(3)],
7958 "spelling {spelling} should parse per margin"
7959 );
7960 }
7961
7962 options.insert("degree".to_string(), "[1, none]".to_string());
7963 assert_eq!(
7964 parse_tensor_per_axis_usize(&options, "degree", 2).expect("none keeps the default"),
7965 vec![Some(1), None]
7966 );
7967
7968 options.insert("degree".to_string(), "[1, 2, 3]".to_string());
7969 let err = parse_tensor_per_axis_usize(&options, "degree", 2)
7970 .expect_err("a length mismatch must be refused");
7971 assert!(err.contains("3 entries") && err.contains("2 margins"), "got: {err}");
7972
7973 options.insert("degree".to_string(), "[1, banana]".to_string());
7974 let err = parse_tensor_per_axis_usize(&options, "degree", 2)
7975 .expect_err("a non-integer entry must be refused");
7976 assert!(err.contains("banana"), "got: {err}");
7977
7978 let mut placement = BTreeMap::new();
7979 assert!(
7980 explicit_knot_placement(&placement)
7981 .expect("absent")
7982 .is_none()
7983 );
7984 placement.insert("knot_placement".to_string(), "uniform".to_string());
7985 assert_eq!(
7986 explicit_knot_placement(&placement).expect("explicit uniform"),
7987 Some(crate::basis::BSplineKnotPlacement::Uniform)
7988 );
7989 }
7990
7991 #[test]
7994 fn tensor_margin_leaves_cr_only_when_the_request_needs_a_bspline() {
7995 let ds = continuous_dataset(
7996 &["y", "x", "z"],
7997 (0..200)
7998 .map(|i| {
7999 let x = (i % 20) as f64 / 19.0;
8000 let z = (i / 20) as f64 / 9.0;
8001 vec![x + z, x, z]
8002 })
8003 .collect(),
8004 );
8005 let col_map = ds.column_map();
8006 let policy = gam_runtime::resource::ResourcePolicy::default_library();
8007 let margins = |formula: &str| -> Vec<BSplineKnotSpec> {
8008 let parsed = parse_formula(formula).expect("parse");
8009 let mut notes = Vec::new();
8010 let terms = build_termspec(&parsed.terms, &ds, &col_map, &mut notes, &policy)
8011 .unwrap_or_else(|e| panic!("{formula} should build: {e}"));
8012 let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
8013 panic!("expected a tensor spec for {formula}");
8014 };
8015 spec.marginalspecs
8016 .iter()
8017 .map(|m| m.knotspec.clone())
8018 .collect()
8019 };
8020 let is_cr = |k: &BSplineKnotSpec| {
8021 matches!(k, BSplineKnotSpec::NaturalCubicRegression { .. })
8022 };
8023
8024 for formula in [
8027 "y ~ te(x, z, k=5)",
8028 "y ~ te(x, z, k=5, degree=3)",
8029 "y ~ te(x, z, k=5, penalty_order=2)",
8030 "y ~ te(x, z, k=5, degree=3, penalty_order=2)",
8031 ] {
8032 assert!(
8033 margins(formula).iter().all(is_cr),
8034 "{formula} must keep both cr margins"
8035 );
8036 }
8037
8038 for formula in [
8040 "y ~ te(x, z, k=5, degree=1)",
8041 "y ~ te(x, z, k=5, degree=4)",
8042 "y ~ te(x, z, k=5, penalty_order=1)",
8043 "y ~ te(x, z, k=5, penalty_order=3)",
8044 "y ~ te(x, z, k=5, knot_placement='uniform')",
8045 "y ~ te(x, z, k=5, knot_placement='quantile')",
8046 ] {
8047 assert!(
8048 margins(formula).iter().all(|k| !is_cr(k)),
8049 "{formula} must move both margins off the cr basis"
8050 );
8051 }
8052
8053 let per_margin = margins("y ~ te(x, z, k=5, degree=[1, 3])");
8055 assert!(!is_cr(&per_margin[0]), "the degree=1 margin must be a B-spline");
8056 assert!(is_cr(&per_margin[1]), "the degree=3 margin must stay cr");
8057
8058 let periodic = margins("y ~ te(x, z, k=5, periods=[1, None])");
8060 assert!(matches!(
8061 periodic[0],
8062 BSplineKnotSpec::PeriodicUniform { .. }
8063 ));
8064 assert!(is_cr(&periodic[1]));
8065 }
8066
8067 #[test]
8087 fn no_whitelisted_smooth_option_is_accepted_and_inert() {
8088 let structurally_inert = |kind: &str, key: &str| -> Option<&'static str> {
8091 match (kind, key) {
8092 (_, "type" | "bs") => Some("selects the arm; covered by the other rows"),
8096 (_, "by" | "__by_col") => Some("consumed by the BySmooth wrapper, not the arm"),
8100 (_, "ordered") => Some("qualifies a factor by=, read by the BySmooth wrapper"),
8103 (_, "id") => Some("shares a smoothing parameter; does not touch the basis"),
8107 ("cyclic", "double_penalty") => {
8116 Some("no null space survives the periodic sum-to-zero chart (#874)")
8117 }
8118 _ => None,
8119 }
8120 };
8121
8122 let probe = |kind: &str, key: &str| -> &'static str {
8125 match (kind, key) {
8126 ("tensor", "period" | "periods") => "[1.0, None]",
8129 ("tensor", "origin" | "origins" | "period_origin" | "period-origin"
8130 | "domain_origin") => "[0.0, None]",
8131 ("tensor", "periodic" | "cyclic") => "[0]",
8132 ("tensor", "boundary" | "bc") => "['periodic', 'natural']",
8133 (_, "periodic" | "cyclic") => "true",
8134 (_, "period" | "periods") => "0.7",
8135 (_, "period_start" | "start") => "0.05",
8136 (_, "period_end" | "end") => "0.7",
8137 (_, "origin" | "origins" | "period_origin" | "period-origin"
8138 | "domain_origin") => "0.1",
8139 (_, "boundary" | "bc" | "boundary_conditions") => "clamped",
8140 (_, "bc_left" | "left_bc" | "start_bc" | "bc_right" | "right_bc"
8141 | "end_bc") => "clamped",
8142 (_, "side") => "left",
8143 (_, "anchor" | "anchor_value" | "value" | "anchor_left"
8144 | "left_anchor" | "anchor_right" | "right_anchor") => "0.0",
8145 (_, "k" | "basis_dim" | "basis-dim" | "basisdim") => "6",
8147 (_, "centers") => "6",
8148 (_, "knots") => "13",
8149 (_, "knot_placement" | "knot-placement" | "knotplacement") => "quantile",
8150 (_, "degree") => "2",
8151 (_, "penalty_order" | "m") => "1",
8152 (_, "l" | "l_max" | "l-max" | "lmax" | "max_degree" | "max-degree") => "2",
8153 (_, "rank") => "5",
8154 (_, "order" | "nullspace_order") => "3",
8155 (_, "p" | "power") => "1.5",
8156 (_, "nu") => "1.5",
8157 (_, "kappa") => "0.5",
8158 (_, "alpha") => "0.5",
8159 (_, "tau") => "0.5",
8160 (_, "s" | "scales") => "3",
8161 (_, "length_scale") => "0.4",
8162 (_, "chunk_size") => "64",
8163 (_, "double_penalty") => "false",
8165 (_, "identifiability") => "none",
8166 (_, "include_intercept") => "true",
8167 (_, "scale_dims") => "true",
8168 (_, "multiscale") => "true",
8169 (_, "learn_length_scale") => "false",
8170 (_, "centered") => "false",
8171 (_, "smooth_penalty") => "false",
8172 (_, "lazy_path") => "true",
8173 (_, "radians") => "true",
8174 (_, "units") => "radians",
8175 (_, "kernel") => "pseudo",
8176 (_, "method") => "harmonic",
8177 (_, "path" | "pca_basis_path") => "'/nonexistent/pca.npy'",
8178 other => panic!(
8179 "no probe value for {other:?}; add one (or an exemption with a \
8180 reason) so the guard stays exhaustive"
8181 ),
8182 }
8183 };
8184
8185 let ds = continuous_dataset(
8190 &["y", "x", "z", "zbig", "lat", "lon"],
8191 (0..240)
8192 .map(|i| {
8193 let t = i as f64;
8194 let x = (i % 24) as f64 / 23.0;
8195 let z = (i / 24) as f64 / 9.0;
8196 vec![
8197 (t * 0.13).sin() + x + z,
8198 x,
8199 z,
8200 500.0 * z + 3.0,
8201 -80.0 + 160.0 * x,
8202 -170.0 + 340.0 * z,
8203 ]
8204 })
8205 .collect(),
8206 );
8207 let col_map = ds.column_map();
8208 let policy = gam_runtime::resource::ResourcePolicy::default_library();
8209 let build = |formula: &str| -> Result<String, String> {
8210 let parsed = parse_formula(formula)?;
8211 let mut notes = Vec::new();
8212 let spec = build_termspec(&parsed.terms, &ds, &col_map, &mut notes, &policy)
8213 .map_err(|err| err.to_string())?;
8214 let design = crate::smooth::build_term_collection_design(ds.values.view(), &spec)
8223 .map_err(|err| err.to_string())?;
8224 let dense = design.design.to_dense();
8225 let mut fingerprint = format!("design {}x{}", dense.nrows(), dense.ncols());
8226 for column in dense.columns() {
8227 let sum: f64 = column.iter().sum();
8228 let energy: f64 = column.iter().map(|v| v * v).sum();
8229 fingerprint.push_str(&format!(" |{sum:.10e},{energy:.10e}"));
8230 }
8231 for penalty in &design.smooth.penalties {
8232 let block = &penalty.local;
8233 let trace: f64 = (0..block.nrows()).map(|i| block[[i, i]]).sum();
8234 let energy: f64 = block.iter().map(|v| v * v).sum();
8235 fingerprint.push_str(&format!(
8236 " S[{}..{}]{}x{}:{trace:.10e},{energy:.10e}",
8237 penalty.col_range.start,
8238 penalty.col_range.end,
8239 block.nrows(),
8240 block.ncols()
8241 ));
8242 }
8243 Ok(fingerprint)
8244 };
8245
8246 let kinds: &[(&str, &str, &[&str])] = &[
8249 ("bspline", "s(x", BSPLINE_SMOOTH_OPTION_KEYS),
8250 ("cyclic", "cyclic(x", CYCLIC_SMOOTH_OPTION_KEYS),
8251 ("thinplate", "thinplate(x, zbig", THINPLATE_SMOOTH_OPTION_KEYS),
8252 ("matern", "matern(x, zbig", MATERN_SMOOTH_OPTION_KEYS),
8253 ("duchon", "duchon(x, zbig", DUCHON_SMOOTH_OPTION_KEYS),
8254 ("sphere", "sphere(lat, lon", SPHERE_SMOOTH_OPTION_KEYS),
8255 ("curvature", "curv(x, zbig", CURVATURE_SMOOTH_OPTION_KEYS),
8256 ("measurejet", "mjs(x, zbig", MEASURE_JET_SMOOTH_OPTION_KEYS),
8257 ("tensor", "te(x, z", TENSOR_SMOOTH_OPTION_KEYS),
8258 ];
8259
8260 let known_inert: &[(&str, &str)] = &[ (
8268 "y ~ thinplate(x, zbig, include_intercept=true)",
8269 "parsed into the spec, but the built radial design is unchanged",
8270 ),
8271 (
8272 "y ~ thinplate(x, zbig, scale_dims=true)",
8273 "parsed into the spec, but the built design is unchanged even on axes 500x apart in scale",
8274 ),
8275 (
8276 "y ~ matern(x, zbig, double_penalty=false)",
8277 "the flag does not change the shipped penalty set",
8278 ),
8279 (
8280 "y ~ curv(x, zbig, double_penalty=false)",
8281 "the flag does not change the shipped penalty set",
8282 ),
8283 ("y ~ mjs(x, zbig, tau=0.5)", "parsed, but the built design is unchanged"),
8284 (
8285 "y ~ mjs(x, zbig, learn_length_scale=false)",
8286 "parsed, but the built design is unchanged",
8287 ),
8288 ];
8289
8290 let mut inert = Vec::<String>::new();
8291 let mut honoured = 0usize;
8292 let mut refused = 0usize;
8293 for (kind, term, keys) in kinds {
8294 let baseline = match build(&format!("y ~ {term})")) {
8295 Ok(spec) => spec,
8296 Err(err) => panic!("baseline `y ~ {term})` must build, got: {err}"),
8297 };
8298 for key in *keys {
8299 if structurally_inert(kind, key).is_some() {
8300 continue;
8301 }
8302 let formula = format!("y ~ {term}, {key}={})", probe(kind, key));
8303 match build(&formula) {
8304 Err(_) => refused += 1,
8307 Ok(spec) if spec != baseline => honoured += 1,
8308 Ok(_) => inert.push(formula),
8309 }
8310 }
8311 }
8312
8313 let probed = honoured + refused + inert.len();
8317 assert!(
8318 probed >= 150,
8319 "the sweep should cover the whole option surface, only reached {probed} probes"
8320 );
8321 assert!(
8322 honoured * 2 > probed,
8323 "most probes should be HONOURED rather than refused, otherwise this \
8324 guard is testing error paths instead of option wiring \
8325 (honoured={honoured}, refused={refused}, inert={})",
8326 inert.len()
8327 );
8328
8329 let fixed: Vec<&str> = known_inert
8333 .iter()
8334 .map(|(formula, _)| *formula)
8335 .filter(|formula| !inert.iter().any(|found| found == formula))
8336 .collect();
8337 assert!(
8338 fixed.is_empty(),
8339 "these options are listed in `known_inert` but are no longer inert — \
8340 delete their entries so the list keeps telling the truth:\n {}",
8341 fixed.join("\n ")
8342 );
8343 inert.retain(|formula| {
8344 !known_inert
8345 .iter()
8346 .any(|(known, _)| known == formula)
8347 });
8348
8349 assert!(
8350 inert.is_empty(),
8351 "these formula options were accepted and produced a bit-identical \
8352 smooth design — each is either unwired (wire it), unsatisfiable in \
8353 this configuration (refuse it), or genuinely inert (exempt it in \
8354 `structurally_inert` with a reason). If it is a defect you are not \
8355 fixing right now, add it to `known_inert` WITH ITS REASON so the \
8356 ratchet still holds:\n {}",
8357 inert.join("\n ")
8358 );
8359 }
8360
8361 #[test]
8362 fn sz_factor_smooth_low_cardinality_uses_bspline_marginal() {
8363 let ds = ternary_factor_dataset();
8372 let col_map = ds.column_map();
8373 let parsed = parse_formula("y ~ s(x, g, bs=sz, k=10)").expect("parse sz factor smooth");
8374 let mut notes = Vec::new();
8375 let terms = build_termspec(
8376 &parsed.terms,
8377 &ds,
8378 &col_map,
8379 &mut notes,
8380 &gam_runtime::resource::ResourcePolicy::default_library(),
8381 )
8382 .expect("sz on a ternary covariate must build (B-spline marginal), not hard-fail");
8383 let SmoothBasisSpec::FactorSmooth { spec } = &terms.smooth_terms[0].basis else {
8384 panic!("expected FactorSmooth for s(x, g, bs=sz)");
8385 };
8386 assert!(
8387 !matches!(
8388 spec.marginal.knotspec,
8389 BSplineKnotSpec::NaturalCubicRegression { .. }
8390 ),
8391 "sz marginal must be a B-spline (curvature-capable), not the \
8392 natural-BC cr basis; got {:?}",
8393 spec.marginal.knotspec
8394 );
8395 }
8396
8397 fn continuous_x_factor_dataset(n: usize, n_groups: usize) -> Dataset {
8402 let rows = (0..n)
8403 .map(|i| {
8404 let x = i as f64 / (n as f64 - 1.0);
8405 let g = (i % n_groups) as f64;
8406 vec![x + g, x, g]
8407 })
8408 .collect::<Vec<_>>();
8409 let levels: Vec<String> = (0..n_groups).map(|k| format!("g{k}")).collect();
8410 Dataset {
8411 headers: vec!["y".into(), "x".into(), "g".into()],
8412 values: Array2::from_shape_vec(
8413 (rows.len(), 3),
8414 rows.into_iter().flat_map(|row| row.into_iter()).collect(),
8415 )
8416 .expect("rectangular continuous-x factor data"),
8417 schema: DataSchema {
8418 columns: vec![
8419 SchemaColumn {
8420 name: "y".into(),
8421 kind: ColumnKindTag::Continuous,
8422 levels: vec![],
8423 },
8424 SchemaColumn {
8425 name: "x".into(),
8426 kind: ColumnKindTag::Continuous,
8427 levels: vec![],
8428 },
8429 SchemaColumn {
8430 name: "g".into(),
8431 kind: ColumnKindTag::Categorical,
8432 levels,
8433 },
8434 ],
8435 },
8436 column_kinds: vec![
8437 ColumnKindTag::Continuous,
8438 ColumnKindTag::Continuous,
8439 ColumnKindTag::Categorical,
8440 ],
8441 }
8442 }
8443
8444 fn factor_smooth_spec_for(formula: &str, ds: &Dataset) -> FactorSmoothSpec {
8445 let col_map = ds.column_map();
8446 let parsed = parse_formula(formula).expect("parse factor smooth formula");
8447 let mut notes = Vec::new();
8448 let terms = build_termspec(
8449 &parsed.terms,
8450 ds,
8451 &col_map,
8452 &mut notes,
8453 &gam_runtime::resource::ResourcePolicy::default_library(),
8454 )
8455 .expect("build factor smooth term");
8456 let SmoothBasisSpec::FactorSmooth { spec } = &terms.smooth_terms[0].basis else {
8457 panic!("expected FactorSmooth basis for `{formula}`");
8458 };
8459 spec.clone()
8460 }
8461
8462 #[test]
8481 fn sz_factor_smooth_carries_null_space_ridge_like_fs() {
8482 let ds = continuous_x_factor_dataset(180, 4);
8483 let mut workspace = crate::basis::BasisWorkspace::new();
8484
8485 let sz_spec = factor_smooth_spec_for("y ~ s(x, g, bs=sz, k=8)", &ds);
8486 let sz_built = crate::smooth::build_factor_smooth(
8487 ds.values.view(),
8488 &sz_spec,
8489 "sz_term",
8490 &mut workspace,
8491 )
8492 .expect("build sz factor smooth");
8493
8494 let fs_spec = factor_smooth_spec_for("y ~ s(x, g, bs=fs, k=8)", &ds);
8495 let fs_built = crate::smooth::build_factor_smooth(
8496 ds.values.view(),
8497 &fs_spec,
8498 "fs_term",
8499 &mut workspace,
8500 )
8501 .expect("build fs factor smooth");
8502
8503 let n_levels = sz_spec
8522 .group_frozen_levels
8523 .as_ref()
8524 .map(|l| l.len())
8525 .unwrap_or(4);
8526 assert!(n_levels >= 3, "test needs >=3 groups, got {n_levels}");
8527
8528 let nw = 1usize; let expected_sz = fs_built.active_penalties.len() + (n_levels - 1) * nw;
8534 assert_eq!(
8535 sz_built.active_penalties.len(),
8536 expected_sz,
8537 "sz must split its wiggliness penalty per level (#1074): expected \
8538 fs_count {} + (L-1)·nw {} = {}, but sz had {}",
8539 fs_built.active_penalties.len(),
8540 (n_levels - 1) * nw,
8541 expected_sz,
8542 sz_built.active_penalties.len(),
8543 );
8544 assert!(
8545 sz_built.active_penalties.len() > fs_built.active_penalties.len(),
8546 "sz must carry strictly more penalties than fs after the per-group \
8547 split (sz={}, fs={})",
8548 sz_built.active_penalties.len(),
8549 fs_built.active_penalties.len(),
8550 );
8551
8552 let n_wiggliness = n_levels * nw; assert!(
8559 sz_built.active_penalties.len() > n_wiggliness,
8560 "sz deviation block carries no null-space ridge (penalties={}, \
8561 wiggliness blocks={}); the null space is unpenalized and REML \
8562 over-smooths the deviations",
8563 sz_built.active_penalties.len(),
8564 n_wiggliness,
8565 );
8566
8567 assert!(
8572 sz_built.dim < fs_built.dim,
8573 "sz design width {} must be strictly less than fs width {} \
8574 (zero-sum contrast drops one level block)",
8575 sz_built.dim,
8576 fs_built.dim,
8577 );
8578
8579 for penalty in &sz_built.active_penalties {
8580 assert_eq!(
8581 penalty
8582 .null_eigenvectors
8583 .as_ref()
8584 .map_or(0, |basis| basis.ncols()),
8585 penalty.nullity
8586 );
8587 }
8588 }
8589
8590 #[test]
8591 fn sz_penalty_metadata_is_emitted_in_matrix_order_2289() {
8592 let ds = continuous_x_factor_dataset(180, 4);
8593 let mut workspace = crate::basis::BasisWorkspace::new();
8594 let spec = factor_smooth_spec_for("y ~ s(x, g, bs=sz, k=8, double_penalty=true)", &ds);
8595 let built = crate::smooth::build_factor_smooth(
8596 ds.values.view(),
8597 &spec,
8598 "sz_metadata_order",
8599 &mut workspace,
8600 )
8601 .expect("build multi-penalty sz smooth");
8602 let n_levels = spec.group_frozen_levels.as_ref().map(Vec::len).unwrap_or(4);
8603
8604 assert!(built.active_penalties.len() >= 2 * n_levels);
8605 for (idx, penalty) in built.active_penalties.iter().enumerate() {
8606 let analysis =
8607 crate::basis::analyze_penalty_block(&penalty.matrix).expect("PSD penalty");
8608 assert_eq!(penalty.info.original_index, idx);
8609 assert_eq!(penalty.info.effective_rank, analysis.rank, "penalty {idx}");
8610 assert_eq!(penalty.nullity, analysis.nullity, "penalty {idx}");
8611 }
8612 assert!(
8613 built.active_penalties[..n_levels]
8614 .iter()
8615 .all(|penalty| matches!(penalty.info.source, PenaltySource::Primary))
8616 );
8617 assert!(
8618 built.active_penalties[n_levels..2 * n_levels]
8619 .iter()
8620 .all(|penalty| matches!(
8621 penalty.info.source,
8622 PenaltySource::DoublePenaltyNullspace
8623 ))
8624 );
8625 }
8626
8627 fn factor_dataset_l3() -> Dataset {
8638 let rows = (0..30)
8640 .map(|i| {
8641 let x = i as f64 / 29.0;
8642 let g = (i % 3) as f64;
8643 vec![x + g, x, g]
8644 })
8645 .collect::<Vec<_>>();
8646 Dataset {
8647 headers: vec!["y".into(), "x".into(), "g".into()],
8648 values: Array2::from_shape_vec(
8649 (rows.len(), 3),
8650 rows.into_iter().flat_map(|row| row.into_iter()).collect(),
8651 )
8652 .expect("rectangular L=3 factor test data"),
8653 schema: DataSchema {
8654 columns: vec![
8655 SchemaColumn {
8656 name: "y".into(),
8657 kind: ColumnKindTag::Continuous,
8658 levels: vec![],
8659 },
8660 SchemaColumn {
8661 name: "x".into(),
8662 kind: ColumnKindTag::Continuous,
8663 levels: vec![],
8664 },
8665 SchemaColumn {
8666 name: "g".into(),
8667 kind: ColumnKindTag::Categorical,
8668 levels: vec!["a".into(), "b".into(), "c".into()],
8669 },
8670 ],
8671 },
8672 column_kinds: vec![
8673 ColumnKindTag::Continuous,
8674 ColumnKindTag::Continuous,
8675 ColumnKindTag::Categorical,
8676 ],
8677 }
8678 }
8679
8680 #[test]
8681 fn factor_by_smooth_plus_bare_categorical_does_not_duplicate_factor_block() {
8682 let ds = factor_dataset_l3();
8683 let col_map = ds.column_map();
8684
8685 let g_blocks = |formula: &str| -> usize {
8686 let parsed = parse_formula(formula).expect("parse by-smooth formula");
8687 let mut notes = Vec::new();
8688 let terms = build_termspec(
8689 &parsed.terms,
8690 &ds,
8691 &col_map,
8692 &mut notes,
8693 &ResourcePolicy::default_library(),
8694 )
8695 .unwrap_or_else(|err| panic!("`{formula}` must build, got: {err:?}"));
8696 terms
8697 .random_effect_terms
8698 .iter()
8699 .filter(|rt| rt.name == "g")
8700 .count()
8701 };
8702
8703 let by_only = g_blocks("y ~ s(x, by=g, k=10)");
8707 assert_eq!(
8708 by_only, 1,
8709 "`y ~ s(x, by=g)` must produce exactly one `g` design block"
8710 );
8711
8712 let by_plus_bare = g_blocks("y ~ s(x, by=g, k=10) + g");
8716 assert_eq!(
8717 by_plus_bare, 1,
8718 "`y ~ s(x, by=g) + g` must collapse to ONE `g` block (#1457): the bare \
8719 `+ g` already owns the factor's level offsets, so the `by=` branch \
8720 must not add a second, treatment-coded main effect"
8721 );
8722
8723 assert_eq!(
8725 by_plus_bare, by_only,
8726 "the bare `+ g` collision must add zero extra `g` blocks (#1457)"
8727 );
8728 }
8729
8730 #[test]
8731 fn factor_by_penalties_carry_full_expanded_null_geometry_2293() {
8732 let ds = factor_dataset_l3();
8733 let col_map = ds.column_map();
8734 let parsed =
8739 parse_formula("y ~ s(x, by=g, k=8, double_penalty=false)").expect("parse by smooth");
8740 let mut notes = Vec::new();
8741 let terms = build_termspec(
8742 &parsed.terms,
8743 &ds,
8744 &col_map,
8745 &mut notes,
8746 &ResourcePolicy::default_library(),
8747 )
8748 .expect("build by smooth spec");
8749 assert_eq!(terms.smooth_terms.len(), 3, "one smooth per factor level");
8750
8751 for term in &terms.smooth_terms {
8756 assert!(matches!(
8757 &term.basis,
8758 SmoothBasisSpec::ByVariable {
8759 by: ByVariableSpec::Level { .. },
8760 ..
8761 }
8762 ));
8763 let mut workspace = crate::basis::BasisWorkspace::new();
8764 let built = crate::smooth::build_single_local_smooth_term(
8765 ds.values.view(),
8766 term,
8767 &mut workspace,
8768 )
8769 .expect("build level-gated factor-by smooth");
8770
8771 for (idx, penalty) in built.active_penalties.iter().enumerate() {
8772 let analysis =
8773 crate::basis::analyze_penalty_block(&penalty.matrix).expect("PSD block");
8774 assert_eq!(analysis.rank + penalty.nullity, built.dim, "penalty {idx}");
8775 assert_eq!(analysis.nullity, penalty.nullity, "penalty {idx}");
8776 assert_eq!(penalty.info.effective_rank, analysis.rank);
8777 let basis = penalty
8778 .null_eigenvectors
8779 .as_ref()
8780 .expect("nontrivial factor-level null basis");
8781 assert_eq!(basis.nrows(), built.dim);
8782 assert_eq!(basis.ncols(), penalty.nullity);
8783 }
8784 let joint = built
8785 .joint_null_rotation
8786 .as_ref()
8787 .expect("factor-level joint null geometry");
8788 assert!(joint.joint_nullity > 0);
8789 assert_eq!(joint.rotation.nrows(), built.dim);
8790 assert_eq!(joint.rotation.ncols(), built.dim);
8791 }
8792 }
8793
8794 #[test]
8795 fn parse_tensor_periods_and_origins_aliases() {
8796 let mut opts = BTreeMap::new();
8797 opts.insert(
8798 "boundary".to_string(),
8799 "['periodic', 'periodic']".to_string(),
8800 );
8801 opts.insert("periods".to_string(), "[7, 24]".to_string());
8802 opts.insert("origins".to_string(), "[0, -12]".to_string());
8803 let axes = parse_periodic_axes(&opts, 2).expect("axes");
8804 let periods = parse_periods(&opts, &axes).expect("periods");
8805 let origins = parse_period_origins(&opts, &axes).expect("origins");
8806 assert_eq!(axes, vec![true, true]);
8807 assert_eq!(periods, vec![Some(7.0), Some(24.0)]);
8808 assert_eq!(origins, vec![Some(0.0), Some(-12.0)]);
8809 }
8810
8811 #[test]
8812 fn tensor_smooth_honors_per_margin_k_list() {
8813 let ds = continuous_dataset(
8814 &["y", "theta", "h"],
8815 (0..20)
8816 .map(|i| {
8817 let theta = std::f64::consts::TAU * i as f64 / 20.0;
8818 let h = -1.0 + 2.0 * (i % 5) as f64 / 4.0;
8819 vec![theta.cos() + h, theta, h]
8820 })
8821 .collect(),
8822 );
8823 let parsed = parse_formula(
8824 "y ~ te(theta, h, periodic=[0], period=[2*pi, None], origin=[0, None], k=[9,5])",
8825 )
8826 .expect("parse tensor formula");
8827 let col_map = ds.column_map();
8828 let mut notes = Vec::new();
8829 let terms = build_termspec(
8830 &parsed.terms,
8831 &ds,
8832 &col_map,
8833 &mut notes,
8834 &gam_runtime::resource::ResourcePolicy::default_library(),
8835 )
8836 .expect("build tensor terms");
8837 let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
8838 panic!("expected tensor B-spline");
8839 };
8840 let dims = spec
8841 .marginalspecs
8842 .iter()
8843 .map(|m| match m.knotspec {
8844 BSplineKnotSpec::PeriodicUniform { num_basis, .. } => num_basis,
8845 BSplineKnotSpec::Generate {
8846 num_internal_knots, ..
8847 } => num_internal_knots + m.degree + 1,
8848 BSplineKnotSpec::NaturalCubicRegression { ref knots } => knots.len(),
8851 _ => panic!("unexpected tensor marginal knotspec"),
8852 })
8853 .collect::<Vec<_>>();
8854 assert_eq!(dims, vec![9, 5]);
8855 }
8856
8857 #[test]
8858 fn tensor_smooth_honors_per_margin_k_axis_aliases() {
8859 let ds = continuous_dataset(
8860 &["resp", "x", "y"],
8861 (0..12)
8862 .map(|i| {
8863 let t = i as f64 / 11.0;
8864 vec![t, t, 1.0 - t]
8865 })
8866 .collect(),
8867 );
8868 assert_eq!(
8869 tensor_margin_basis_sizes(&ds, "resp ~ te(x, y, k_x=9, k_y=5)"),
8870 vec![9, 5],
8871 "k_<margin> aliases should materialize requested per-margin values"
8872 );
8873 }
8874
8875 #[test]
8876 fn tensor_smooth_low_cardinality_axis_falls_back_to_lower_degree_basis() {
8877 let ds = continuous_dataset(
8884 &["y", "x", "b"],
8885 (0..40)
8886 .map(|i| {
8887 let x = i as f64 / 39.0;
8888 let b = (i % 2) as f64;
8889 vec![x.sin() + 0.5 * b, x, b]
8890 })
8891 .collect(),
8892 );
8893 let parsed = parse_formula("y ~ te(x, b, k=[5, 2])").expect("parse tensor with k=[5,2]");
8894 let col_map = ds.column_map();
8895 let mut notes = Vec::new();
8896 let terms = build_termspec(
8897 &parsed.terms,
8898 &ds,
8899 &col_map,
8900 &mut notes,
8901 &gam_runtime::resource::ResourcePolicy::default_library(),
8902 )
8903 .expect("build tensor with binary margin");
8904 let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
8905 panic!("expected tensor B-spline for te(x, b)");
8906 };
8907 let continuous = &spec.marginalspecs[0];
8911 let binary = &spec.marginalspecs[1];
8912 assert_eq!(continuous.degree, 3);
8913 assert_eq!(binary.degree, 1);
8914 assert!(
8915 binary.penalty_order >= 1 && binary.penalty_order <= binary.degree,
8916 "binary margin penalty_order {} must satisfy 1 <= order <= degree={}",
8917 binary.penalty_order,
8918 binary.degree
8919 );
8920 let basis_size = |m: &BSplineBasisSpec| match m.knotspec {
8921 BSplineKnotSpec::PeriodicUniform { num_basis, .. } => num_basis,
8922 BSplineKnotSpec::Generate {
8923 num_internal_knots, ..
8924 } => num_internal_knots + m.degree + 1,
8925 BSplineKnotSpec::Automatic {
8926 num_internal_knots: Some(n),
8927 ..
8928 } => n + m.degree + 1,
8929 BSplineKnotSpec::NaturalCubicRegression { ref knots } => knots.len(),
8932 _ => panic!("unexpected tensor marginal knotspec"),
8933 };
8934 assert_eq!(basis_size(continuous), 5);
8935 assert_eq!(basis_size(binary), 2);
8936 }
8937
8938 #[test]
8939 fn tensor_smooth_uniform_k_is_capped_to_a_low_cardinality_margins_distinct_values() {
8940 let ds = continuous_dataset(
8948 &["y", "x", "b"],
8949 (0..40)
8950 .map(|i| {
8951 let x = i as f64 / 39.0;
8952 let b = (i % 2) as f64;
8953 vec![x.sin() + 0.5 * b, x, b]
8954 })
8955 .collect(),
8956 );
8957 let parsed = parse_formula("y ~ te(x, b, k=5)").expect("parse tensor with uniform k=5");
8958 let col_map = ds.column_map();
8959 let mut notes = Vec::new();
8960 let terms = build_termspec(
8961 &parsed.terms,
8962 &ds,
8963 &col_map,
8964 &mut notes,
8965 &gam_runtime::resource::ResourcePolicy::default_library(),
8966 )
8967 .expect("uniform k=5 must auto-cap the binary margin instead of erroring");
8968 let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
8969 panic!("expected tensor B-spline for te(x, b)");
8970 };
8971 let basis_size = |m: &BSplineBasisSpec| match &m.knotspec {
8972 BSplineKnotSpec::PeriodicUniform { num_basis, .. } => *num_basis,
8973 BSplineKnotSpec::Generate {
8974 num_internal_knots, ..
8975 } => num_internal_knots + m.degree + 1,
8976 BSplineKnotSpec::Automatic {
8977 num_internal_knots: Some(n),
8978 ..
8979 } => n + m.degree + 1,
8980 BSplineKnotSpec::NaturalCubicRegression { knots } => knots.len(),
8981 other => panic!("unexpected tensor marginal knotspec: {other:?}"),
8982 };
8983 let binary = &spec.marginalspecs[1];
8984 assert_eq!(basis_size(binary), 2);
8987 assert_eq!(binary.degree, 1);
8988 assert_eq!(basis_size(&spec.marginalspecs[0]), 5);
8990 }
8991
8992 #[test]
8993 fn tensor_all_tp_margins_with_per_margin_k_routes_to_bspline_tensor() {
8994 let ds = continuous_dataset(
9003 &["y", "x1", "x2"],
9004 (0..32)
9005 .map(|i| {
9006 let t = i as f64 / 31.0;
9007 vec![t.sin(), t, 1.0 - t]
9008 })
9009 .collect(),
9010 );
9011 let parsed =
9012 parse_formula("y ~ te(x1, x2, bs=c('tp','tp'), k=c(5,5))").expect("parse tensor");
9013 let col_map = ds.column_map();
9014 let mut notes = Vec::new();
9015 let terms = build_termspec(
9016 &parsed.terms,
9017 &ds,
9018 &col_map,
9019 &mut notes,
9020 &gam_runtime::resource::ResourcePolicy::default_library(),
9021 )
9022 .expect("build tensor terms with per-margin k");
9023 let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
9024 panic!(
9025 "expected B-spline tensor when k=c(5,5) is supplied with bs=c('tp','tp'), got {:?}",
9026 terms.smooth_terms[0].basis
9027 );
9028 };
9029 let dims = spec
9039 .marginalspecs
9040 .iter()
9041 .map(|m| match m.knotspec {
9042 BSplineKnotSpec::Generate {
9043 num_internal_knots, ..
9044 } => num_internal_knots + m.degree + 1,
9045 BSplineKnotSpec::Automatic {
9046 num_internal_knots: Some(num_internal_knots),
9047 ..
9048 } => num_internal_knots + m.degree + 1,
9049 BSplineKnotSpec::PeriodicUniform { num_basis, .. } => num_basis,
9050 BSplineKnotSpec::Provided(ref knots) => knots.len().saturating_sub(m.degree + 1),
9051 BSplineKnotSpec::NaturalCubicRegression { ref knots } => knots.len(),
9052 BSplineKnotSpec::Automatic {
9053 num_internal_knots: None,
9054 ..
9055 } => panic!("test cannot infer automatic knot count"),
9056 })
9057 .collect::<Vec<_>>();
9058 assert_eq!(dims, vec![5, 5]);
9059 }
9060
9061 #[test]
9062 fn tensor_all_tp_margins_without_per_margin_k_builds_anisotropic_tensor() {
9063 let ds = continuous_dataset(
9071 &["y", "x1", "x2"],
9072 (0..32)
9073 .map(|i| {
9074 let t = i as f64 / 31.0;
9075 vec![t.sin(), t, 1.0 - t]
9076 })
9077 .collect(),
9078 );
9079 let parsed = parse_formula("y ~ te(x1, x2, bs=c('tp','tp'))").expect("parse tensor");
9080 let col_map = ds.column_map();
9081 let mut notes = Vec::new();
9082 let terms = build_termspec(
9083 &parsed.terms,
9084 &ds,
9085 &col_map,
9086 &mut notes,
9087 &gam_runtime::resource::ResourcePolicy::default_library(),
9088 )
9089 .expect("build tensor terms without per-margin k");
9090 let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
9091 panic!(
9092 "te(...,bs=c('tp','tp')) must route to an anisotropic tensor product, not a \
9093 silent isotropic thin-plate substitution; got {:?}",
9094 terms.smooth_terms[0].basis
9095 );
9096 };
9097 assert_eq!(
9098 spec.marginalspecs.len(),
9099 2,
9100 "tp tensor must carry one penalized B-spline margin per axis"
9101 );
9102 }
9103
9104 #[test]
9105 fn explicit_basis_sizes_are_not_small_n_clamped() {
9106 let ds = continuous_dataset(
9107 &["y", "x1", "x2", "x3", "x4", "x5"],
9108 (0..12)
9109 .map(|i| {
9110 let x = i as f64 / 11.0;
9111 vec![x.sin(), x, x * x, x + 0.1, 1.0 - x, (2.0 * x).sin()]
9112 })
9113 .collect(),
9114 );
9115 let parsed = parse_formula("y ~ s(x1, k=10) + s(x2) + s(x3) + s(x4) + s(x5)")
9116 .expect("parse multi-smooth formula");
9117 let col_map = ds.column_map();
9118 let mut notes = Vec::new();
9119 let terms = build_termspec(
9120 &parsed.terms,
9121 &ds,
9122 &col_map,
9123 &mut notes,
9124 &gam_runtime::resource::ResourcePolicy::default_library(),
9125 )
9126 .expect("build multi-smooth terms");
9127 let SmoothBasisSpec::BSpline1D { spec, .. } = &terms.smooth_terms[0].basis else {
9128 panic!("expected first smooth to be B-spline");
9129 };
9130 assert!(matches!(
9131 &spec.knotspec,
9132 BSplineKnotSpec::Generate {
9133 num_internal_knots: 6,
9134 ..
9135 }
9136 ));
9137 }
9138
9139 #[test]
9140 fn explicit_duchon_centers_are_not_small_n_bumped() {
9141 let ds = continuous_dataset(
9142 &["y", "x1", "x2", "x3", "x4", "x5"],
9143 (0..12)
9144 .map(|i| {
9145 let x = i as f64 / 11.0;
9146 vec![x.sin(), x, x * x, x + 0.1, 1.0 - x, (2.0 * x).sin()]
9147 })
9148 .collect(),
9149 );
9150 let parsed = parse_formula("y ~ duchon(x1, centers=3) + s(x2) + s(x3) + s(x4) + s(x5)")
9157 .expect("parse multi-smooth formula");
9158 let col_map = ds.column_map();
9159 let mut notes = Vec::new();
9160 let terms = build_termspec(
9161 &parsed.terms,
9162 &ds,
9163 &col_map,
9164 &mut notes,
9165 &gam_runtime::resource::ResourcePolicy::default_library(),
9166 )
9167 .expect("build multi-smooth terms");
9168 let SmoothBasisSpec::Duchon { spec, .. } = &terms.smooth_terms[0].basis else {
9169 panic!("expected first smooth to be Duchon");
9170 };
9171 assert!(matches!(
9172 spec.center_strategy,
9173 CenterStrategy::UniformGrid { points_per_dim: 3 }
9174 ));
9175 }
9176
9177 #[test]
9178 fn inferred_tensor_basis_cap_uses_coordinate_support_not_duplicate_rows() {
9179 let mut unique_rows = Vec::new();
9180 for i in 0..50 {
9181 let theta = i as f64 / 50.0;
9182 for j in 0..16 {
9183 let h = -1.0 + 2.0 * (j as f64) / 15.0;
9184 let y = theta.cos() + h;
9185 unique_rows.push(vec![y, theta, h]);
9186 }
9187 }
9188 let mut repeated_rows = Vec::new();
9189 for _ in 0..12 {
9190 repeated_rows.extend(unique_rows.iter().cloned());
9191 }
9192
9193 let unique = continuous_dataset(&["y", "theta", "h"], unique_rows);
9194 let repeated = continuous_dataset(&["y", "theta", "h"], repeated_rows);
9195
9196 let unique_basis = inferred_tensor_basis_product(&unique);
9197 let repeated_basis = inferred_tensor_basis_product(&repeated);
9198
9199 assert_eq!(
9200 unique_basis, repeated_basis,
9201 "duplicating existing tensor coordinates must not inflate inferred basis width"
9202 );
9203 }
9204
9205 #[test]
9206 fn inferred_three_dim_tensor_basis_stays_bounded_for_reml_selection() {
9207 let make = |n: usize| -> usize {
9215 let mut rows = Vec::with_capacity(n);
9216 for i in 0..n {
9217 let f = i as f64 / n as f64;
9218 rows.push(vec![f.sin(), f, (2.0 * f).cos(), (3.0 * f) % 1.0]);
9219 }
9220 let ds = continuous_dataset(&["y", "x1", "x2", "x3"], rows);
9221 let parsed = parse_formula("y ~ te(x1, x2, x3)").expect("parse 3-D tensor");
9222 let col_map = ds.column_map();
9223 let mut notes = Vec::new();
9224 let terms = build_termspec(
9225 &parsed.terms,
9226 &ds,
9227 &col_map,
9228 &mut notes,
9229 &ResourcePolicy::default_library(),
9230 )
9231 .expect("build 3-D tensor termspec");
9232 let SmoothBasisSpec::TensorBSpline { spec, .. } = &terms.smooth_terms[0].basis else {
9233 panic!("expected tensor smooth");
9234 };
9235 spec.marginalspecs
9236 .iter()
9237 .map(|m| match m.knotspec {
9238 BSplineKnotSpec::Generate {
9239 num_internal_knots, ..
9240 } => num_internal_knots + m.degree + 1,
9241 BSplineKnotSpec::Automatic {
9242 num_internal_knots: Some(num_internal_knots),
9243 ..
9244 } => num_internal_knots + m.degree + 1,
9245 BSplineKnotSpec::NaturalCubicRegression { ref knots } => knots.len(),
9248 _ => panic!("unexpected tensor margin knotspec"),
9249 })
9250 .product()
9251 };
9252
9253 assert!(
9255 make(60) <= 216,
9256 "3-D te at small n must stay near the mgcv te default, got {}",
9257 make(60)
9258 );
9259 assert!(
9261 make(2000) <= 216,
9262 "3-D te at large n must not blow ∏k toward the data size, got {}",
9263 make(2000)
9264 );
9265 }
9266
9267 #[test]
9268 fn parse_bspline_boundary_conditions_and_side_selector() {
9269 let mut opts = BTreeMap::new();
9273 opts.insert("boundary_conditions".to_string(), "anchored".to_string());
9274 opts.insert("side".to_string(), "left".to_string());
9275 opts.insert("anchor".to_string(), "2.5".to_string());
9276 let parsed = parse_bspline_boundary_conditions(&opts).expect("left anchor parses");
9277 assert!(matches!(
9278 parsed.left,
9279 BSplineEndpointBoundaryCondition::Anchored { value } if value == 2.5
9280 ));
9281 assert!(matches!(
9282 parsed.right,
9283 BSplineEndpointBoundaryCondition::Free
9284 ));
9285
9286 let mut opts = BTreeMap::new();
9290 opts.insert("start_bc".to_string(), "clamped".to_string());
9291 opts.insert("end_bc".to_string(), "zero".to_string());
9292 opts.insert("right_anchor".to_string(), "-1.0".to_string());
9293 let parsed = parse_bspline_boundary_conditions(&opts).expect("right anchor parses");
9294 assert!(matches!(
9295 parsed.left,
9296 BSplineEndpointBoundaryCondition::Clamped
9297 ));
9298 assert!(matches!(
9299 parsed.right,
9300 BSplineEndpointBoundaryCondition::Anchored { value } if value == -1.0
9301 ));
9302
9303 let mut opts = BTreeMap::new();
9307 opts.insert("start_bc".to_string(), "clamped".to_string());
9308 opts.insert("end_bc".to_string(), "zero".to_string());
9309 let parsed = parse_bspline_boundary_conditions(&opts).expect("boundary conditions");
9310 assert!(matches!(
9311 parsed.left,
9312 BSplineEndpointBoundaryCondition::Clamped
9313 ));
9314 assert!(matches!(
9315 parsed.right,
9316 BSplineEndpointBoundaryCondition::Anchored { value } if value.abs() < 1e-12
9317 ));
9318 }
9319
9320 #[test]
9321 fn one_sided_anchor_owns_level_without_sum_to_zero_constraint_1867() {
9322 let ds = continuous_dataset(
9323 &["y", "x"],
9324 (0..32)
9325 .map(|i| {
9326 let x = i as f64 / 31.0;
9327 vec![x * (1.0 - x), x]
9328 })
9329 .collect(),
9330 );
9331 let col_map = ds.column_map();
9332
9333 let build = |formula: &str| {
9334 let parsed = parse_formula(formula).expect("parse anchored smooth");
9335 let mut notes = Vec::new();
9336 build_termspec(
9337 &parsed.terms,
9338 &ds,
9339 &col_map,
9340 &mut notes,
9341 &ResourcePolicy::default_library(),
9342 )
9343 .expect("build anchored smooth")
9344 };
9345
9346 let one_sided = build("y ~ s(x, bc_left=anchored, anchor_left=0, k=10)");
9347 let SmoothBasisSpec::BSpline1D { spec, .. } = &one_sided.smooth_terms[0].basis else {
9348 panic!("expected one-dimensional B-spline");
9349 };
9350 assert!(matches!(spec.identifiability, BSplineIdentifiability::None));
9351
9352 let two_sided = build("y ~ s(x, bc_left=anchored, bc_right=anchored, k=10)");
9359 let SmoothBasisSpec::BSpline1D { spec, .. } = &two_sided.smooth_terms[0].basis else {
9360 panic!("expected one-dimensional B-spline");
9361 };
9362 assert!(matches!(spec.identifiability, BSplineIdentifiability::None));
9363
9364 let plain = build("y ~ s(x, k=10)");
9367 let SmoothBasisSpec::BSpline1D { spec, .. } = &plain.smooth_terms[0].basis else {
9368 panic!("expected one-dimensional B-spline");
9369 };
9370 assert!(matches!(
9371 spec.identifiability,
9372 BSplineIdentifiability::WeightedSumToZero { .. }
9373 ));
9374 }
9375
9376 #[test]
9377 fn categorical_by_numeric_interaction_expands_treatment_coded_cells() {
9378 let ds = factor_dataset();
9389 let parsed = parse_formula("y ~ x:g").expect("parse `y ~ x:g`");
9391 let col_map = ds.column_map();
9392 let mut notes = Vec::new();
9393 let terms = build_termspec(
9394 &parsed.terms,
9395 &ds,
9396 &col_map,
9397 &mut notes,
9398 &ResourcePolicy::default_library(),
9399 )
9400 .expect("factor-aware `x:g` interaction must build, not error");
9401
9402 assert_eq!(
9403 terms.linear_terms.len(),
9404 2,
9405 "interaction-only `x:g` keeps ALL factor levels (full dummy coding): one slope column per group"
9406 );
9407
9408 let x_col = *col_map.get("x").expect("x column");
9409 let g_col = *col_map.get("g").expect("g column");
9410
9411 let mut seen_bits = std::collections::HashSet::new();
9414 for term in &terms.linear_terms {
9415 assert!(
9416 term.is_interaction(),
9417 "the categorical-by-numeric cell is a Wilkinson-Rogers interaction"
9418 );
9419 assert_eq!(term.feature_cols, vec![x_col]);
9420 assert_eq!(term.categorical_levels.len(), 1);
9421 let (gate_col, gate_bits) = term.categorical_levels[0];
9422 assert_eq!(gate_col, g_col);
9423 assert!(seen_bits.insert(gate_bits), "each level appears once");
9424
9425 let column = term
9427 .realized_design_column(ds.values.view())
9428 .expect("realize cell column");
9429 let n = ds.values.nrows();
9430 assert_eq!(column.len(), n);
9431 for row in 0..n {
9432 let x = ds.values[[row, x_col]];
9433 let g = ds.values[[row, g_col]];
9434 let expected = if g.to_bits() == gate_bits { x } else { 0.0 };
9435 assert!(
9436 (column[row] - expected).abs() < 1e-12,
9437 "row {row}: g={g}, x={x}, expected {expected}, got {}",
9438 column[row]
9439 );
9440 }
9441 }
9442 assert!(seen_bits.contains(&0.0_f64.to_bits()));
9445 assert!(seen_bits.contains(&1.0_f64.to_bits()));
9446 }
9447
9448 #[test]
9449 fn categorical_by_numeric_interaction_keeps_treatment_coding_with_parent() {
9450 let ds = factor_dataset();
9458 let parsed = parse_formula("y ~ x + x:g").expect("parse `y ~ x + x:g`");
9459 let col_map = ds.column_map();
9460 let mut notes = Vec::new();
9461 let terms = build_termspec(
9462 &parsed.terms,
9463 &ds,
9464 &col_map,
9465 &mut notes,
9466 &ResourcePolicy::default_library(),
9467 )
9468 .expect("`x + x:g` must build");
9469
9470 let x_col = *col_map.get("x").expect("x column");
9472 let g_col = *col_map.get("g").expect("g column");
9473 let interaction_cells: Vec<_> = terms
9474 .linear_terms
9475 .iter()
9476 .filter(|t| t.is_interaction())
9477 .collect();
9478 assert_eq!(
9479 interaction_cells.len(),
9480 1,
9481 "with `x` present, `x:g` is treatment-coded → one cell (reference dropped)"
9482 );
9483 let term = interaction_cells[0];
9484 assert_eq!(term.feature_cols, vec![x_col]);
9485 assert_eq!(term.categorical_levels.len(), 1);
9486 let (gate_col, gate_bits) = term.categorical_levels[0];
9487 assert_eq!(gate_col, g_col);
9488 assert_eq!(gate_bits, 1.0_f64.to_bits());
9490 }
9491
9492 #[test]
9493 fn categorical_by_categorical_interaction_expands_full_cross_cells() {
9494 let n = 30usize;
9505 let mut rows = Vec::with_capacity(n);
9506 for i in 0..n {
9507 let y = (i as f64).sin();
9508 let f = (i % 3) as f64; let g = (i % 2) as f64; rows.push(vec![y, f, g]);
9511 }
9512 let values = Array2::from_shape_vec(
9513 (n, 3),
9514 rows.into_iter().flat_map(|row| row.into_iter()).collect(),
9515 )
9516 .expect("rectangular cross-factor data");
9517 let ds = Dataset {
9518 headers: vec!["y".into(), "f".into(), "g".into()],
9519 values,
9520 schema: DataSchema {
9521 columns: vec![
9522 SchemaColumn {
9523 name: "y".into(),
9524 kind: ColumnKindTag::Continuous,
9525 levels: vec![],
9526 },
9527 SchemaColumn {
9528 name: "f".into(),
9529 kind: ColumnKindTag::Categorical,
9530 levels: vec!["f0".into(), "f1".into(), "f2".into()],
9531 },
9532 SchemaColumn {
9533 name: "g".into(),
9534 kind: ColumnKindTag::Categorical,
9535 levels: vec!["g0".into(), "g1".into()],
9536 },
9537 ],
9538 },
9539 column_kinds: vec![
9540 ColumnKindTag::Continuous,
9541 ColumnKindTag::Categorical,
9542 ColumnKindTag::Categorical,
9543 ],
9544 };
9545
9546 let parsed = parse_formula("y ~ f:g").expect("parse `y ~ f:g`");
9547 let col_map = ds.column_map();
9548 let mut notes = Vec::new();
9549 let terms = build_termspec(
9550 &parsed.terms,
9551 &ds,
9552 &col_map,
9553 &mut notes,
9554 &ResourcePolicy::default_library(),
9555 )
9556 .expect("factor-by-factor `f:g` interaction must build, not error");
9557
9558 assert_eq!(
9559 terms.linear_terms.len(),
9560 5,
9561 "saturated 3*2 = 6 cross cells minus one reference cell (f0:g0) = 5"
9562 );
9563
9564 let f_col = *col_map.get("f").expect("f column");
9565 let g_col = *col_map.get("g").expect("g column");
9566 let f0 = 0.0_f64.to_bits();
9570 let g0 = 0.0_f64.to_bits();
9571 let mut emitted = std::collections::HashSet::new();
9572 for term in &terms.linear_terms {
9573 assert!(term.feature_cols.is_empty());
9575 assert_eq!(term.categorical_levels.len(), 2);
9576 let mut gates = std::collections::HashMap::new();
9577 for &(col, bits) in &term.categorical_levels {
9578 gates.insert(col, bits);
9579 }
9580 let f_bits = *gates.get(&f_col).expect("f gate present");
9581 let g_bits = *gates.get(&g_col).expect("g gate present");
9582 assert!(
9584 !(f_bits == f0 && g_bits == g0),
9585 "the reference cell f0:g0 must be absorbed by the intercept, not emitted"
9586 );
9587 emitted.insert((f_bits, g_bits));
9588
9589 let column = term
9590 .realized_design_column(ds.values.view())
9591 .expect("realize cross cell");
9592 for row in 0..n {
9593 let f = ds.values[[row, f_col]];
9594 let g = ds.values[[row, g_col]];
9595 let expected = if f.to_bits() == f_bits && g.to_bits() == g_bits {
9596 1.0
9597 } else {
9598 0.0
9599 };
9600 assert!(
9601 (column[row] - expected).abs() < 1e-12,
9602 "row {row}: expected {expected}, got {}",
9603 column[row]
9604 );
9605 }
9606 assert!(
9607 column.iter().any(|&v| v == 1.0),
9608 "each cross cell must be observed in the data"
9609 );
9610 }
9611 let f_levels = [0.0_f64.to_bits(), 1.0_f64.to_bits(), 2.0_f64.to_bits()];
9614 let g_levels = [0.0_f64.to_bits(), 1.0_f64.to_bits()];
9615 for &fb in &f_levels {
9616 for &gb in &g_levels {
9617 if fb == f0 && gb == g0 {
9618 continue;
9619 }
9620 assert!(
9621 emitted.contains(&(fb, gb)),
9622 "saturated cross cell must be present"
9623 );
9624 }
9625 }
9626 }
9627
9628 #[test]
9634 fn by_level_thin_plate_sizes_default_centers_from_the_smallest_level() {
9635 let n_a = 60usize;
9636 let n_b = 180usize;
9637 let rows: Vec<Vec<f64>> = (0..(n_a + n_b))
9638 .map(|i| {
9639 let in_a = i < n_a;
9640 let x = if in_a {
9641 i as f64 / (n_a - 1) as f64
9642 } else {
9643 (i - n_a) as f64 / (n_b - 1) as f64
9644 };
9645 let g = if in_a { 0.0 } else { 1.0 };
9646 vec![x + g, x, g]
9647 })
9648 .collect();
9649 let ds = Dataset {
9650 headers: vec!["y".into(), "x".into(), "g".into()],
9651 values: Array2::from_shape_vec(
9652 (rows.len(), 3),
9653 rows.into_iter().flat_map(|row| row.into_iter()).collect(),
9654 )
9655 .expect("rectangular by-level test data"),
9656 schema: DataSchema {
9657 columns: vec![
9658 SchemaColumn {
9659 name: "y".into(),
9660 kind: ColumnKindTag::Continuous,
9661 levels: vec![],
9662 },
9663 SchemaColumn {
9664 name: "x".into(),
9665 kind: ColumnKindTag::Continuous,
9666 levels: vec![],
9667 },
9668 SchemaColumn {
9669 name: "g".into(),
9670 kind: ColumnKindTag::Categorical,
9671 levels: vec!["a".into(), "b".into()],
9672 },
9673 ],
9674 },
9675 column_kinds: vec![
9676 ColumnKindTag::Continuous,
9677 ColumnKindTag::Continuous,
9678 ColumnKindTag::Categorical,
9679 ],
9680 };
9681 let build_tp = |with_by: bool| -> SmoothBasisSpec {
9682 let mut options = BTreeMap::new();
9683 options.insert("bs".to_string(), "tps".to_string());
9684 if with_by {
9685 options.insert("by".to_string(), "g".to_string());
9686 options.insert("__by_col".to_string(), "2".to_string());
9687 }
9688 let mut notes = Vec::new();
9689 build_smooth_basis(
9690 SmoothKind::S,
9691 &["x".to_string()],
9692 &[1],
9693 &options,
9694 &ds,
9695 &mut notes,
9696 &ResourcePolicy::default_library(),
9697 1,
9698 )
9699 .expect("thin-plate basis builds")
9700 };
9701 let pooled = build_tp(false);
9702 let by_level = build_tp(true);
9703 let tp_centers = |basis: &SmoothBasisSpec| -> usize {
9704 match basis {
9705 SmoothBasisSpec::ThinPlate { spec, .. } => {
9706 spec.center_strategy.planned_num_centers(1)
9707 }
9708 SmoothBasisSpec::BySmooth { smooth, .. } => match smooth.as_ref() {
9709 SmoothBasisSpec::ThinPlate { spec, .. } => {
9710 spec.center_strategy.planned_num_centers(1)
9711 }
9712 other => panic!("expected ThinPlate inside BySmooth, got {other:?}"),
9713 },
9714 other => panic!("expected ThinPlate, got {other:?}"),
9715 }
9716 };
9717 let pooled_centers = tp_centers(&pooled);
9718 let by_centers = tp_centers(&by_level);
9719 assert!(
9720 by_centers < pooled_centers,
9721 "by-level default centers must size from the smallest level: \
9722 by={by_centers} pooled={pooled_centers}"
9723 );
9724 let ds_small = continuous_dataset(
9727 &["y", "x"],
9728 (0..n_a)
9729 .map(|i| {
9730 let x = i as f64 / (n_a - 1) as f64;
9731 vec![x, x]
9732 })
9733 .collect(),
9734 );
9735 let mut small_options = BTreeMap::new();
9736 small_options.insert("bs".to_string(), "tps".to_string());
9737 let mut notes = Vec::new();
9738 let small = build_smooth_basis(
9739 SmoothKind::S,
9740 &["x".to_string()],
9741 &[1],
9742 &small_options,
9743 &ds_small,
9744 &mut notes,
9745 &ResourcePolicy::default_library(),
9746 1,
9747 )
9748 .expect("small-level thin-plate basis builds");
9749 assert_eq!(
9750 by_centers,
9751 tp_centers(&small),
9752 "by-level default must equal the smallest level's own default"
9753 );
9754 }
9755}