1use coefficient_transforms::{
2 convex_divided_difference_transform_matrix, cumulative_exp, cumulative_sum_transform_matrix,
3 second_cumulative_exp,
4};
5
6pub use error::SmoothError;
7
8use input_standardization::{
9 apply_input_standardization, compensate_length_scale_for_standardization,
10 compensate_optional_length_scale_for_standardization, compute_spatial_input_scales,
11};
12
13use shape_constraints::{
14 build_shape_constraint_design_1d, build_shape_linear_constraints_1d,
15 merge_linear_constraints_global, shape_lower_bounds_local, shape_order_and_sign,
16 shape_supports_basis, shape_uses_box_reparameterization,
17};
18
19pub fn describe_thin_plate_center_request(strategy: &CenterStrategy) -> String {
20 match strategy {
21 CenterStrategy::Auto(inner) => describe_thin_plate_center_request(inner),
22 CenterStrategy::UserProvided(centers) => format!("{} centers", centers.nrows()),
23 CenterStrategy::EqualMass { num_centers }
24 | CenterStrategy::EqualMassCovarRepresentative { num_centers }
25 | CenterStrategy::FarthestPoint { num_centers }
26 | CenterStrategy::KMeans { num_centers, .. } => format!("{num_centers} centers"),
27 CenterStrategy::UniformGrid { points_per_dim } => {
28 format!("uniform grid with {points_per_dim} points per dimension")
29 }
30 }
31}
32
33pub fn rewrite_thin_plate_knots_error(
34 err: BasisError,
35 termname: &str,
36 feature_count: usize,
37 spec: &ThinPlateBasisSpec,
38) -> BasisError {
39 match err {
40 BasisError::InvalidInput(msg)
43 if msg.contains("thin-plate spline requires at least")
44 && (msg.contains("centers to span") || msg.contains("knots to span")) =>
45 {
46 let min_centers = crate::basis::thin_plate_polynomial_basis_dimension(feature_count);
47 let requested = describe_thin_plate_center_request(&spec.center_strategy);
48 BasisError::InvalidInput(format!(
49 "joint TPS term '{termname}' over {feature_count} covariates with {requested} is invalid; minimum centers is {min_centers}"
50 ))
51 }
52 BasisError::InvalidInput(msg)
57 if msg.starts_with("requested ") && msg.contains(" knots but only ") =>
58 {
59 let min_centers = crate::basis::thin_plate_polynomial_basis_dimension(feature_count);
60 let requested = describe_thin_plate_center_request(&spec.center_strategy);
61 BasisError::InvalidInput(format!(
62 "joint TPS term '{termname}' over {feature_count} covariates with {requested} is invalid; minimum centers is {min_centers}"
63 ))
64 }
65 other => other,
66 }
67}
68
69#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
70pub enum ShapeConstraint {
71 None,
72 MonotoneIncreasing,
73 MonotoneDecreasing,
74 Convex,
75 Concave,
76}
77
78pub fn parse_shape_constraint(raw: &str) -> Result<ShapeConstraint, String> {
89 let normalized = raw.trim().to_ascii_lowercase().replace('-', "_");
90 match normalized.as_str() {
91 "" | "none" => Ok(ShapeConstraint::None),
92 "monotone_increasing" | "monotonic_increasing" | "increasing" | "mono_inc" | "mpi" => {
93 Ok(ShapeConstraint::MonotoneIncreasing)
94 }
95 "monotone_decreasing" | "monotonic_decreasing" | "decreasing" | "mono_dec" | "mpd" => {
96 Ok(ShapeConstraint::MonotoneDecreasing)
97 }
98 "convex" | "cvx" => Ok(ShapeConstraint::Convex),
99 "concave" | "ccv" => Ok(ShapeConstraint::Concave),
100 other => Err(format!(
101 "unknown shape constraint {other:?}; expected one of \
102 \"none\", \"monotone_increasing\", \"monotone_decreasing\", \
103 \"convex\", \"concave\""
104 )),
105 }
106}
107
108impl ShapeConstraint {
109 pub fn dsl_str(&self) -> &'static str {
112 match self {
113 ShapeConstraint::None => "none",
114 ShapeConstraint::MonotoneIncreasing => "monotone_increasing",
115 ShapeConstraint::MonotoneDecreasing => "monotone_decreasing",
116 ShapeConstraint::Convex => "convex",
117 ShapeConstraint::Concave => "concave",
118 }
119 }
120}
121
122pub const SMOOTH_HEAD_KEYWORDS: [&str; 11] = [
125 "s",
126 "smooth",
127 "te",
128 "tensor",
129 "thinplate",
130 "tps",
131 "duchon",
132 "matern",
133 "sphere",
134 "bs",
135 "bspline",
136];
137
138pub fn apply_shape_constraints_to_formula(
151 formula: &str,
152 constraints: &[(String, String)],
153) -> Result<String, String> {
154 use std::collections::{BTreeMap, BTreeSet};
155
156 if constraints.is_empty() {
157 return Ok(formula.to_string());
158 }
159 let strip_ws = |s: &str| -> String { s.chars().filter(|c| !c.is_whitespace()).collect() };
160
161 let mut wanted: BTreeMap<String, &'static str> = BTreeMap::new();
163 let mut originals: BTreeMap<String, String> = BTreeMap::new();
165 for (key, kind_raw) in constraints {
166 let kind = parse_shape_constraint(kind_raw)?;
167 let nk = strip_ws(key);
168 originals.entry(nk.clone()).or_insert_with(|| key.clone());
169 if kind != ShapeConstraint::None {
170 wanted.insert(nk, kind.dsl_str());
171 }
172 }
173 if wanted.is_empty() {
174 return Ok(formula.to_string());
175 }
176
177 let chars: Vec<char> = formula.chars().collect();
178 let n = chars.len();
179 let is_ident = |c: char| c.is_ascii_alphanumeric() || c == '_';
180
181 let mut out = String::with_capacity(formula.len() + 32);
182 let mut matched: BTreeSet<String> = BTreeSet::new();
183 let mut i = 0usize;
184 while i < n {
185 let mut head: Option<(usize, usize)> = None; let mut p = i;
189 while p < n {
190 let boundary = p == 0 || !is_ident(chars[p - 1]);
191 if boundary {
192 for kw in SMOOTH_HEAD_KEYWORDS.iter() {
193 let klen = kw.chars().count();
194 if p + klen > n || chars[p..p + klen].iter().collect::<String>() != **kw {
195 continue;
196 }
197 let mut q = p + klen;
198 while q < n && chars[q].is_whitespace() {
199 q += 1;
200 }
201 if q < n && chars[q] == '(' {
202 head = Some((p, q));
203 break;
204 }
205 }
206 }
207 if head.is_some() {
208 break;
209 }
210 p += 1;
211 }
212 let (head_start, paren_open) = match head {
213 Some(h) => h,
214 None => {
215 out.extend(chars[i..].iter());
216 break;
217 }
218 };
219 out.extend(chars[i..head_start].iter());
220
221 let body_start = paren_open + 1;
223 let mut depth = 1i32;
224 let mut j = body_start;
225 let mut in_str: Option<char> = None;
226 let mut closed = false;
227 while j < n {
228 let ch = chars[j];
229 if let Some(quote) = in_str {
230 if ch == quote {
231 in_str = None;
232 }
233 } else if ch == '\'' || ch == '"' {
234 in_str = Some(ch);
235 } else if ch == '(' {
236 depth += 1;
237 } else if ch == ')' {
238 depth -= 1;
239 if depth == 0 {
240 closed = true;
241 break;
242 }
243 }
244 j += 1;
245 }
246
247 if !closed {
248 out.extend(chars[head_start..].iter());
251 break;
252 }
253
254 let term_text: String = chars[head_start..=j].iter().collect();
255
256 let key_norm = strip_ws(&term_text);
257
258 match wanted.get(&key_norm) {
259 None => out.extend(chars[head_start..=j].iter()),
260 Some(kind) => {
261 let head_paren: String = chars[head_start..body_start].iter().collect();
262 let inside: String = chars[body_start..j].iter().collect();
263 let inside = inside.trim();
264 if inside.is_empty() {
265 out.push_str(&format!("{head_paren}shape={kind})"));
266 } else {
267 out.push_str(&format!("{head_paren}{inside}, shape={kind})"));
268 }
269 matched.insert(key_norm);
270 }
271 }
272
273 i = j + 1;
274 }
275
276 let mut missing: Vec<String> = wanted
277 .keys()
278 .filter(|k| !matched.contains(*k))
279 .map(|k| originals.get(k).cloned().unwrap_or_else(|| k.clone()))
280 .collect();
281
282 if !missing.is_empty() {
283 missing.sort();
284 return Err(format!(
285 "shape constraints referenced smooth term(s) not found in formula: {}",
286 missing.join(", ")
287 ));
288 }
289
290 Ok(out)
291}
292
293#[derive(Debug, Clone, Serialize, Deserialize)]
294pub enum BySmoothKind {
295 Numeric,
296 Level { level_bits: u64 },
297}
298
299#[derive(Debug, Clone, Serialize, Deserialize)]
300pub enum SmoothBasisSpec {
301 ByVariable {
311 inner: Box<SmoothBasisSpec>,
312 by_col: usize,
313 kind: BySmoothKind,
314 by: ByVariableSpec,
315 },
316 FactorSumToZero {
320 inner: Box<SmoothBasisSpec>,
321 by_col: usize,
322 levels: Vec<u64>,
323 #[serde(default)]
334 frozen_global_orthogonality: Option<Array2<f64>>,
335 },
336 BSpline1D {
337 feature_col: usize,
338 spec: BSplineBasisSpec,
339 },
340 BySmooth {
343 smooth: Box<SmoothBasisSpec>,
344 by_kind: ByVarKind,
345 },
346 FactorSmooth { spec: FactorSmoothSpec },
349 ThinPlate {
350 feature_cols: Vec<usize>,
351 spec: ThinPlateBasisSpec,
352 #[serde(default)]
356 input_scales: Option<Vec<f64>>,
357 },
358 Sphere {
359 feature_cols: Vec<usize>,
360 spec: SphericalSplineBasisSpec,
361 },
362 ConstantCurvature {
368 feature_cols: Vec<usize>,
369 spec: ConstantCurvatureBasisSpec,
370 },
371 Matern {
372 feature_cols: Vec<usize>,
373 spec: MaternBasisSpec,
374 #[serde(default)]
375 input_scales: Option<Vec<f64>>,
376 },
377 MeasureJet {
383 feature_cols: Vec<usize>,
384 spec: MeasureJetBasisSpec,
385 #[serde(default)]
386 input_scales: Option<Vec<f64>>,
387 },
388 Duchon {
389 feature_cols: Vec<usize>,
390 spec: DuchonBasisSpec,
391 #[serde(default)]
392 input_scales: Option<Vec<f64>>,
393 },
394 Pca {
395 feature_cols: Vec<usize>,
396 basis_matrix: Array2<f64>,
397 centered: bool,
398 #[serde(default = "default_pca_smooth_penalty")]
399 smooth_penalty: f64,
400 #[serde(default)]
401 center_mean: Option<Array1<f64>>,
402 #[serde(default)]
403 pca_basis_path: Option<PathBuf>,
404 #[serde(default = "default_pca_chunk_size")]
405 chunk_size: usize,
406 },
407 TensorBSpline {
412 feature_cols: Vec<usize>,
413 spec: TensorBSplineSpec,
414 },
415}
416
417impl SmoothBasisSpec {
418 pub fn min_sample_rows(&self) -> usize {
435 const RADIAL_FLOOR: usize = 5;
440
441 match self {
442 Self::ByVariable { inner, .. } => inner.min_sample_rows(),
443 Self::FactorSumToZero { inner, levels, .. } => {
444 let inner_min = inner.min_sample_rows();
448 let lvls = levels.len().saturating_sub(1).max(1);
449 inner_min.saturating_mul(lvls)
450 }
451 Self::BSpline1D { spec, .. } => bspline_basis_min_rows(spec),
452 Self::BySmooth { smooth, .. } => smooth.min_sample_rows(),
453 Self::FactorSmooth { spec } => {
454 bspline_basis_min_rows(&spec.marginal)
458 }
459 Self::ThinPlate { .. }
460 | Self::Sphere { .. }
461 | Self::ConstantCurvature { .. }
462 | Self::Matern { .. }
463 | Self::MeasureJet { .. }
464 | Self::Duchon { .. } => RADIAL_FLOOR,
465 Self::Pca { basis_matrix, .. } => basis_matrix.ncols().max(1),
466 Self::TensorBSpline { spec, .. } => {
467 let mut total: usize = 0;
513 for marginal in &spec.marginalspecs {
514 let m = bspline_basis_min_rows(marginal);
515 total = total.saturating_add(m.max(1));
516 }
517 total.max(RADIAL_FLOOR)
518 }
519 }
520 }
521
522 pub fn structural_kind(&self) -> &'static str {
533 match self {
534 Self::ByVariable { .. } => "by_variable",
535 Self::FactorSumToZero { .. } => "factor_sum_to_zero",
536 Self::BSpline1D { .. } => "bspline_1d",
537 Self::BySmooth { .. } => "by_smooth",
538 Self::FactorSmooth { .. } => "factor_smooth",
539 Self::ThinPlate { .. } => "thin_plate",
540 Self::Sphere { .. } => "sphere",
541 Self::ConstantCurvature { .. } => "constant_curvature",
542 Self::Matern { .. } => "matern",
543 Self::MeasureJet { .. } => "measurejet",
544 Self::Duchon { .. } => "duchon",
545 Self::Pca { .. } => "pca",
546 Self::TensorBSpline { .. } => "tensor_bspline",
547 }
548 }
549
550 pub fn is_marginally_centered_tensor(&self) -> bool {
559 matches!(
560 self,
561 Self::TensorBSpline { spec, .. }
562 if matches!(spec.identifiability, TensorBSplineIdentifiability::MarginalSumToZero)
563 )
564 }
565
566 pub fn is_sum_to_zero_factor_smooth(&self) -> bool {
583 matches!(
584 self,
585 Self::FactorSumToZero { .. }
586 | Self::FactorSmooth {
587 spec: FactorSmoothSpec {
588 flavour: FactorSmoothFlavour::Sz,
589 ..
590 }
591 }
592 )
593 }
594
595 pub fn structural_feature_cols(&self) -> Vec<usize> {
599 match self {
600 Self::ByVariable { inner, .. } | Self::FactorSumToZero { inner, .. } => {
601 inner.structural_feature_cols()
602 }
603 Self::BySmooth { smooth, .. } => smooth.structural_feature_cols(),
604 Self::FactorSmooth { .. } => Vec::new(),
605 Self::BSpline1D { feature_col, .. } => vec![*feature_col],
606 Self::ThinPlate { feature_cols, .. }
607 | Self::Sphere { feature_cols, .. }
608 | Self::ConstantCurvature { feature_cols, .. }
609 | Self::Matern { feature_cols, .. }
610 | Self::MeasureJet { feature_cols, .. }
611 | Self::Duchon { feature_cols, .. }
612 | Self::Pca { feature_cols, .. }
613 | Self::TensorBSpline { feature_cols, .. } => feature_cols.clone(),
614 }
615 }
616}
617
618pub fn bspline_basis_min_rows(spec: &crate::basis::BSplineBasisSpec) -> usize {
643 use crate::basis::BSplineKnotSpec;
644 let columns = match &spec.knotspec {
645 BSplineKnotSpec::Generate {
646 num_internal_knots, ..
647 } => *num_internal_knots + spec.degree + 1,
648 BSplineKnotSpec::Automatic {
649 num_internal_knots: Some(k),
650 ..
651 } => *k + spec.degree + 1,
652 BSplineKnotSpec::Automatic {
653 num_internal_knots: None,
654 ..
655 } => {
656 spec.degree + 2
660 }
661 BSplineKnotSpec::Provided(knots) => knots.len().saturating_sub(spec.degree + 1).max(1),
662 BSplineKnotSpec::NaturalCubicRegression { knots } => knots.len(),
664 BSplineKnotSpec::PeriodicUniform { num_basis, .. } => *num_basis,
665 };
666 let columns = columns.max(spec.degree + 2);
667
668 if spec.double_penalty {
669 const DOUBLE_PENALTY_FLOOR: usize = 2;
672 DOUBLE_PENALTY_FLOOR.min(columns).max(1)
673 } else {
674 columns
675 }
676}
677
678#[derive(Debug, Clone, Serialize, Deserialize)]
679pub enum ByVariableSpec {
680 Numeric,
681 Level { value_bits: u64, label: String },
682}
683
684
685#[derive(Debug, Clone, Serialize, Deserialize)]
686pub enum ByVarKind {
687 Numeric {
688 feature_col: usize,
689 },
690 Factor {
691 feature_col: usize,
692 ordered: bool,
693 frozen_levels: Option<Vec<u64>>,
694 },
695}
696
697#[derive(Debug, Clone, Serialize, Deserialize)]
698pub struct FactorSmoothSpec {
699 pub continuous_cols: Vec<usize>,
700 pub group_col: usize,
701 pub marginal: BSplineBasisSpec,
702 pub flavour: FactorSmoothFlavour,
703 pub group_frozen_levels: Option<Vec<u64>>,
704 #[serde(default)]
710 pub frozen_global_orthogonality: Option<Array2<f64>>,
711}
712
713#[derive(Debug, Clone, Serialize, Deserialize)]
714pub enum FactorSmoothFlavour {
715 Fs { m_null_penalty_orders: Vec<usize> },
716 Sz,
717 Re,
718}
719
720#[derive(Debug, Default, Clone, Serialize, Deserialize)]
721pub struct TensorBSplineSpec {
722 pub marginalspecs: Vec<BSplineBasisSpec>,
723 #[serde(default)]
724 pub periods: Vec<Option<f64>>,
725 pub double_penalty: bool,
726 #[serde(default)]
727 pub identifiability: TensorBSplineIdentifiability,
728 #[serde(default)]
729 pub penalty_decomposition: TensorBSplinePenaltyDecomposition,
730}
731
732#[derive(Debug, Default, Clone, Serialize, Deserialize)]
733pub enum TensorBSplineIdentifiability {
734 None,
735 #[default]
736 SumToZero,
737 MarginalSumToZero,
747 FrozenTransform {
748 transform: Array2<f64>,
749 },
750}
751
752#[derive(Debug, Default, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
753pub enum TensorBSplinePenaltyDecomposition {
754 #[default]
757 MarginalKroneckerSum,
758 Separable,
762}
763
764#[derive(Debug, Clone, Serialize, Deserialize)]
765pub struct SmoothTermSpec {
766 pub name: String,
767 pub basis: SmoothBasisSpec,
768 pub shape: ShapeConstraint,
769 #[serde(default)]
778 pub joint_null_rotation: Option<crate::basis::JointNullRotation>,
779}
780
781#[derive(Debug, Clone)]
782pub struct SmoothTerm {
783 pub name: String,
784 pub coeff_range: Range<usize>,
785 pub shape: ShapeConstraint,
786 pub penalties_local: Vec<Array2<f64>>,
787 pub nullspace_dims: Vec<usize>,
788 pub penaltyinfo_local: Vec<PenaltyInfo>,
789 pub metadata: BasisMetadata,
790 pub lower_bounds_local: Option<Array1<f64>>,
793 pub linear_constraints_local: Option<LinearInequalityConstraints>,
796 pub kronecker_factored: Option<KroneckerFactoredBasis>,
799 pub joint_null_rotation: Option<crate::basis::JointNullRotation>,
822 pub unabsorbed_global_orthogonality: Option<Array2<f64>>,
832}
833
834impl SmoothTerm {
835 pub fn apply_rotation_to_predict(
851 &self,
852 x_new_raw: Array2<f64>,
853 ) -> Result<Array2<f64>, BasisError> {
854 let Some(rot) = self.joint_null_rotation.as_ref() else {
855 return Ok(x_new_raw);
856 };
857 let p_local = rot.rotation.nrows();
858 if x_new_raw.ncols() != p_local {
859 crate::bail_dim_basis!(
860 "joint-null rotation replay for term '{}': raw design has {} columns, \
861 rotation expects {} (the raw basis builder must emit the same column \
862 count as at fit time)",
863 self.name,
864 x_new_raw.ncols(),
865 p_local,
866 );
867 }
868 Ok(gam_linalg::faer_ndarray::fast_ab(
869 &x_new_raw,
870 &rot.rotation,
871 ))
872 }
873
874 pub fn wald_unpenalized_dim(&self) -> usize {
897 joint_unpenalized_dim(
898 self.coeff_range.len(),
899 &self.penalties_local,
900 &self.nullspace_dims,
901 )
902 }
903}
904
905pub fn joint_unpenalized_dim(
910 p_local: usize,
911 penalties_local: &[Array2<f64>],
912 nullspace_dims: &[usize],
913) -> usize {
914 use gam_linalg::faer_ndarray::FaerEigh;
915 if p_local == 0 {
916 return 0;
917 }
918 if penalties_local.is_empty() {
919 return p_local;
921 }
922 let mut s_total = Array2::<f64>::zeros((p_local, p_local));
927 let mut materialized = 0usize;
928 for s in penalties_local {
929 if s.nrows() == p_local && s.ncols() == p_local {
930 s_total += s;
931 materialized += 1;
932 }
933 }
934 if materialized == penalties_local.len() {
935 let symmetric = {
936 let transpose = s_total.t().to_owned();
937 (&s_total + &transpose) * 0.5
938 };
939 if let Ok((evals, _)) = symmetric.eigh(faer::Side::Lower) {
940 let max_abs = evals.iter().fold(0.0_f64, |acc, &v| acc.max(v.abs()));
941 if max_abs == 0.0 {
942 return p_local;
944 }
945 let tol = max_abs * (p_local as f64) * 1e-12;
946 let rank = evals.iter().filter(|&&v| v > tol).count();
947 return p_local.saturating_sub(rank);
948 }
949 }
950 if penalties_local.len() >= 2 {
955 0
956 } else {
957 nullspace_dims
958 .iter()
959 .copied()
960 .min()
961 .unwrap_or(0)
962 .min(p_local)
963 }
964}
965
966#[derive(Debug, Clone, Serialize, Deserialize)]
967pub struct PenaltyBlockInfo {
968 pub global_index: usize,
969 pub termname: Option<String>,
970 pub penalty: PenaltyInfo,
971}
972
973#[derive(Debug, Clone, Serialize, Deserialize)]
974pub struct DroppedPenaltyBlockInfo {
975 pub termname: Option<String>,
976 pub penalty: PenaltyInfo,
977}
978
979#[derive(Debug, Clone)]
980pub struct SmoothDesign {
981 pub term_designs: Vec<DesignMatrix>,
982 pub penalties: Vec<BlockwisePenalty>,
985 pub nullspace_dims: Vec<usize>,
986 pub penaltyinfo: Vec<PenaltyBlockInfo>,
987 pub dropped_penaltyinfo: Vec<DroppedPenaltyBlockInfo>,
988 pub terms: Vec<SmoothTerm>,
989 pub coefficient_lower_bounds: Option<Array1<f64>>,
992 pub linear_constraints: Option<LinearInequalityConstraints>,
995}
996
997impl SmoothDesign {
998 pub fn total_smooth_cols(&self) -> usize {
999 self.term_designs.iter().map(DesignMatrix::ncols).sum()
1000 }
1001 pub fn nrows(&self) -> usize {
1002 self.term_designs.first().map_or(0, DesignMatrix::nrows)
1003 }
1004}
1005
1006#[derive(Debug, Clone)]
1007pub struct RawSmoothDesign {
1008 pub term_designs: Vec<DesignMatrix>,
1009 pub penalties: Vec<BlockwisePenalty>,
1012 pub nullspace_dims: Vec<usize>,
1013 pub penaltyinfo: Vec<PenaltyBlockInfo>,
1014 pub dropped_penaltyinfo: Vec<DroppedPenaltyBlockInfo>,
1015 pub terms: Vec<SmoothTerm>,
1016 pub coefficient_lower_bounds: Option<Array1<f64>>,
1017 pub linear_constraints: Option<LinearInequalityConstraints>,
1018}
1019
1020impl RawSmoothDesign {
1021 pub fn total_smooth_cols(&self) -> usize {
1022 self.term_designs.iter().map(DesignMatrix::ncols).sum()
1023 }
1024 pub fn nrows(&self) -> usize {
1025 self.term_designs.first().map_or(0, DesignMatrix::nrows)
1026 }
1027}
1028
1029impl From<RawSmoothDesign> for SmoothDesign {
1030 fn from(value: RawSmoothDesign) -> Self {
1031 Self {
1032 term_designs: value.term_designs,
1033 penalties: value.penalties,
1034 nullspace_dims: value.nullspace_dims,
1035 penaltyinfo: value.penaltyinfo,
1036 dropped_penaltyinfo: value.dropped_penaltyinfo,
1037 terms: value.terms,
1038 coefficient_lower_bounds: value.coefficient_lower_bounds,
1039 linear_constraints: value.linear_constraints,
1040 }
1041 }
1042}
1043
1044#[derive(Debug, Default, Clone, Serialize, Deserialize)]
1045pub enum BoundedCoefficientPriorSpec {
1046 #[default]
1047 None,
1048 Uniform,
1049 Beta {
1050 a: f64,
1051 b: f64,
1052 },
1053}
1054
1055#[derive(Debug, Clone, Serialize, Deserialize, Default)]
1056pub enum LinearCoefficientGeometry {
1057 #[default]
1058 Unconstrained,
1059 Bounded {
1060 min: f64,
1061 max: f64,
1062 #[serde(default)]
1063 prior: BoundedCoefficientPriorSpec,
1064 },
1065}
1066
1067#[derive(Debug, Clone, Serialize, Deserialize)]
1068pub struct LinearTermSpec {
1069 pub name: String,
1070 pub feature_col: usize,
1076 #[serde(default)]
1079 pub feature_cols: Vec<usize>,
1080 #[serde(default)]
1094 pub categorical_levels: Vec<(usize, u64)>,
1095 #[serde(default = "default_linear_term_double_penalty")]
1102 pub double_penalty: bool,
1103 #[serde(default)]
1104 pub coefficient_geometry: LinearCoefficientGeometry,
1105 #[serde(default)]
1106 pub coefficient_min: Option<f64>,
1107 #[serde(default)]
1108 pub coefficient_max: Option<f64>,
1109}
1110
1111impl LinearTermSpec {
1112 pub fn effective_feature_cols(&self) -> Vec<usize> {
1115 if self.feature_cols.is_empty() {
1116 vec![self.feature_col]
1117 } else {
1118 self.feature_cols.clone()
1119 }
1120 }
1121
1122 pub fn is_interaction(&self) -> bool {
1124 self.feature_cols.len() > 1 || !self.categorical_levels.is_empty()
1125 }
1126
1127 pub fn realized_design_column(&self, data: ArrayView2<'_, f64>) -> Result<Array1<f64>, String> {
1140 let n = data.nrows();
1141 let p = data.ncols();
1142 let bounds = |col: usize| -> Result<(), String> {
1143 if col >= p {
1144 Err(format!(
1145 "linear term '{}' feature column {} out of bounds for {} columns",
1146 self.name, col, p
1147 ))
1148 } else {
1149 Ok(())
1150 }
1151 };
1152
1153 let mut column = if self.categorical_levels.is_empty() {
1158 let cols = self.effective_feature_cols();
1159 for &c in &cols {
1160 bounds(c)?;
1161 }
1162 let mut acc = data.column(cols[0]).to_owned();
1163 for &c in cols.iter().skip(1) {
1164 acc *= &data.column(c);
1165 }
1166 acc
1167 } else {
1168 let mut acc = Array1::<f64>::ones(n);
1169 for &c in &self.feature_cols {
1170 bounds(c)?;
1171 acc *= &data.column(c);
1172 }
1173 acc
1174 };
1175
1176 for &(col, level_bits) in &self.categorical_levels {
1177 bounds(col)?;
1178 let level_bits = gam_data::canonical_level_bits(f64::from_bits(level_bits));
1182 let gate = data.column(col);
1183 for (out, &v) in column.iter_mut().zip(gate.iter()) {
1184 if gam_data::canonical_level_bits(v) != level_bits {
1185 *out = 0.0;
1186 }
1187 }
1188 }
1189
1190 Ok(column)
1191 }
1192}
1193
1194pub const fn default_linear_term_double_penalty() -> bool {
1195 false
1203}
1204
1205pub const fn default_pca_smooth_penalty() -> f64 {
1206 1.0
1207}
1208
1209pub const fn default_pca_chunk_size() -> usize {
1210 4096
1211}
1212
1213#[derive(Debug, Clone, Serialize, Deserialize)]
1219pub struct RandomEffectTermSpec {
1220 pub name: String,
1221 pub feature_col: usize,
1222 pub drop_first_level: bool,
1225 #[serde(default = "default_random_effect_penalized")]
1229 pub penalized: bool,
1230 #[serde(default)]
1233 pub frozen_levels: Option<Vec<u64>>,
1234 #[serde(default = "default_random_effect_lenient_unseen")]
1251 pub lenient_unseen: bool,
1252}
1253
1254pub fn default_random_effect_penalized() -> bool {
1255 true
1256}
1257
1258pub fn default_random_effect_lenient_unseen() -> bool {
1259 true
1260}
1261
1262pub fn validate_measure_jet_positive_vec_len(
1263 label: &str,
1264 term_name: &str,
1265 field: &str,
1266 values: &[f64],
1267 expected: usize,
1268) -> Result<(), String> {
1269 if values.len() != expected {
1270 return Err(SmoothError::invalid_config(format!(
1271 "{label} term '{term_name}' frozen MeasureJet {field} has length {}, expected {expected}",
1272 values.len()
1273 ))
1274 .into());
1275 }
1276 if values
1277 .iter()
1278 .any(|value| !(value.is_finite() && *value > 0.0))
1279 {
1280 return Err(SmoothError::invalid_config(format!(
1281 "{label} term '{term_name}' frozen MeasureJet {field} values must be positive and finite"
1282 ))
1283 .into());
1284 }
1285 Ok(())
1286}
1287
1288#[derive(Debug, Clone, Serialize, Deserialize)]
1289pub struct TermCollectionSpec {
1290 pub linear_terms: Vec<LinearTermSpec>,
1291 pub random_effect_terms: Vec<RandomEffectTermSpec>,
1292 pub smooth_terms: Vec<SmoothTermSpec>,
1293}
1294
1295pub fn validate_smooth_basis_frozen(
1296 basis: &SmoothBasisSpec,
1297 label: &str,
1298 term_name: &str,
1299) -> Result<(), String> {
1300 match basis {
1301 SmoothBasisSpec::ByVariable { inner, .. }
1302 | SmoothBasisSpec::FactorSumToZero { inner, .. } => {
1303 validate_smooth_basis_frozen(inner, label, term_name)
1304 }
1305 SmoothBasisSpec::BSpline1D { spec, .. } => {
1306 if !matches!(
1307 spec.knotspec,
1308 BSplineKnotSpec::Provided(_)
1309 | BSplineKnotSpec::PeriodicUniform { .. }
1310 | BSplineKnotSpec::NaturalCubicRegression { .. }
1311 ) {
1312 return Err(format!(
1313 "{label} term '{term_name}' is not frozen: BSpline knotspec must be Provided, PeriodicUniform, or NaturalCubicRegression"
1314 ));
1315 }
1316 Ok(())
1317 }
1318 SmoothBasisSpec::ThinPlate { spec, .. } => {
1319 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1320 return Err(format!(
1321 "{label} term '{term_name}' is not frozen: ThinPlate centers must be UserProvided"
1322 ));
1323 }
1324 if matches!(
1325 spec.identifiability,
1326 SpatialIdentifiability::OrthogonalToParametric
1327 ) {
1328 return Err(format!(
1329 "{label} term '{term_name}' is not frozen: ThinPlate identifiability must be FrozenTransform or None"
1330 ));
1331 }
1332 Ok(())
1333 }
1334 _ => Ok(()),
1335 }
1336}
1337
1338impl TermCollectionSpec {
1339 pub fn write_structural_shape_hash(&self, h: &mut gam_runtime::warm_start::Fingerprinter) {
1353 h.write_str("term-collection");
1354 h.write_usize(self.linear_terms.len());
1355 for linear in &self.linear_terms {
1356 h.write_str(&linear.name);
1357 }
1358 h.write_usize(self.random_effect_terms.len());
1359 h.write_usize(self.smooth_terms.len());
1360 for smooth in &self.smooth_terms {
1361 h.write_str(&smooth.name);
1362 h.write_str(smooth.basis.structural_kind());
1363 for col in smooth.basis.structural_feature_cols() {
1364 h.write_usize(col);
1365 }
1366 }
1367 }
1368
1369 pub fn validate_frozen(&self, label: &str) -> Result<(), String> {
1373 for linear in &self.linear_terms {
1374 if let (Some(min), Some(max)) = (linear.coefficient_min, linear.coefficient_max)
1375 && (!min.is_finite() || !max.is_finite() || min > max)
1376 {
1377 return Err(SmoothError::invalid_config(format!(
1378 "{label} linear term '{}' has invalid coefficient constraint [{min}, {max}]",
1379 linear.name
1380 ))
1381 .into());
1382 }
1383 if let Some(min) = linear.coefficient_min
1384 && !min.is_finite()
1385 {
1386 return Err(SmoothError::invalid_config(format!(
1387 "{label} linear term '{}' has non-finite coefficient minimum {min}",
1388 linear.name
1389 ))
1390 .into());
1391 }
1392 if let Some(max) = linear.coefficient_max
1393 && !max.is_finite()
1394 {
1395 return Err(SmoothError::invalid_config(format!(
1396 "{label} linear term '{}' has non-finite coefficient maximum {max}",
1397 linear.name
1398 ))
1399 .into());
1400 }
1401 if let LinearCoefficientGeometry::Bounded { min, max, prior } =
1402 &linear.coefficient_geometry
1403 {
1404 if !min.is_finite() || !max.is_finite() || min >= max {
1405 return Err(SmoothError::invalid_config(format!(
1406 "{label} bounded term '{}' has invalid bounds [{min}, {max}]",
1407 linear.name
1408 ))
1409 .into());
1410 }
1411 match prior {
1412 BoundedCoefficientPriorSpec::None | BoundedCoefficientPriorSpec::Uniform => {}
1413 BoundedCoefficientPriorSpec::Beta { a, b } => {
1414 if !a.is_finite() || !b.is_finite() || *a < 1.0 || *b < 1.0 {
1415 return Err(SmoothError::invalid_config(format!(
1416 "{label} bounded term '{}' has invalid Beta prior ({a}, {b})",
1417 linear.name
1418 ))
1419 .into());
1420 }
1421 }
1422 }
1423 }
1424 }
1425 for st in &self.smooth_terms {
1426 match &st.basis {
1427 SmoothBasisSpec::ByVariable { inner, .. } => {
1428 validate_smooth_basis_frozen(inner, label, &st.name)?;
1429 let nested = SmoothTermSpec {
1430 name: st.name.clone(),
1431 basis: (**inner).clone(),
1432 shape: st.shape,
1433 joint_null_rotation: None,
1434 };
1435 TermCollectionSpec {
1436 linear_terms: Vec::new(),
1437 random_effect_terms: Vec::new(),
1438 smooth_terms: vec![nested],
1439 }
1440 .validate_frozen(label)?;
1441 }
1442 SmoothBasisSpec::FactorSumToZero { inner, levels, .. } => {
1443 if levels.len() < 2 {
1444 return Err(format!(
1445 "{label} term '{}' has invalid frozen sz levels",
1446 st.name
1447 ));
1448 }
1449 validate_smooth_basis_frozen(inner, label, &st.name)?;
1450 }
1451 SmoothBasisSpec::BSpline1D { spec, .. } => {
1452 if !matches!(
1453 spec.knotspec,
1454 BSplineKnotSpec::Provided(_)
1455 | BSplineKnotSpec::PeriodicUniform { .. }
1456 | BSplineKnotSpec::NaturalCubicRegression { .. }
1457 ) {
1458 return Err(SmoothError::invalid_config(format!(
1459 "{label} term '{}' is not frozen: BSpline knotspec must be Provided, PeriodicUniform, or NaturalCubicRegression",
1460 st.name
1461 ))
1462 .into());
1463 }
1464 }
1465 SmoothBasisSpec::ThinPlate { spec, .. } => {
1466 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1467 return Err(SmoothError::invalid_config(format!(
1468 "{label} term '{}' is not frozen: ThinPlate centers must be UserProvided",
1469 st.name
1470 ))
1471 .into());
1472 }
1473 if matches!(
1474 spec.identifiability,
1475 SpatialIdentifiability::OrthogonalToParametric
1476 ) {
1477 return Err(SmoothError::invalid_config(format!(
1478 "{label} term '{}' is not frozen: ThinPlate identifiability must be FrozenTransform or None",
1479 st.name
1480 ))
1481 .into());
1482 }
1483 }
1484 SmoothBasisSpec::Sphere { spec, .. } => {
1485 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1486 return Err(SmoothError::invalid_config(format!(
1487 "{label} term '{}' is not frozen: Sphere centers must be UserProvided",
1488 st.name
1489 ))
1490 .into());
1491 }
1492 if matches!(spec.method, crate::basis::SphereMethod::Harmonic)
1493 && spec.max_degree.is_none_or(|d| d == 0)
1494 {
1495 return Err(format!(
1496 "{label} term '{}' is not frozen: sphere max_degree must be positive",
1497 st.name
1498 ));
1499 }
1500 }
1501 SmoothBasisSpec::ConstantCurvature { spec, .. } => {
1502 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1503 return Err(SmoothError::invalid_config(format!(
1504 "{label} term '{}' is not frozen: ConstantCurvature centers must be UserProvided",
1505 st.name
1506 ))
1507 .into());
1508 }
1509 if !(spec.length_scale.is_finite() && spec.length_scale > 0.0) {
1510 return Err(SmoothError::invalid_config(format!(
1511 "{label} term '{}' is not frozen: ConstantCurvature length_scale must be the realized positive value",
1512 st.name
1513 ))
1514 .into());
1515 }
1516 }
1517 SmoothBasisSpec::MeasureJet { spec, .. } => {
1518 let centers = match &spec.center_strategy {
1519 CenterStrategy::UserProvided(centers) => centers,
1520 _ => {
1521 return Err(SmoothError::invalid_config(format!(
1522 "{label} term '{}' is not frozen: MeasureJet centers must be UserProvided",
1523 st.name
1524 ))
1525 .into());
1526 }
1527 };
1528 if centers.nrows() == 0 {
1529 return Err(SmoothError::invalid_config(format!(
1530 "{label} term '{}' is not frozen: MeasureJet centers are empty",
1531 st.name
1532 ))
1533 .into());
1534 }
1535 if !(spec.length_scale.is_finite() && spec.length_scale > 0.0) {
1536 return Err(SmoothError::invalid_config(format!(
1537 "{label} term '{}' is not frozen: MeasureJet length_scale must be the realized positive value",
1538 st.name
1539 ))
1540 .into());
1541 }
1542 let frozen = spec.frozen_quadrature.as_ref().ok_or_else(|| {
1545 SmoothError::invalid_config(format!(
1546 "{label} term '{}' is not frozen: MeasureJet frozen_quadrature payload is missing",
1547 st.name
1548 ))
1549 })?;
1550 if frozen.masses.len() != centers.nrows() {
1551 return Err(SmoothError::invalid_config(format!(
1552 "{label} term '{}' frozen MeasureJet has {} masses for {} centers",
1553 st.name,
1554 frozen.masses.len(),
1555 centers.nrows()
1556 ))
1557 .into());
1558 }
1559 let total_mass = frozen.masses.sum();
1560 if frozen
1561 .masses
1562 .iter()
1563 .any(|mass| !(mass.is_finite() && *mass >= 0.0))
1564 || !(total_mass.is_finite() && total_mass > 0.0)
1565 {
1566 return Err(SmoothError::invalid_config(format!(
1567 "{label} term '{}' frozen MeasureJet masses must be finite, nonnegative, and have positive total mass",
1568 st.name
1569 ))
1570 .into());
1571 }
1572 let n_levels = frozen.eps_band.len();
1573 if n_levels == 0
1574 || frozen
1575 .eps_band
1576 .iter()
1577 .any(|eps| !(eps.is_finite() && *eps > 0.0))
1578 {
1579 return Err(SmoothError::invalid_config(format!(
1580 "{label} term '{}' frozen MeasureJet eps_band must be nonempty, finite, and positive",
1581 st.name
1582 ))
1583 .into());
1584 }
1585 for (idx, pair) in frozen.eps_band.windows(2).enumerate() {
1586 if pair[1] <= pair[0] {
1587 return Err(SmoothError::invalid_config(format!(
1588 "{label} term '{}' frozen MeasureJet eps_band is not strictly ascending at {idx}: {} then {}",
1589 st.name,
1590 pair[0],
1591 pair[1]
1592 ))
1593 .into());
1594 }
1595 }
1596 validate_measure_jet_positive_vec_len(
1597 label,
1598 &st.name,
1599 "support_means",
1600 &frozen.support_means,
1601 n_levels,
1602 )?;
1603 let per_level = crate::basis::measure_jet_multiscale_mode(spec);
1611 if per_level {
1612 validate_measure_jet_positive_vec_len(
1613 label,
1614 &st.name,
1615 "penalty_normalization_scales",
1616 &frozen.penalty_normalization_scales,
1617 n_levels,
1618 )?;
1619 validate_measure_jet_positive_vec_len(
1620 label,
1621 &st.name,
1622 "raw_penalty_normalization_scales",
1623 &frozen.raw_penalty_normalization_scales,
1624 n_levels,
1625 )?;
1626 if frozen.fused_penalty_normalization_scale.is_some() {
1627 return Err(SmoothError::invalid_config(format!(
1628 "{label} term '{}' per-level MeasureJet must not carry a fused penalty normalization scale",
1629 st.name
1630 ))
1631 .into());
1632 }
1633 } else {
1634 if !frozen.penalty_normalization_scales.is_empty()
1635 || !frozen.raw_penalty_normalization_scales.is_empty()
1636 {
1637 return Err(SmoothError::invalid_config(format!(
1638 "{label} term '{}' fused MeasureJet must not carry per-level penalty normalization scales",
1639 st.name
1640 ))
1641 .into());
1642 }
1643 match frozen.fused_penalty_normalization_scale {
1644 Some(scale) if scale.is_finite() && scale > 0.0 => {}
1645 Some(scale) => {
1646 return Err(SmoothError::invalid_config(format!(
1647 "{label} term '{}' fused MeasureJet penalty normalization scale must be positive and finite, got {scale}",
1648 st.name
1649 ))
1650 .into());
1651 }
1652 None => {
1653 return Err(SmoothError::invalid_config(format!(
1654 "{label} term '{}' fused MeasureJet is missing its penalty normalization scale",
1655 st.name
1656 ))
1657 .into());
1658 }
1659 }
1660 }
1661 }
1662 SmoothBasisSpec::Matern { spec, .. } => {
1663 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1664 return Err(SmoothError::invalid_config(format!(
1665 "{label} term '{}' is not frozen: Matern centers must be UserProvided",
1666 st.name
1667 ))
1668 .into());
1669 }
1670 }
1671 SmoothBasisSpec::Duchon { spec, .. } => {
1672 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1673 return Err(SmoothError::invalid_config(format!(
1674 "{label} term '{}' is not frozen: Duchon centers must be UserProvided",
1675 st.name
1676 ))
1677 .into());
1678 }
1679 if matches!(
1680 spec.identifiability,
1681 SpatialIdentifiability::OrthogonalToParametric
1682 ) {
1683 return Err(SmoothError::invalid_config(format!(
1684 "{label} term '{}' is not frozen: Duchon identifiability must be FrozenTransform or None",
1685 st.name
1686 ))
1687 .into());
1688 }
1689 }
1690 SmoothBasisSpec::Pca {
1691 centered,
1692 center_mean,
1693 pca_basis_path,
1694 ..
1695 } => {
1696 if *centered && center_mean.is_none() && pca_basis_path.is_none() {
1697 return Err(SmoothError::invalid_config(format!(
1698 "{label} term '{}' is not frozen: centered Pca missing center_mean",
1699 st.name
1700 ))
1701 .into());
1702 }
1703 }
1704 SmoothBasisSpec::BySmooth { smooth, by_kind } => {
1705 if let SmoothBasisSpec::BySmooth { .. } = smooth.as_ref() {
1706 return Err(format!("{label} term '{}' has nested by-smooths", st.name));
1707 }
1708 match by_kind {
1709 ByVarKind::Numeric { .. } => {}
1710 ByVarKind::Factor { frozen_levels, .. } if frozen_levels.is_none() => {
1711 return Err(format!(
1712 "{label} term '{}' is not frozen: by-factor levels missing",
1713 st.name
1714 ));
1715 }
1716 ByVarKind::Factor { .. } => {}
1717 }
1718 let nested = TermCollectionSpec {
1719 linear_terms: vec![],
1720 random_effect_terms: vec![],
1721 smooth_terms: vec![SmoothTermSpec {
1722 name: st.name.clone(),
1723 basis: (**smooth).clone(),
1724 shape: st.shape,
1725 joint_null_rotation: None,
1726 }],
1727 };
1728 nested.validate_frozen(label)?;
1729 }
1730 SmoothBasisSpec::FactorSmooth { spec } => {
1731 if spec.group_frozen_levels.is_none() {
1732 return Err(format!(
1733 "{label} term '{}' is not frozen: factor-smooth levels missing",
1734 st.name
1735 ));
1736 }
1737 if !matches!(
1738 spec.marginal.knotspec,
1739 BSplineKnotSpec::Provided(_)
1740 | BSplineKnotSpec::PeriodicUniform { .. }
1741 | BSplineKnotSpec::NaturalCubicRegression { .. }
1753 ) {
1754 return Err(format!(
1755 "{label} term '{}' is not frozen: factor-smooth marginal knots missing",
1756 st.name
1757 ));
1758 }
1759 }
1760 SmoothBasisSpec::TensorBSpline { spec, .. } => {
1761 for (dim, marginal) in spec.marginalspecs.iter().enumerate() {
1762 if !matches!(
1763 marginal.knotspec,
1764 BSplineKnotSpec::Provided(_)
1765 | BSplineKnotSpec::PeriodicUniform { .. }
1766 | BSplineKnotSpec::NaturalCubicRegression { .. }
1767 ) {
1768 return Err(SmoothError::invalid_config(format!(
1769 "{label} term '{}' dim {} is not frozen: tensor marginal knotspec must be Provided, PeriodicUniform, or NaturalCubicRegression",
1770 st.name, dim
1771 ))
1772 .into());
1773 }
1774 }
1775 if matches!(
1776 spec.identifiability,
1777 TensorBSplineIdentifiability::SumToZero
1778 | TensorBSplineIdentifiability::MarginalSumToZero
1779 ) {
1780 return Err(SmoothError::invalid_config(format!(
1781 "{label} term '{}' is not frozen: tensor identifiability must be FrozenTransform or None",
1782 st.name
1783 ))
1784 .into());
1785 }
1786 }
1787 }
1788 }
1789
1790 for rt in &self.random_effect_terms {
1791 if rt.frozen_levels.is_none() {
1792 return Err(SmoothError::invalid_config(format!(
1793 "{label} random-effect term '{}' is not frozen: missing frozen_levels",
1794 rt.name
1795 ))
1796 .into());
1797 }
1798 }
1799
1800 Ok(())
1801 }
1802
1803 pub fn remap_feature_columns<E, F>(&self, mut remap: F) -> Result<TermCollectionSpec, E>
1822 where
1823 F: FnMut(usize) -> Result<usize, E>,
1824 {
1825 let mut out = self.clone();
1826 for lt in &mut out.linear_terms {
1827 lt.feature_col = remap(lt.feature_col)?;
1828 for fc in lt.feature_cols.iter_mut() {
1838 *fc = remap(*fc)?;
1839 }
1840 for (col, _bits) in lt.categorical_levels.iter_mut() {
1845 *col = remap(*col)?;
1846 }
1847 }
1848 for rt in &mut out.random_effect_terms {
1849 rt.feature_col = remap(rt.feature_col)?;
1850 }
1851 for st in &mut out.smooth_terms {
1852 remap_smooth_basis_feature_columns(&mut st.basis, &mut remap)?;
1853 }
1854 Ok(out)
1855 }
1856}
1857
1858pub fn remap_smooth_basis_feature_columns<E, F>(
1863 basis: &mut SmoothBasisSpec,
1864 remap: &mut F,
1865) -> Result<(), E>
1866where
1867 F: FnMut(usize) -> Result<usize, E>,
1868{
1869 match basis {
1870 SmoothBasisSpec::ByVariable { inner, by_col, .. }
1871 | SmoothBasisSpec::FactorSumToZero { inner, by_col, .. } => {
1872 *by_col = remap(*by_col)?;
1873 remap_smooth_basis_feature_columns(inner, remap)?;
1874 }
1875 SmoothBasisSpec::BSpline1D { feature_col, .. } => {
1876 *feature_col = remap(*feature_col)?;
1877 }
1878 SmoothBasisSpec::BySmooth { smooth, by_kind } => {
1879 let by_feature_col = match by_kind {
1880 ByVarKind::Numeric { feature_col } | ByVarKind::Factor { feature_col, .. } => {
1881 feature_col
1882 }
1883 };
1884 *by_feature_col = remap(*by_feature_col)?;
1885 remap_smooth_basis_feature_columns(smooth, remap)?;
1886 }
1887 SmoothBasisSpec::FactorSmooth { spec } => {
1888 for fc in spec.continuous_cols.iter_mut() {
1889 *fc = remap(*fc)?;
1890 }
1891 spec.group_col = remap(spec.group_col)?;
1892 }
1893 SmoothBasisSpec::ThinPlate { feature_cols, .. }
1894 | SmoothBasisSpec::Sphere { feature_cols, .. }
1895 | SmoothBasisSpec::ConstantCurvature { feature_cols, .. }
1896 | SmoothBasisSpec::Matern { feature_cols, .. }
1897 | SmoothBasisSpec::MeasureJet { feature_cols, .. }
1898 | SmoothBasisSpec::Duchon { feature_cols, .. }
1899 | SmoothBasisSpec::Pca { feature_cols, .. }
1900 | SmoothBasisSpec::TensorBSpline { feature_cols, .. } => {
1901 for fc in feature_cols.iter_mut() {
1902 *fc = remap(*fc)?;
1903 }
1904 }
1905 }
1906 Ok(())
1907}
1908
1909#[derive(Debug, Clone)]
1910pub enum PenaltyStructureHint {
1911 Ridge(f64),
1912 Kronecker(Vec<Array2<f64>>),
1913}
1914
1915#[derive(Clone)]
1922pub struct BlockwisePenalty {
1923 pub col_range: Range<usize>,
1925 pub local: Array2<f64>,
1928 pub prior_mean: gam_problem::CoefficientPriorMean,
1930 pub structure_hint: Option<PenaltyStructureHint>,
1933 pub op: Option<std::sync::Arc<dyn crate::analytic_penalties::PenaltyOp>>,
1938}
1939
1940impl std::fmt::Debug for BlockwisePenalty {
1941 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
1942 f.debug_struct("BlockwisePenalty")
1943 .field("col_range", &self.col_range)
1944 .field(
1945 "local",
1946 &format_args!("{}×{}", self.local.nrows(), self.local.ncols()),
1947 )
1948 .field("prior_mean", &self.prior_mean)
1949 .field("structure_hint", &self.structure_hint)
1950 .field("op", &self.op.as_ref().map(|o| o.dim()))
1951 .finish()
1952 }
1953}
1954
1955impl BlockwisePenalty {
1956 pub fn new(col_range: Range<usize>, local: Array2<f64>) -> Self {
1958 assert_eq!(col_range.len(), local.nrows());
1959 assert_eq!(col_range.len(), local.ncols());
1960 Self {
1961 col_range,
1962 local,
1963 prior_mean: gam_problem::CoefficientPriorMean::Zero,
1964 structure_hint: None,
1965 op: None,
1966 }
1967 }
1968
1969 pub fn with_prior_mean(
1970 mut self,
1971 prior_mean: gam_problem::CoefficientPriorMean,
1972 ) -> Self {
1973 self.prior_mean = prior_mean;
1974 self
1975 }
1976
1977 pub fn with_op(
1979 mut self,
1980 op: Option<std::sync::Arc<dyn crate::analytic_penalties::PenaltyOp>>,
1981 ) -> Self {
1982 self.op = op;
1983 self
1984 }
1985
1986 pub fn ridge(col_range: Range<usize>, scale: f64) -> Self {
1987 let block_size = col_range.len();
1988 let mut local = Array2::<f64>::zeros((block_size, block_size));
1989 for i in 0..block_size {
1990 local[[i, i]] = scale;
1991 }
1992 Self {
1993 col_range,
1994 local,
1995 prior_mean: gam_problem::CoefficientPriorMean::Zero,
1996 structure_hint: Some(PenaltyStructureHint::Ridge(scale)),
1997 op: None,
1998 }
1999 }
2000
2001 pub fn kronecker(
2002 col_range: Range<usize>,
2003 local: Array2<f64>,
2004 factors: Vec<Array2<f64>>,
2005 ) -> Self {
2006 assert_eq!(col_range.len(), local.nrows());
2007 assert_eq!(col_range.len(), local.ncols());
2008 Self {
2009 col_range,
2010 local,
2011 prior_mean: gam_problem::CoefficientPriorMean::Zero,
2012 structure_hint: Some(PenaltyStructureHint::Kronecker(factors)),
2013 op: None,
2014 }
2015 }
2016
2017 pub fn to_global(&self, p_total: usize) -> Array2<f64> {
2021 let mut g = Array2::<f64>::zeros((p_total, p_total));
2022 let r = &self.col_range;
2023 assert!(
2024 r.end <= p_total && self.local.nrows() == r.len() && self.local.ncols() == r.len(),
2025 "BlockwisePenalty::to_global shape invariant violated: \
2026 col_range={}..{}, local={}x{}, p_total={}",
2027 r.start,
2028 r.end,
2029 self.local.nrows(),
2030 self.local.ncols(),
2031 p_total,
2032 );
2033 g.slice_mut(s![r.start..r.end, r.start..r.end])
2034 .assign(&self.local);
2035 g
2036 }
2037
2038 pub fn to_penalty_matrix(
2041 &self,
2042 total_dim: usize,
2043 ) -> gam_problem::PenaltyMatrix {
2044 gam_problem::PenaltyMatrix::Blockwise {
2045 local: self.local.clone(),
2046 col_range: self.col_range.clone(),
2047 total_dim,
2048 }
2049 }
2050
2051 #[inline]
2053 pub fn block_size(&self) -> usize {
2054 self.col_range.len()
2055 }
2056}
2057
2058pub fn weighted_blockwise_penalty_sum(
2062 penalties: &[BlockwisePenalty],
2063 lambdas: &[f64],
2064 p_total: usize,
2065) -> Array2<f64> {
2066 assert_eq!(penalties.len(), lambdas.len());
2067 for (idx, &lam) in lambdas.iter().enumerate() {
2074 assert!(
2075 lam.is_finite() && lam >= 0.0,
2076 "weighted_blockwise_penalty_sum: lambdas[{idx}] = {lam} is invalid (must be finite and non-negative; negative smoothing parameters violate S_λ ⪰ 0)",
2077 );
2078 }
2079 for (idx, bp) in penalties.iter().enumerate() {
2083 let r = &bp.col_range;
2084 assert!(
2085 r.end <= p_total,
2086 "weighted_blockwise_penalty_sum: penalties[{idx}] col_range {:?} exceeds p_total = {p_total}",
2087 r,
2088 );
2089 }
2090 let mut out = Array2::<f64>::zeros((p_total, p_total));
2091 for (bp, &lam) in penalties.iter().zip(lambdas.iter()) {
2092 let r = &bp.col_range;
2093 let mut slice = out.slice_mut(s![r.start..r.end, r.start..r.end]);
2094 slice.scaled_add(lam, &bp.local);
2095 }
2096 out
2097}
2098
2099#[derive(Debug, Clone)]
2106pub struct KroneckerPenaltySystem {
2107 pub marginal_penalties: Vec<Array2<f64>>,
2109 pub marginal_eigensystems: Vec<(Array1<f64>, Array2<f64>)>,
2111 pub marginal_dims: Vec<usize>,
2113 pub has_double_penalty: bool,
2115}
2116
2117impl KroneckerPenaltySystem {
2118 pub fn new(
2119 marginal_penalties: Vec<Array2<f64>>,
2120 marginal_dims: Vec<usize>,
2121 has_double_penalty: bool,
2122 ) -> Result<Self, BasisError> {
2123 if marginal_penalties.len() != marginal_dims.len() {
2124 crate::bail_dim_basis!(
2125 "KroneckerPenaltySystem: {} penalties vs {} dims",
2126 marginal_penalties.len(),
2127 marginal_dims.len()
2128 );
2129 }
2130 let eigensystems =
2131 kronecker_marginal_eigensystems(&marginal_penalties, "KroneckerPenaltySystem")
2132 .map_err(|e| BasisError::InvalidInput(e.to_string()))?;
2133 Ok(Self {
2134 marginal_penalties,
2135 marginal_eigensystems: eigensystems,
2136 marginal_dims,
2137 has_double_penalty,
2138 })
2139 }
2140
2141 pub fn p_total(&self) -> usize {
2142 self.marginal_dims.iter().copied().product()
2143 }
2144
2145 pub fn ndim(&self) -> usize {
2146 self.marginal_dims.len()
2147 }
2148
2149 pub fn num_penalties(&self) -> usize {
2150 self.marginal_dims.len() + if self.has_double_penalty { 1 } else { 0 }
2151 }
2152
2153 pub fn logdet_and_derivatives(
2157 &self,
2158 lambdas: &[f64],
2159 ridge: f64,
2160 ) -> (f64, Array1<f64>, Array2<f64>) {
2161 let n_pen = self.num_penalties();
2162 assert_eq!(lambdas.len(), n_pen, "lambda count mismatch");
2163 let marginal_evals: Vec<_> = self
2164 .marginal_eigensystems
2165 .iter()
2166 .map(|(evals, _)| evals.view())
2167 .collect();
2168 kronecker_logdet_and_derivatives(
2169 &marginal_evals,
2170 &self.marginal_dims,
2171 lambdas,
2172 self.has_double_penalty,
2173 ridge,
2174 )
2175 }
2176
2177 pub fn logdet_rank_and_derivatives(
2178 &self,
2179 lambdas: &[f64],
2180 ridge: f64,
2181 ) -> (f64, usize, Array1<f64>, Array2<f64>) {
2182 let n_pen = self.num_penalties();
2183 assert_eq!(lambdas.len(), n_pen, "lambda count mismatch");
2184 let d = self.marginal_dims.len();
2185 let mut logdet = 0.0;
2186 let mut rank = 0usize;
2187 let mut grad = Array1::<f64>::zeros(n_pen);
2188 let mut hess = Array2::<f64>::zeros((n_pen, n_pen));
2189 const EIGENVALUE_POSITIVITY_FLOOR: f64 = 1e-12;
2193 const STRUCTURAL_ZERO_FLOOR: f64 = 1e-12;
2197 let mut multi_idx = vec![0usize; d];
2198 loop {
2199 let mut sigma = 0.0;
2200 let mut structural_sigma = 0.0;
2201 for k in 0..d {
2202 let marginal_eigenvalue = self.marginal_eigensystems[k].0[multi_idx[k]];
2203 structural_sigma += marginal_eigenvalue;
2204 sigma += lambdas[k] * marginal_eigenvalue;
2205 }
2206 let joint_null = structural_sigma <= STRUCTURAL_ZERO_FLOOR;
2207 if self.has_double_penalty && joint_null {
2208 sigma += lambdas[d];
2209 }
2210 if structural_sigma > STRUCTURAL_ZERO_FLOOR {
2211 sigma += ridge;
2212 }
2213
2214 if sigma > EIGENVALUE_POSITIVITY_FLOOR {
2215 rank += 1;
2216 logdet += sigma.ln();
2217 let inv_sigma = 1.0 / sigma;
2218 let inv_sigma2 = inv_sigma * inv_sigma;
2219 for k in 0..n_pen {
2220 let ck = if k < d {
2221 lambdas[k] * self.marginal_eigensystems[k].0[multi_idx[k]]
2222 } else if joint_null {
2223 lambdas[d]
2224 } else {
2225 0.0
2226 };
2227 grad[k] += ck * inv_sigma;
2228 hess[[k, k]] += ck * inv_sigma - ck * ck * inv_sigma2;
2229 for l in (k + 1)..n_pen {
2230 let cl = if l < d {
2231 lambdas[l] * self.marginal_eigensystems[l].0[multi_idx[l]]
2232 } else if joint_null {
2233 lambdas[d]
2234 } else {
2235 0.0
2236 };
2237 let off = -ck * cl * inv_sigma2;
2238 hess[[k, l]] += off;
2239 hess[[l, k]] += off;
2240 }
2241 }
2242 }
2243
2244 let mut carry = true;
2245 for dim in (0..d).rev() {
2246 if carry {
2247 multi_idx[dim] += 1;
2248 if multi_idx[dim] < self.marginal_dims[dim] {
2249 carry = false;
2250 } else {
2251 multi_idx[dim] = 0;
2252 }
2253 }
2254 }
2255 if carry {
2256 break;
2257 }
2258 }
2259 (logdet, rank, grad, hess)
2260 }
2261}
2262
2263#[cfg(test)]
2264mod joint_unpenalized_dim_tests {
2265 use super::joint_unpenalized_dim;
2266 use ndarray::{Array2, array};
2267
2268 #[test]
2269 fn no_penalty_is_fully_unpenalized() {
2270 assert_eq!(joint_unpenalized_dim(4, &[], &[]), 4);
2271 }
2272
2273 #[test]
2274 fn single_penalty_returns_its_own_null_space() {
2275 let s = array![[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 5.0]];
2278 assert_eq!(joint_unpenalized_dim(3, std::slice::from_ref(&s), &[2]), 2);
2279 }
2280
2281 #[test]
2282 fn complementary_double_penalty_has_empty_joint_null_space() {
2283 let bending = array![[0.0, 0.0, 0.0], [0.0, 4.0, 0.0], [0.0, 0.0, 4.0]];
2290 let ridge = array![[2.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]];
2291 assert_eq!(joint_unpenalized_dim(3, &[bending, ridge], &[1, 2]), 0);
2292 }
2293
2294 #[test]
2295 fn partial_overlap_keeps_shared_null_direction() {
2296 let a = array![[0.0, 0.0, 0.0], [0.0, 3.0, 0.0], [0.0, 0.0, 0.0]];
2300 let b = array![[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 3.0]];
2301 assert_eq!(joint_unpenalized_dim(3, &[a, b], &[2, 2]), 1);
2302 }
2303
2304 #[test]
2305 fn non_materialized_penalty_falls_back_conservatively() {
2306 let full: Array2<f64> = array![[0.0, 0.0], [0.0, 1.0]];
2310 let factor: Array2<f64> = array![[1.0]]; assert_eq!(
2312 joint_unpenalized_dim(2, &[full, factor.clone()], &[1, 0]),
2313 0
2314 );
2315 assert_eq!(joint_unpenalized_dim(4, std::slice::from_ref(&factor), &[2]), 2);
2317 }
2318}
2319
2320#[cfg(test)]
2321mod kronecker_penalty_system_tests {
2322 use super::KroneckerPenaltySystem;
2323 use ndarray::array;
2324
2325 #[test]
2326 fn double_penalty_rank_derivatives_use_only_joint_null_space() {
2327 let penalties = vec![
2328 array![[0.0, 0.0], [0.0, 2.0]],
2329 array![[0.0, 0.0], [0.0, 3.0]],
2330 ];
2331 let system = KroneckerPenaltySystem::new(penalties, vec![2usize, 2usize], true).unwrap();
2332 let lambdas = vec![5.0, 7.0, 11.0];
2333
2334 let (logdet, rank, grad, hess) = system.logdet_rank_and_derivatives(&lambdas, 0.0);
2335
2336 let expected_diag = [11.0_f64, 21.0, 10.0, 31.0];
2337 let expected_logdet: f64 = expected_diag.iter().map(|v| v.ln()).sum();
2338 assert_eq!(rank, 4);
2339 assert!((logdet - expected_logdet).abs() <= 1e-12);
2340 assert!(
2341 (grad[2] - 1.0).abs() <= 1e-12,
2342 "double-penalty rank derivative must count only the joint null mode, got {}",
2343 grad[2]
2344 );
2345 assert!(hess[[2, 2]].abs() <= 1e-12);
2346 }
2347}
2348
2349#[derive(Clone, Debug)]
2350pub struct TermCollectionDesign {
2351 pub design: DesignMatrix,
2360 pub penalties: Vec<BlockwisePenalty>,
2361 pub nullspace_dims: Vec<usize>,
2362 pub penaltyinfo: Vec<PenaltyBlockInfo>,
2363 pub dropped_penaltyinfo: Vec<DroppedPenaltyBlockInfo>,
2364 pub coefficient_lower_bounds: Option<Array1<f64>>,
2367 pub linear_constraints: Option<LinearInequalityConstraints>,
2370 pub intercept_range: Range<usize>,
2371 pub linear_ranges: Vec<(String, Range<usize>)>,
2372 pub random_effect_ranges: Vec<(String, Range<usize>)>,
2373 pub random_effect_levels: Vec<(String, Vec<u64>)>,
2374 pub smooth: SmoothDesign,
2375}
2376
2377impl TermCollectionDesign {
2378 pub fn leading_penalty_blocks_before_smooth(&self) -> usize {
2386 self.penaltyinfo
2387 .iter()
2388 .take_while(|info| {
2389 matches!(
2390 &info.penalty.source,
2391 crate::basis::PenaltySource::Other(source)
2392 if source == "LinearTermRidge"
2393 || source.starts_with("RandomEffectRidge(")
2394 )
2395 })
2396 .count()
2397 }
2398
2399 pub fn penalties_as_penalty_matrix(&self) -> Vec<gam_problem::PenaltyMatrix> {
2403 let p = self.design.ncols();
2404 self.penalties
2405 .iter()
2406 .map(|bp| bp.to_penalty_matrix(p))
2407 .collect()
2408 }
2409
2410 #[inline]
2412 pub fn num_penalties(&self) -> usize {
2413 self.penalties.len()
2414 }
2415
2416 pub fn realize_coefficient_groups(
2419 &self,
2420 groups: &[CoefficientGroupSpec],
2421 base_prior: &gam_spec::RhoPrior,
2422 ) -> Result<RealizedCoefficientGroups, BasisError> {
2423 realize_coefficient_groups(self, groups, base_prior)
2424 }
2425
2426 pub fn kronecker_penalty_system(&self) -> Option<KroneckerPenaltySystem> {
2437 let [only_term] = self.smooth.terms.as_slice() else {
2438 return None;
2439 };
2440 let kron = only_term.kronecker_factored.as_ref()?;
2441 if kron.marginal_dims.len() < 2
2447 || kron.marginal_penalties.len() != kron.marginal_dims.len()
2448 || kron.marginal_designs.len() != kron.marginal_dims.len()
2449 {
2450 return None;
2451 }
2452 KroneckerPenaltySystem::new(
2453 kron.marginal_penalties.clone(),
2454 kron.marginal_dims.clone(),
2455 kron.has_double_penalty,
2456 )
2457 .ok()
2458 }
2459}
2460
2461#[derive(Clone)]
2467pub struct StandardLatentCoordConfig {
2468 pub values: std::sync::Arc<crate::latent::LatentCoordValues>,
2469 pub term_index: gam_problem::types::SmoothTermIdx,
2470 pub feature_cols: Vec<usize>,
2471 pub manifold: crate::latent::LatentManifold,
2472 pub manifold_auto: bool,
2473 pub retraction_registry: gam_problem::LatentRetractionRegistry,
2474 pub analytic_penalties: Option<std::sync::Arc<crate::AnalyticPenaltyRegistry>>,
2475}
2476
2477#[derive(Clone, Debug, Serialize, Deserialize)]
2478pub struct AdaptiveSpatialMap {
2479 pub termname: String,
2480 pub feature_cols: Vec<usize>,
2481 pub collocation_points: Array2<f64>,
2482 pub inv_magweight: Array1<f64>,
2483 pub invgradweight: Array1<f64>,
2484 pub inv_lapweight: Array1<f64>,
2485}
2486
2487#[derive(Clone, Debug, Serialize, Deserialize)]
2488pub struct AdaptiveRegularizationDiagnostics {
2489 pub epsilon_0: f64,
2490 pub epsilon_g: f64,
2491 pub epsilon_c: f64,
2492 pub epsilon_outer_iterations: usize,
2493 pub mm_iterations: usize,
2494 pub converged: bool,
2495 pub maps: Vec<AdaptiveSpatialMap>,
2496}
2497
2498#[derive(Debug, Clone)]
2499pub struct LinearColumnConditioning {
2500 col_idx: usize,
2501 mean: f64,
2502 scale: f64,
2503}
2504
2505#[derive(Debug, Clone, Default)]
2506pub struct LinearFitConditioning {
2507 pub intercept_idx: usize,
2508 pub columns: Vec<LinearColumnConditioning>,
2509}
2510
2511#[derive(Clone)]
2512pub struct SpatialPsiDerivative {
2513 pub penalty_index: usize,
2515 pub penalty_indices: Vec<usize>,
2516 pub global_range: Range<usize>,
2517 pub total_p: usize,
2518 pub x_psi_local: Array2<f64>,
2519 pub s_psi_components_local: Vec<Array2<f64>>,
2520 pub x_psi_psi_local: Array2<f64>,
2521 pub s_psi_psi_components_local: Vec<Array2<f64>>,
2522 pub aniso_group_id: Option<usize>,
2523 pub aniso_cross_designs: Option<Vec<(usize, Array2<f64>)>>,
2526 pub aniso_cross_penalty_provider: Option<
2530 std::sync::Arc<
2531 dyn Fn(usize) -> Result<Vec<Array2<f64>>, EstimationError> + Send + Sync + 'static,
2532 >,
2533 >,
2534 pub implicit_operator: Option<std::sync::Arc<crate::basis::ImplicitDesignPsiDerivative>>,
2539 pub implicit_axis: usize,
2541}
2542
2543#[derive(Debug, Clone)]
2544pub struct SpatialLogKappaCoords {
2545 pub values: Array1<f64>,
2548 pub dims_per_term: Vec<usize>,
2550}
2551
2552#[derive(Clone, Copy)]
2557pub enum AnisoBoundEnd {
2558 Lower,
2559 Upper,
2560}
2561
2562impl SpatialLogKappaCoords {
2563 pub fn new_with_dims(values: Array1<f64>, dims_per_term: Vec<usize>) -> Self {
2565 assert_eq!(
2566 values.len(),
2567 dims_per_term.iter().sum::<usize>(),
2568 "SpatialLogKappaCoords: values length {} != sum of dims_per_term {}",
2569 values.len(),
2570 dims_per_term.iter().sum::<usize>(),
2571 );
2572 Self {
2573 values,
2574 dims_per_term,
2575 }
2576 }
2577
2578 pub fn from_length_scales(
2580 spec: &TermCollectionSpec,
2581 term_indices: &[usize],
2582 options: &SpatialLengthScaleOptimizationOptions,
2583 ) -> Self {
2584 let mut out = Array1::<f64>::zeros(term_indices.len());
2585 for (slot, &term_idx) in term_indices.iter().enumerate() {
2586 if let Some(cc) = constant_curvature_term_spec(spec, term_idx) {
2592 out[slot] = cc.kappa;
2593 continue;
2594 }
2595 let length_scale = get_spatial_length_scale(spec, term_idx)
2596 .unwrap_or(options.min_length_scale)
2597 .clamp(options.min_length_scale, options.max_length_scale);
2598 out[slot] = -length_scale.ln();
2599 }
2600 Self {
2601 values: out,
2602 dims_per_term: vec![1; term_indices.len()],
2603 }
2604 }
2605
2606 pub fn from_length_scales_aniso(
2624 spec: &TermCollectionSpec,
2625 term_indices: &[usize],
2626 options: &SpatialLengthScaleOptimizationOptions,
2627 ) -> Self {
2628 let mut vals = Vec::new();
2629 let mut dims = Vec::new();
2630 for &term_idx in term_indices {
2631 if let Some(mj) = measure_jet_term_spec(spec, term_idx) {
2635 let seed = measure_jet_psi_seed(mj);
2636 dims.push(seed.len());
2637 vals.extend(seed);
2638 continue;
2639 }
2640 if let Some(cc) = constant_curvature_term_spec(spec, term_idx) {
2646 vals.push(cc.kappa);
2647 dims.push(1);
2648 continue;
2649 }
2650 let length_scale = get_spatial_length_scale(spec, term_idx)
2651 .unwrap_or(options.min_length_scale)
2652 .clamp(options.min_length_scale, options.max_length_scale);
2653 let psi_bar = -length_scale.ln(); if spatial_term_uses_per_axis_psi(spec, term_idx) {
2656 let d = get_spatial_feature_dim(spec, term_idx).unwrap_or(1);
2661 let eta_raw = get_spatial_aniso_log_scales(spec, term_idx)
2662 .expect("predicate guarantees aniso_log_scales is Some");
2663 let eta = center_aniso_log_scales(&eta_raw);
2664 for &eta_a in &eta {
2665 vals.push(psi_bar + eta_a);
2666 }
2667 dims.push(d);
2668 } else {
2669 vals.push(psi_bar);
2676 dims.push(1);
2677 }
2678 }
2679 Self {
2680 values: Array1::from_vec(vals),
2681 dims_per_term: dims,
2682 }
2683 }
2684
2685 pub fn lower_bounds_from_data(
2689 data: ArrayView2<'_, f64>,
2690 spec: &TermCollectionSpec,
2691 term_indices: &[usize],
2692 options: &SpatialLengthScaleOptimizationOptions,
2693 ) -> Self {
2694 let mut values = Array1::<f64>::zeros(term_indices.len());
2695 for (slot, &term_idx) in term_indices.iter().enumerate() {
2696 values[slot] = spatial_term_psi_bounds(data, spec, term_idx, options).0;
2697 }
2698 Self {
2699 values,
2700 dims_per_term: vec![1; term_indices.len()],
2701 }
2702 }
2703
2704 pub fn upper_bounds_from_data(
2706 data: ArrayView2<'_, f64>,
2707 spec: &TermCollectionSpec,
2708 term_indices: &[usize],
2709 options: &SpatialLengthScaleOptimizationOptions,
2710 ) -> Self {
2711 let mut values = Array1::<f64>::zeros(term_indices.len());
2712 for (slot, &term_idx) in term_indices.iter().enumerate() {
2713 values[slot] = spatial_term_psi_bounds(data, spec, term_idx, options).1;
2714 }
2715 Self {
2716 values,
2717 dims_per_term: vec![1; term_indices.len()],
2718 }
2719 }
2720
2721 pub fn lower_bounds_aniso_from_data(
2738 data: ArrayView2<'_, f64>,
2739 spec: &TermCollectionSpec,
2740 term_indices: &[usize],
2741 dims_per_term: &[usize],
2742 options: &SpatialLengthScaleOptimizationOptions,
2743 ) -> Self {
2744 Self::aniso_bounds_from_data(
2745 data,
2746 spec,
2747 term_indices,
2748 dims_per_term,
2749 options,
2750 AnisoBoundEnd::Lower,
2751 )
2752 }
2753
2754 pub fn upper_bounds_aniso_from_data(
2758 data: ArrayView2<'_, f64>,
2759 spec: &TermCollectionSpec,
2760 term_indices: &[usize],
2761 dims_per_term: &[usize],
2762 options: &SpatialLengthScaleOptimizationOptions,
2763 ) -> Self {
2764 Self::aniso_bounds_from_data(
2765 data,
2766 spec,
2767 term_indices,
2768 dims_per_term,
2769 options,
2770 AnisoBoundEnd::Upper,
2771 )
2772 }
2773
2774 fn aniso_bounds_from_data(
2780 data: ArrayView2<'_, f64>,
2781 spec: &TermCollectionSpec,
2782 term_indices: &[usize],
2783 dims_per_term: &[usize],
2784 options: &SpatialLengthScaleOptimizationOptions,
2785 end: AnisoBoundEnd,
2786 ) -> Self {
2787 assert_eq!(term_indices.len(), dims_per_term.len());
2788 let total: usize = dims_per_term.iter().sum();
2789 let mut values = Array1::<f64>::zeros(total);
2790 let mut cursor = 0;
2791 for (slot, &term_idx) in term_indices.iter().enumerate() {
2792 let d = dims_per_term[slot];
2793 if let Some(mj) = measure_jet_term_spec(spec, term_idx) {
2796 let bounds = measure_jet_psi_bound_values(mj, matches!(end, AnisoBoundEnd::Upper));
2797 for (offset, bound) in bounds.into_iter().enumerate() {
2798 if offset < d {
2799 values[cursor + offset] = bound;
2800 }
2801 }
2802 cursor += d;
2803 continue;
2804 }
2805 if constant_curvature_term_spec(spec, term_idx).is_some() {
2808 let (lo, hi) = constant_curvature_kappa_bounds(data, spec, term_idx);
2809 if d >= 1 {
2810 values[cursor] = match end {
2811 AnisoBoundEnd::Lower => lo,
2812 AnisoBoundEnd::Upper => hi,
2813 };
2814 }
2815 cursor += d;
2816 continue;
2817 }
2818 let psi_bound = {
2819 let (lo, hi) = spatial_term_psi_bounds(data, spec, term_idx, options);
2820 match end {
2821 AnisoBoundEnd::Lower => lo,
2822 AnisoBoundEnd::Upper => hi,
2823 }
2824 };
2825 let axis_offsets = if d <= 1 {
2826 vec![0.0; d]
2827 } else {
2828 get_spatial_aniso_log_scales(spec, term_idx)
2829 .filter(|eta| eta.len() == d)
2830 .map(|eta| center_aniso_log_scales(&eta))
2831 .unwrap_or_else(|| vec![0.0; d])
2832 };
2833 for offset in 0..d {
2834 values[cursor + offset] = psi_bound + axis_offsets[offset];
2835 }
2836 cursor += d;
2837 }
2838 Self {
2839 values,
2840 dims_per_term: dims_per_term.to_vec(),
2841 }
2842 }
2843
2844 pub fn reseed_from_data(
2853 mut self,
2854 data: ArrayView2<'_, f64>,
2855 spec: &TermCollectionSpec,
2856 term_indices: &[usize],
2857 options: &SpatialLengthScaleOptimizationOptions,
2858 ) -> Self {
2859 assert_eq!(term_indices.len(), self.dims_per_term.len());
2860 let mut cursor = 0;
2861 for (slot, &term_idx) in term_indices.iter().enumerate() {
2862 let d = self.dims_per_term[slot];
2863 if measure_jet_term_spec(spec, term_idx).is_some() {
2866 cursor += d;
2867 continue;
2868 }
2869 if constant_curvature_term_spec(spec, term_idx).is_some() {
2873 cursor += d;
2874 continue;
2875 }
2876 let Some(psi_bar_new) = spatial_term_psi_seed(data, spec, term_idx, options) else {
2877 cursor += d;
2878 continue;
2879 };
2880 if d == 0 {
2881 continue;
2882 }
2883 let current: Vec<f64> = self.values.slice(s![cursor..cursor + d]).to_vec();
2884 let psi_bar_old = current.iter().sum::<f64>() / d as f64;
2885 for (offset, &old_value) in current.iter().enumerate() {
2886 self.values[cursor + offset] = psi_bar_new + (old_value - psi_bar_old);
2887 }
2888 cursor += d;
2889 }
2890 self
2891 }
2892
2893 pub fn clamp_to_bounds(
2904 mut self,
2905 lower: &SpatialLogKappaCoords,
2906 upper: &SpatialLogKappaCoords,
2907 ) -> Self {
2908 assert_eq!(self.values.len(), lower.values.len());
2909 assert_eq!(self.values.len(), upper.values.len());
2910 let mut n_projected = 0usize;
2911 let mut worst_delta = 0.0_f64;
2912 for idx in 0..self.values.len() {
2913 let lo = lower.values[idx];
2914 let hi = upper.values[idx];
2915 if !(lo.is_finite() && hi.is_finite()) {
2916 continue;
2917 }
2918 let v = self.values[idx];
2919 if v < lo {
2920 worst_delta = worst_delta.max(lo - v);
2921 self.values[idx] = lo;
2922 n_projected += 1;
2923 } else if v > hi {
2924 worst_delta = worst_delta.max(v - hi);
2925 self.values[idx] = hi;
2926 n_projected += 1;
2927 }
2928 }
2929 if n_projected > 0 {
2930 log::info!(
2931 "[spatial-kappa] projected {n_projected}/{} ψ seed coords into data-derived bounds \
2932 (worst excess={worst_delta:.3} log units); user length_scale falls outside \
2933 [{KERNEL_RANGE_MIN_DIAMETER_FRACTION}/r_max, {KERNEL_RANGE_MAX_SPACING_MULTIPLE}/r_min] geometry window",
2934 self.values.len()
2935 );
2936 }
2937 self
2938 }
2939
2940 pub fn from_theta_tail_with_dims(
2942 theta: &Array1<f64>,
2943 start: usize,
2944 dims_per_term: Vec<usize>,
2945 ) -> Self {
2946 let total: usize = dims_per_term.iter().sum();
2947 Self {
2948 values: theta.slice(s![start..start + total]).to_owned(),
2949 dims_per_term,
2950 }
2951 }
2952
2953 pub fn len(&self) -> usize {
2955 self.values.len()
2956 }
2957
2958 pub fn dims_per_term(&self) -> &[usize] {
2960 &self.dims_per_term
2961 }
2962
2963 fn term_offset(&self, term_idx: usize) -> usize {
2965 self.dims_per_term[..term_idx].iter().sum()
2966 }
2967
2968 pub fn term_slice(&self, term_idx: usize) -> &[f64] {
2970 let offset = self.term_offset(term_idx);
2971 let d = self.dims_per_term[term_idx];
2972 &self.values.as_slice().unwrap()[offset..offset + d]
2973 }
2974
2975 pub fn as_array(&self) -> &Array1<f64> {
2976 &self.values
2977 }
2978
2979 pub fn set_scalar_slot(&mut self, slot: usize, value: f64) -> bool {
2985 if slot >= self.dims_per_term.len() || self.dims_per_term[slot] != 1 {
2986 return false;
2987 }
2988 let offset = self.term_offset(slot);
2989 self.values[offset] = value;
2990 true
2991 }
2992
2993 pub fn split_at(&self, mid: usize) -> (Self, Self) {
2996 let flat_mid: usize = self.dims_per_term[..mid].iter().sum();
2997 (
2998 Self {
2999 values: self.values.slice(s![0..flat_mid]).to_owned(),
3000 dims_per_term: self.dims_per_term[..mid].to_vec(),
3001 },
3002 Self {
3003 values: self.values.slice(s![flat_mid..]).to_owned(),
3004 dims_per_term: self.dims_per_term[mid..].to_vec(),
3005 },
3006 )
3007 }
3008
3009 pub fn apply_tospec(
3016 &self,
3017 spec: &TermCollectionSpec,
3018 term_indices: &[usize],
3019 ) -> Result<TermCollectionSpec, EstimationError> {
3020 if term_indices.len() != self.dims_per_term.len() {
3021 crate::bail_invalid_estim!(
3022 "SpatialLogKappaCoords::apply_tospec: term count mismatch: \
3023 term_indices={} dims_per_term={}",
3024 term_indices.len(),
3025 self.dims_per_term.len()
3026 );
3027 }
3028 let mut updated = spec.clone();
3029 for (slot, &term_idx) in term_indices.iter().enumerate() {
3030 let psi = self.term_slice(slot);
3031 let d = self.dims_per_term[slot];
3032 if measure_jet_term_spec(&updated, term_idx).is_some() {
3035 set_measure_jet_psi_dials(&mut updated, term_idx, psi)?;
3036 continue;
3037 }
3038 if constant_curvature_term_spec(&updated, term_idx).is_some() {
3042 set_constant_curvature_kappa(&mut updated, term_idx, psi)?;
3043 continue;
3044 }
3045 let (next_length_scale, next_aniso) = spatial_term_psi_to_length_scale_and_aniso(psi);
3046 if (d == 1 || next_length_scale.is_some())
3047 && let Some(length_scale) = next_length_scale
3048 {
3049 set_spatial_length_scale(&mut updated, term_idx, length_scale)?;
3050 }
3051 if let Some(eta) = next_aniso {
3052 set_spatial_aniso_log_scales(&mut updated, term_idx, eta)?;
3053 }
3054 }
3055 Ok(updated)
3056 }
3057}
3058
3059pub fn center_aniso_log_scales(eta: &[f64]) -> Vec<f64> {
3060 if eta.len() <= 1 {
3061 return eta.to_vec();
3062 }
3063 let mean = eta.iter().sum::<f64>() / eta.len() as f64;
3064 eta.iter()
3065 .map(|&v| {
3066 let centered = v - mean;
3067 if centered.abs() <= 1e-15 {
3068 0.0
3069 } else {
3070 centered
3071 }
3072 })
3073 .collect()
3074}
3075
3076pub fn spatial_term_uses_per_axis_psi(resolvedspec: &TermCollectionSpec, term_idx: usize) -> bool {
3079 if let Some(mj) = measure_jet_term_spec(resolvedspec, term_idx) {
3080 return measure_jet_enrolls_psi(mj);
3081 }
3082 let Some(d) = get_spatial_feature_dim(resolvedspec, term_idx) else {
3083 return false;
3084 };
3085 if d <= 1 {
3086 return false;
3087 }
3088 let Some(eta) = get_spatial_aniso_log_scales(resolvedspec, term_idx) else {
3089 return false;
3090 };
3091 if eta.len() != d {
3092 return false;
3093 }
3094 !matches!(
3095 resolvedspec.smooth_terms.get(term_idx).map(|term| &term.basis),
3096 Some(SmoothBasisSpec::Duchon { .. })
3097 )
3098}
3099
3100pub fn set_spatial_length_scale(
3101 spec: &mut TermCollectionSpec,
3102 term_idx: usize,
3103 length_scale: f64,
3104) -> Result<(), EstimationError> {
3105 let Some(term) = spec.smooth_terms.get_mut(term_idx) else {
3106 crate::bail_invalid_estim!("spatial length-scale term index {term_idx} out of range");
3107 };
3108 match &mut term.basis {
3109 SmoothBasisSpec::ThinPlate { spec, .. } => {
3110 spec.length_scale = length_scale;
3111 Ok(())
3112 }
3113 SmoothBasisSpec::Matern { spec, .. } => {
3114 spec.length_scale = length_scale;
3115 Ok(())
3116 }
3117 SmoothBasisSpec::Duchon { spec, .. } => {
3118 spec.length_scale = Some(length_scale);
3119 Ok(())
3120 }
3121 _ => Err(EstimationError::InvalidInput(format!(
3122 "term '{}' does not expose a spatial length scale",
3123 term.name
3124 ))),
3125 }
3126}
3127
3128pub fn get_spatial_length_scale(spec: &TermCollectionSpec, term_idx: usize) -> Option<f64> {
3129 spec.smooth_terms
3130 .get(term_idx)
3131 .and_then(|term| match &term.basis {
3132 SmoothBasisSpec::ThinPlate { spec, .. } => Some(spec.length_scale),
3133 SmoothBasisSpec::Matern { spec, .. } => Some(spec.length_scale),
3134 SmoothBasisSpec::Duchon { spec, .. } => spec.length_scale,
3135 _ => None,
3136 })
3137}
3138
3139pub fn spatial_term_supports_hyper_optimization(spec: &TermCollectionSpec, term_idx: usize) -> bool {
3140 if let Some(term) = spec.smooth_terms.get(term_idx)
3146 && let SmoothBasisSpec::ThinPlate { .. } = &term.basis
3147 {
3148 return false;
3149 }
3150
3151 if let Some(term) = spec.smooth_terms.get(term_idx)
3176 && let SmoothBasisSpec::Matern { .. } = &term.basis
3177 {
3178 return true;
3179 }
3180
3181 if let Some(mj) = measure_jet_term_spec(spec, term_idx) {
3184 return measure_jet_enrolls_psi(mj);
3185 }
3186
3187 if constant_curvature_term_spec(spec, term_idx).is_some() {
3194 return true;
3195 }
3196
3197 get_spatial_length_scale(spec, term_idx).is_some()
3198}
3199
3200pub fn measure_jet_term_spec(
3203 spec: &TermCollectionSpec,
3204 term_idx: usize,
3205) -> Option<&crate::basis::MeasureJetBasisSpec> {
3206 spec.smooth_terms
3207 .get(term_idx)
3208 .and_then(|term| match &term.basis {
3209 SmoothBasisSpec::MeasureJet { spec, .. } => Some(spec),
3210 _ => None,
3211 })
3212}
3213
3214pub fn measure_jet_enrolls_psi(mj: &crate::basis::MeasureJetBasisSpec) -> bool {
3221 measure_jet_learns_length_scale(mj)
3230 || (mj.tau0 > 0.0 && crate::basis::measure_jet_multiscale_mode(mj))
3231}
3232
3233pub fn measure_jet_learns_length_scale(mj: &crate::basis::MeasureJetBasisSpec) -> bool {
3236 mj.learn_length_scale
3237}
3238
3239pub fn freeze_measure_jet_length_scale_learning(spec: &mut TermCollectionSpec) -> usize {
3240 let mut frozen = 0;
3241 for term in spec.smooth_terms.iter_mut() {
3242 if let SmoothBasisSpec::MeasureJet { spec: mj, .. } = &mut term.basis
3243 && mj.learn_length_scale
3244 {
3245 mj.learn_length_scale = false;
3246 frozen += 1;
3247 }
3248 }
3249 frozen
3250}
3251
3252pub const MEASURE_JET_PSI_ALPHA_BOUNDS: (f64, f64) = (-1.0, 3.0);
3260
3261pub const MEASURE_JET_PSI_LN_TAU_BOUNDS: (f64, f64) = (-18.420680743952367, 4.605170185988092);
3262
3263pub const MEASURE_JET_PSI_LN_LENGTH_SCALE_BOUNDS: (f64, f64) = (-6.907755278982137, 4.605170185988092);
3269
3270pub fn measure_jet_penalty_psi_dim(mj: &crate::basis::MeasureJetBasisSpec) -> usize {
3278 if crate::basis::measure_jet_multiscale_mode(mj) {
3279 2
3280 } else {
3281 0
3282 }
3283}
3284
3285pub fn measure_jet_psi_dim(mj: &crate::basis::MeasureJetBasisSpec) -> usize {
3289 usize::from(measure_jet_learns_length_scale(mj)) + measure_jet_penalty_psi_dim(mj)
3290}
3291
3292pub fn measure_jet_psi_seed(mj: &crate::basis::MeasureJetBasisSpec) -> Vec<f64> {
3297 let mut seed = Vec::with_capacity(measure_jet_psi_dim(mj));
3298 if measure_jet_learns_length_scale(mj) {
3299 let ell = if mj.length_scale > 0.0 {
3303 mj.length_scale
3304 } else {
3305 1.0
3306 };
3307 seed.push(ell.ln());
3308 }
3309 if measure_jet_penalty_psi_dim(mj) > 0 {
3310 let ln_tau = mj.tau0.max(f64::MIN_POSITIVE).ln();
3312 seed.extend_from_slice(&[mj.alpha, ln_tau]);
3313 }
3314 seed
3315}
3316
3317pub fn measure_jet_psi_bound_values(mj: &crate::basis::MeasureJetBasisSpec, upper: bool) -> Vec<f64> {
3320 let pick = |b: (f64, f64)| if upper { b.1 } else { b.0 };
3321 let mut bounds = Vec::with_capacity(measure_jet_psi_dim(mj));
3322 if measure_jet_learns_length_scale(mj) {
3323 bounds.push(pick(MEASURE_JET_PSI_LN_LENGTH_SCALE_BOUNDS));
3324 }
3325 if measure_jet_penalty_psi_dim(mj) > 0 {
3326 bounds.push(pick(MEASURE_JET_PSI_ALPHA_BOUNDS));
3328 bounds.push(pick(MEASURE_JET_PSI_LN_TAU_BOUNDS));
3329 }
3330 bounds
3331}
3332
3333pub fn apply_measure_jet_psi(
3338 mj: &mut crate::basis::MeasureJetBasisSpec,
3339 psi: &[f64],
3340) -> Result<bool, EstimationError> {
3341 if psi.len() != measure_jet_psi_dim(mj) {
3342 crate::bail_invalid_estim!(
3343 "measure-jet ψ write-back dimension mismatch: got {} values for a {}-dial term",
3344 psi.len(),
3345 measure_jet_psi_dim(mj)
3346 );
3347 }
3348 let mut changed = false;
3349 let mut cursor = 0usize;
3353 if measure_jet_learns_length_scale(mj) {
3354 let next_ell = psi[cursor].exp();
3355 cursor += 1;
3356 if !(next_ell.is_finite() && next_ell > 0.0) {
3357 crate::bail_invalid_estim!(
3358 "measure-jet ψ write-back produced a non-finite/non-positive length_scale (ℓ={next_ell})"
3359 );
3360 }
3361 if next_ell != mj.length_scale {
3362 mj.length_scale = next_ell;
3363 changed = true;
3364 }
3365 }
3366 if measure_jet_penalty_psi_dim(mj) > 0 {
3367 let next_alpha = psi[cursor];
3370 let next_tau = psi[cursor + 1].exp();
3371 if !(next_alpha.is_finite() && next_tau.is_finite() && next_tau > 0.0) {
3372 crate::bail_invalid_estim!(
3373 "measure-jet ψ write-back produced non-finite dials (alpha={next_alpha}, tau={next_tau})"
3374 );
3375 }
3376 if next_alpha != mj.alpha {
3377 mj.alpha = next_alpha;
3378 changed = true;
3379 }
3380 if next_tau != mj.tau0 {
3381 mj.tau0 = next_tau;
3382 changed = true;
3383 }
3384 }
3385 Ok(changed)
3386}
3387
3388pub fn set_measure_jet_psi_dials(
3391 spec: &mut TermCollectionSpec,
3392 term_idx: usize,
3393 psi: &[f64],
3394) -> Result<bool, EstimationError> {
3395 let Some(term) = spec.smooth_terms.get_mut(term_idx) else {
3396 crate::bail_invalid_estim!("measure-jet ψ write-back: term index {term_idx} out of range");
3397 };
3398 set_single_term_measure_jet_psi_dials(term, psi)
3399}
3400
3401pub fn set_single_term_measure_jet_psi_dials(
3406 term: &mut SmoothTermSpec,
3407 psi: &[f64],
3408) -> Result<bool, EstimationError> {
3409 let SmoothBasisSpec::MeasureJet { spec: mj, .. } = &mut term.basis else {
3410 crate::bail_invalid_estim!("measure-jet ψ write-back targeted a non-measure-jet term");
3411 };
3412 apply_measure_jet_psi(mj, psi)
3413}
3414
3415pub fn constant_curvature_term_spec(
3418 spec: &TermCollectionSpec,
3419 term_idx: usize,
3420) -> Option<&crate::basis::ConstantCurvatureBasisSpec> {
3421 spec.smooth_terms
3422 .get(term_idx)
3423 .and_then(|term| match &term.basis {
3424 SmoothBasisSpec::ConstantCurvature { spec, .. } => Some(spec),
3425 _ => None,
3426 })
3427}
3428
3429pub const CONSTANT_CURVATURE_KAPPA_CHART_FRACTION: f64 = 0.5;
3437
3438pub const CONSTANT_CURVATURE_MIN_CHART_RADIUS2: f64 = 1e-8;
3442
3443pub fn constant_curvature_kappa_bounds(
3448 data: ArrayView2<'_, f64>,
3449 spec: &TermCollectionSpec,
3450 term_idx: usize,
3451) -> (f64, f64) {
3452 let feature_cols = match spec.smooth_terms.get(term_idx).map(|t| &t.basis) {
3453 Some(SmoothBasisSpec::ConstantCurvature { feature_cols, .. }) => feature_cols,
3454 _ => return (-1.0, 1.0),
3455 };
3456 let mut max_r2 = CONSTANT_CURVATURE_MIN_CHART_RADIUS2;
3457 for row in data.outer_iter() {
3458 let mut r2 = 0.0_f64;
3459 for &c in feature_cols.iter() {
3460 if let Some(&v) = row.get(c)
3461 && v.is_finite()
3462 {
3463 r2 += v * v;
3464 }
3465 }
3466 if r2 > max_r2 {
3467 max_r2 = r2;
3468 }
3469 }
3470 let half = CONSTANT_CURVATURE_KAPPA_CHART_FRACTION / max_r2;
3471 (-half, half)
3472}
3473
3474pub fn set_constant_curvature_kappa(
3478 spec: &mut TermCollectionSpec,
3479 term_idx: usize,
3480 psi: &[f64],
3481) -> Result<bool, EstimationError> {
3482 let Some(term) = spec.smooth_terms.get_mut(term_idx) else {
3483 crate::bail_invalid_estim!(
3484 "constant-curvature κ write-back: term index {term_idx} out of range"
3485 );
3486 };
3487 set_single_term_constant_curvature_kappa(term, psi)
3488}
3489
3490pub fn set_single_term_constant_curvature_kappa(
3495 term: &mut SmoothTermSpec,
3496 psi: &[f64],
3497) -> Result<bool, EstimationError> {
3498 if psi.len() != 1 {
3499 crate::bail_invalid_estim!(
3500 "constant-curvature κ write-back expects exactly one value, got {}",
3501 psi.len()
3502 );
3503 }
3504 let next_kappa = psi[0];
3505 if !next_kappa.is_finite() {
3506 crate::bail_invalid_estim!(
3507 "constant-curvature κ write-back produced a non-finite κ = {next_kappa}"
3508 );
3509 }
3510 let SmoothBasisSpec::ConstantCurvature { spec: cc, .. } = &mut term.basis else {
3511 crate::bail_invalid_estim!(
3512 "constant-curvature κ write-back targeted a non-constant-curvature term"
3513 );
3514 };
3515 if cc.kappa != next_kappa {
3516 cc.kappa = next_kappa;
3517 Ok(true)
3518 } else {
3519 Ok(false)
3520 }
3521}
3522
3523pub fn spatial_term_has_locked_kappa(spec: &TermCollectionSpec, term_idx: usize) -> bool {
3534 get_spatial_length_scale(spec, term_idx).is_some()
3535 && !spatial_term_uses_per_axis_psi(spec, term_idx)
3536}
3537
3538pub fn all_spatial_terms_kappa_fixed(spec: &TermCollectionSpec) -> bool {
3539 spec.smooth_terms.iter().enumerate().all(|(idx, _)| {
3540 !spatial_term_supports_hyper_optimization(spec, idx)
3541 || spatial_term_has_locked_kappa(spec, idx)
3542 })
3543}
3544
3545pub fn spatial_identifiability_policy(termspec: &SmoothTermSpec) -> Option<&SpatialIdentifiability> {
3546 match &termspec.basis {
3547 SmoothBasisSpec::ThinPlate { spec, .. } => Some(&spec.identifiability),
3548 SmoothBasisSpec::Duchon { spec, .. } => Some(&spec.identifiability),
3549 _ => None,
3550 }
3551}
3552
3553pub const NULLSPACE_WELLDET_DEGENERACY_RHO_SD: f64 = 15.0;
3557
3558pub fn is_nullspace_degeneracy_prior(prior: &gam_spec::RhoPrior) -> bool {
3561 matches!(
3562 prior,
3563 gam_spec::RhoPrior::Normal { mean, sd }
3564 if *mean == 0.0 && *sd == NULLSPACE_WELLDET_DEGENERACY_RHO_SD
3565 )
3566}
3567
3568pub const KERNEL_RANGE_MIN_DIAMETER_FRACTION: f64 = 2.0;
3580
3581pub const KERNEL_RANGE_MAX_SPACING_MULTIPLE: f64 = 1e2;
3586
3587
3588pub fn spatial_term_psi_bounds(
3597 data: ArrayView2<'_, f64>,
3598 spec: &TermCollectionSpec,
3599 term_idx: usize,
3600 options: &SpatialLengthScaleOptimizationOptions,
3601) -> (f64, f64) {
3602 let fallback = (
3603 -options.max_length_scale.ln(),
3604 -options.min_length_scale.ln(),
3605 );
3606 if constant_curvature_term_spec(spec, term_idx).is_some() {
3611 return constant_curvature_kappa_bounds(data, spec, term_idx);
3612 }
3613 let Some(term) = spec.smooth_terms.get(term_idx) else {
3614 return fallback;
3615 };
3616 let aniso = get_spatial_aniso_log_scales(spec, term_idx);
3629 let r_bounds = match spatial_term_center_strategy(term) {
3630 Some(CenterStrategy::UserProvided(centers)) if centers.nrows() >= 2 => {
3631 match aniso.as_deref() {
3632 Some(eta) if eta.len() == centers.ncols() => {
3633 let y = points_in_aniso_y_space(centers.view(), eta);
3634 pairwise_distance_bounds(y.view())
3635 }
3636 _ => pairwise_distance_bounds(centers.view()),
3637 }
3638 }
3639 _ => standardized_spatial_term_data(data, term)
3640 .ok()
3641 .and_then(|x| match aniso.as_deref() {
3642 Some(eta) if eta.len() == x.ncols() => {
3643 let y = points_in_aniso_y_space(x.view(), eta);
3644 pairwise_distance_bounds_sampled(y.view())
3645 }
3646 _ => pairwise_distance_bounds_sampled(x.view()),
3647 }),
3648 };
3649 let Some((r_min, r_max)) = r_bounds else {
3650 return fallback;
3651 };
3652 let psi_lo_data = (KERNEL_RANGE_MIN_DIAMETER_FRACTION / r_max).ln();
3658 let psi_hi_data = (KERNEL_RANGE_MAX_SPACING_MULTIPLE / r_min).ln();
3659 let psi_lo = psi_lo_data.max(fallback.0);
3669 let psi_hi = psi_hi_data.min(fallback.1);
3670 if psi_lo >= psi_hi {
3671 return fallback;
3674 }
3675 (psi_lo, psi_hi)
3676}
3677
3678pub fn spatial_term_psi_seed(
3682 data: ArrayView2<'_, f64>,
3683 spec: &TermCollectionSpec,
3684 term_idx: usize,
3685 options: &SpatialLengthScaleOptimizationOptions,
3686) -> Option<f64> {
3687 if get_spatial_length_scale(spec, term_idx).is_some() {
3688 return None; }
3690 let (psi_lo, psi_hi) = spatial_term_psi_bounds(data, spec, term_idx, options);
3691 Some(0.5 * (psi_lo + psi_hi))
3692}
3693
3694pub fn spatial_term_psi_to_length_scale_and_aniso(psi: &[f64]) -> (Option<f64>, Option<Vec<f64>>) {
3695 if psi.len() <= 1 {
3696 (Some((-psi.first().copied().unwrap_or(0.0)).exp()), None)
3697 } else {
3698 let psi_bar = psi.iter().sum::<f64>() / psi.len() as f64;
3699 (
3700 Some((-psi_bar).exp()),
3701 Some(psi.iter().map(|&value| value - psi_bar).collect()),
3702 )
3703 }
3704}
3705
3706pub fn get_spatial_aniso_log_scales(
3708 spec: &TermCollectionSpec,
3709 term_idx: usize,
3710) -> Option<Vec<f64>> {
3711 spec.smooth_terms
3712 .get(term_idx)
3713 .and_then(|term| match &term.basis {
3714 SmoothBasisSpec::Matern { spec, .. } => spec.aniso_log_scales.clone(),
3715 SmoothBasisSpec::Duchon { spec, .. } => spec.aniso_log_scales.clone(),
3716 _ => None,
3717 })
3718}
3719
3720pub fn response_aware_axis_contrasts(
3740 x: ndarray::ArrayView2<'_, f64>,
3741 y: ndarray::ArrayView1<'_, f64>,
3742) -> Option<Vec<f64>> {
3743 let n = x.nrows();
3744 let d = x.ncols();
3745 if d <= 1 || n < 4 || y.len() != n {
3746 return None;
3747 }
3748 if x.iter().any(|v| !v.is_finite()) || y.iter().any(|v| !v.is_finite()) {
3749 return None;
3750 }
3751 let mut scores = Vec::with_capacity(d);
3752 for a in 0..d {
3753 let mut order: Vec<usize> = (0..n).collect();
3754 let col = x.column(a);
3755 order.sort_by(|&i, &j| {
3756 col[i]
3757 .partial_cmp(&col[j])
3758 .unwrap_or(std::cmp::Ordering::Equal)
3759 });
3760 let mut tv = 0.0_f64;
3761 for w in order.windows(2) {
3762 let diff = y[w[1]] - y[w[0]];
3763 tv += diff * diff;
3764 }
3765 scores.push(-0.5 * (tv + 1e-12).ln());
3767 }
3768 if scores.iter().any(|v| !v.is_finite()) {
3769 return None;
3770 }
3771 let mean = scores.iter().sum::<f64>() / d as f64;
3772 let centered: Vec<f64> = scores.iter().map(|&s| s - mean).collect();
3773 if centered.iter().all(|&v| v.abs() < 1e-9) {
3776 return None;
3777 }
3778 Some(centered)
3779}
3780
3781pub fn apply_response_aware_anisotropy_seed(
3790 data: ArrayView2<'_, f64>,
3791 y: ndarray::ArrayView1<'_, f64>,
3792 spec: &mut TermCollectionSpec,
3793 spatial_terms: &[usize],
3794) {
3795 const MAX_NUDGE: f64 = std::f64::consts::LN_2;
3800 for &term_idx in spatial_terms {
3801 let Some(current_eta) = get_spatial_aniso_log_scales(spec, term_idx) else {
3802 continue;
3803 };
3804 let d = current_eta.len();
3805 if d <= 1 {
3806 continue;
3807 }
3808 let Some(term) = spec.smooth_terms.get(term_idx) else {
3809 continue;
3810 };
3811 let feature_cols = term.basis.structural_feature_cols();
3812 if feature_cols.len() != d {
3813 continue;
3814 }
3815 let Ok(x) = select_columns(data, &feature_cols) else {
3816 continue;
3817 };
3818 let Some(contrast) = response_aware_axis_contrasts(x.view(), y) else {
3819 continue;
3820 };
3821 let nudged: Vec<f64> = current_eta
3822 .iter()
3823 .zip(contrast.iter())
3824 .map(|(&eta_a, &c_a)| eta_a + c_a.clamp(-MAX_NUDGE, MAX_NUDGE))
3825 .collect();
3826 if let Err(err) = set_spatial_aniso_log_scales(spec, term_idx, nudged) {
3829 log::debug!(
3830 "[spatial-kappa] response-aware anisotropy seed skipped for term {term_idx}: {err}"
3831 );
3832 }
3833 }
3834}
3835
3836pub fn get_spatial_feature_dim(spec: &TermCollectionSpec, term_idx: usize) -> Option<usize> {
3838 spec.smooth_terms
3839 .get(term_idx)
3840 .and_then(|term| match &term.basis {
3841 SmoothBasisSpec::ThinPlate { feature_cols, .. } => Some(feature_cols.len()),
3842 SmoothBasisSpec::Matern { feature_cols, .. } => Some(feature_cols.len()),
3843 SmoothBasisSpec::Duchon { feature_cols, .. } => Some(feature_cols.len()),
3844 _ => None,
3845 })
3846}
3847
3848pub fn log_spatial_aniso_scales(spec: &TermCollectionSpec) {
3855 for (term_idx, term) in spec.smooth_terms.iter().enumerate() {
3856 let (aniso, length_scale) = match &term.basis {
3857 SmoothBasisSpec::Matern { spec, .. } => {
3858 (spec.aniso_log_scales.as_ref(), Some(spec.length_scale))
3859 }
3860 SmoothBasisSpec::Duchon { spec, .. } => {
3861 (spec.aniso_log_scales.as_ref(), spec.length_scale)
3862 }
3863 _ => (None, None),
3864 };
3865 let Some(eta) = aniso else { continue };
3866 if eta.is_empty() {
3867 continue;
3868 }
3869 let mut lines = match length_scale {
3870 Some(ls) => format!(
3871 "[spatial-kappa] term {} (\"{}\"): anisotropic length scales optimized (global length_scale={:.4})",
3872 term_idx, term.name, ls
3873 ),
3874 None => format!(
3875 "[spatial-kappa] term {} (\"{}\"): pure Duchon shape anisotropy optimized",
3876 term_idx, term.name
3877 ),
3878 };
3879 for (a, &eta_a) in eta.iter().enumerate() {
3880 if let Some(ls) = length_scale {
3881 let length_a = ls * (-eta_a).exp();
3882 let kappa_a = (1.0 / ls) * eta_a.exp();
3883 lines.push_str(&format!(
3884 "\n axis {}: eta={:+.4}, length={:.4}, kappa={:.4}",
3885 a, eta_a, length_a, kappa_a
3886 ));
3887 } else {
3888 lines.push_str(&format!("\n axis {}: eta={:+.4}", a, eta_a));
3889 }
3890 }
3891 log::info!("{}", lines);
3892 }
3893}
3894
3895pub fn set_spatial_aniso_log_scales(
3897 spec: &mut TermCollectionSpec,
3898 term_idx: usize,
3899 eta: Vec<f64>,
3900) -> Result<(), EstimationError> {
3901 let eta = center_aniso_log_scales(&eta);
3902 let Some(term) = spec.smooth_terms.get_mut(term_idx) else {
3903 crate::bail_invalid_estim!("spatial aniso_log_scales term index {term_idx} out of range");
3904 };
3905 match &mut term.basis {
3906 SmoothBasisSpec::Matern { spec, .. } => {
3907 spec.aniso_log_scales = Some(eta);
3908 Ok(())
3909 }
3910 SmoothBasisSpec::Duchon { spec, .. } => {
3911 spec.aniso_log_scales = Some(eta);
3912 Ok(())
3913 }
3914 _ => Err(EstimationError::InvalidInput(format!(
3915 "term '{}' does not support aniso_log_scales",
3916 term.name
3917 ))),
3918 }
3919}
3920
3921pub fn sync_aniso_contrasts_from_metadata(
3928 spec: &mut TermCollectionSpec,
3929 design: &SmoothDesign,
3930) {
3931 for (term_idx, term) in design.terms.iter().enumerate() {
3932 let meta_aniso = match &term.metadata {
3933 BasisMetadata::Matern {
3934 aniso_log_scales, ..
3935 } => aniso_log_scales.clone(),
3936 BasisMetadata::Duchon {
3937 aniso_log_scales, ..
3938 } => aniso_log_scales.clone(),
3939 _ => None,
3940 };
3941 if let Some(eta) = meta_aniso
3942 && eta.len() > 1
3943 {
3944 set_spatial_aniso_log_scales(spec, term_idx, eta).ok();
3945 }
3946 }
3947}
3948
3949#[derive(Debug, Clone)]
3950pub struct SpatialLengthScaleOptimizationOptions {
3951 pub enabled: bool,
3955 pub max_outer_iter: usize,
3957 pub rel_tol: f64,
3959 pub log_step: f64,
3961 pub min_length_scale: f64,
3963 pub max_length_scale: f64,
3965 pub pilot_subsample_threshold: usize,
3978}
3979
3980impl Default for SpatialLengthScaleOptimizationOptions {
3981 fn default() -> Self {
3982 Self {
3983 enabled: true,
3984 max_outer_iter: 80,
3985 rel_tol: 1e-4,
3986 log_step: std::f64::consts::LN_2,
3987 min_length_scale: 1e-3,
3988 max_length_scale: 1e3,
3989 pilot_subsample_threshold: 10_000,
3990 }
3991 }
3992}
3993
3994impl SpatialLengthScaleOptimizationOptions {
3995 pub fn validate(&self) -> Result<(), String> {
4013 if !self.min_length_scale.is_finite() || self.min_length_scale <= 0.0 {
4014 return Err(SmoothError::invalid_config(format!(
4015 "SpatialLengthScaleOptimizationOptions::min_length_scale must be > 0 and finite, got {}",
4016 self.min_length_scale
4017 ))
4018 .into());
4019 }
4020 if !self.max_length_scale.is_finite() || self.max_length_scale <= 0.0 {
4021 return Err(SmoothError::invalid_config(format!(
4022 "SpatialLengthScaleOptimizationOptions::max_length_scale must be > 0 and finite, got {}",
4023 self.max_length_scale
4024 ))
4025 .into());
4026 }
4027 if self.min_length_scale >= self.max_length_scale {
4028 return Err(SmoothError::invalid_config(format!(
4029 "SpatialLengthScaleOptimizationOptions requires min_length_scale < max_length_scale, got min={} max={}",
4030 self.min_length_scale, self.max_length_scale
4031 ))
4032 .into());
4033 }
4034 if !self.rel_tol.is_finite() || self.rel_tol <= 0.0 {
4035 return Err(SmoothError::invalid_config(format!(
4036 "SpatialLengthScaleOptimizationOptions::rel_tol must be > 0 and finite, got {}",
4037 self.rel_tol
4038 ))
4039 .into());
4040 }
4041 if !self.log_step.is_finite() || self.log_step <= 0.0 {
4042 return Err(SmoothError::invalid_config(format!(
4043 "SpatialLengthScaleOptimizationOptions::log_step must be > 0 and finite, got {}",
4044 self.log_step
4045 ))
4046 .into());
4047 }
4048 Ok(())
4049 }
4050}
4051
4052#[derive(Debug, Clone)]
4053pub struct RandomEffectBlock {
4054 pub name: String,
4055 pub group_ids: Vec<Option<usize>>,
4058 pub num_groups: usize,
4059 pub kept_levels: Vec<u64>,
4060}
4061
4062pub const BLOCK_SPARSE_ZERO_EPS: f64 = 1e-12;
4063
4064pub const BLOCK_SPARSE_MAX_DENSITY: f64 = 0.20;
4065
4066pub fn blocks_have_intrinsic_sparse_structure(blocks: &[DesignBlock]) -> bool {
4067 blocks
4068 .iter()
4069 .any(|block| matches!(block, DesignBlock::Sparse(_) | DesignBlock::RandomEffect(_)))
4070}
4071
4072pub fn sparse_compatible_block_nnz(block: &DesignBlock) -> Option<usize> {
4073 match block {
4074 DesignBlock::Intercept(n) => Some(*n),
4075 DesignBlock::RandomEffect(op) => {
4076 Some(op.group_ids.iter().filter(|gid| gid.is_some()).count())
4077 }
4078 DesignBlock::Sparse(sparse) => Some(sparse.val().len()),
4079 DesignBlock::Dense(dense) => dense.as_dense_ref().map(|matrix| {
4080 matrix
4081 .iter()
4082 .filter(|&&value| value.abs() > BLOCK_SPARSE_ZERO_EPS)
4083 .count()
4084 }),
4085 }
4086}
4087
4088pub fn try_build_sparse_design_from_blocks(
4089 blocks: &[DesignBlock],
4090) -> Result<Option<DesignMatrix>, BasisError> {
4091 if blocks.is_empty() {
4092 return Ok(None);
4093 }
4094 let nrows = blocks[0].nrows();
4095 let ncols: usize = blocks.iter().map(DesignBlock::ncols).sum();
4096 if nrows == 0 || ncols == 0 || ncols <= 32 {
4097 return Ok(None);
4098 }
4099
4100 let preserve_sparse_storage = blocks_have_intrinsic_sparse_structure(blocks);
4101 let sparse_nnz_limit = if preserve_sparse_storage {
4102 usize::MAX
4103 } else {
4104 let total_cells = nrows.saturating_mul(ncols);
4105 ((total_cells as f64) * BLOCK_SPARSE_MAX_DENSITY).floor() as usize
4106 };
4107 let mut nnz = 0usize;
4108 for block in blocks {
4109 let block_nnz = if let Some(block_nnz) = sparse_compatible_block_nnz(block) {
4110 block_nnz
4111 } else {
4112 return Ok(None);
4113 };
4114 nnz = nnz.saturating_add(block_nnz);
4115 if nnz > sparse_nnz_limit {
4116 return Ok(None);
4117 }
4118 }
4119
4120 let mut triplets = Vec::<Triplet<usize, usize, f64>>::with_capacity(nnz);
4121 let mut col_offset = 0usize;
4122 for block in blocks {
4123 match block {
4124 DesignBlock::Intercept(n) => {
4125 for row in 0..*n {
4126 triplets.push(Triplet::new(row, col_offset, 1.0));
4127 }
4128 }
4129 DesignBlock::RandomEffect(op) => {
4130 for (row, group_id) in op.group_ids.iter().enumerate() {
4131 if let Some(group) = group_id {
4132 triplets.push(Triplet::new(row, col_offset + group, 1.0));
4133 }
4134 }
4135 }
4136 DesignBlock::Sparse(sparse) => {
4137 let (symbolic, values) = sparse.parts();
4138 let col_ptr = symbolic.col_ptr();
4139 let row_idx = symbolic.row_idx();
4140 for col in 0..sparse.ncols() {
4141 for idx in col_ptr[col]..col_ptr[col + 1] {
4142 let value = values[idx];
4143 if value.abs() > BLOCK_SPARSE_ZERO_EPS {
4144 triplets.push(Triplet::new(row_idx[idx], col_offset + col, value));
4145 }
4146 }
4147 }
4148 }
4149 DesignBlock::Dense(dense) => {
4150 let matrix = dense.as_dense_ref().ok_or_else(|| {
4151 BasisError::InvalidInput(
4152 "sparse-compatible block assembly requires materialized dense blocks"
4153 .to_string(),
4154 )
4155 })?;
4156 for row in 0..matrix.nrows() {
4157 for col in 0..matrix.ncols() {
4158 let value = matrix[[row, col]];
4159 if value.abs() > BLOCK_SPARSE_ZERO_EPS {
4160 triplets.push(Triplet::new(row, col_offset + col, value));
4161 }
4162 }
4163 }
4164 }
4165 }
4166 col_offset += block.ncols();
4167 }
4168
4169 let sparse = SparseColMat::try_new_from_triplets(nrows, ncols, &triplets).map_err(|_| {
4170 BasisError::SparseCreation("failed to assemble sparse term-collection design".to_string())
4171 })?;
4172 Ok(Some(DesignMatrix::Sparse(
4173 gam_linalg::matrix::SparseDesignMatrix::new(sparse),
4174 )))
4175}
4176
4177pub fn assemble_term_collection_design_matrix(
4178 blocks: Vec<DesignBlock>,
4179) -> Result<DesignMatrix, BasisError> {
4180 if let Some(sparse) = try_build_sparse_design_from_blocks(&blocks)? {
4181 return Ok(sparse);
4182 }
4183 let block_op = BlockDesignOperator::new(blocks).map_err(|e| {
4184 BasisError::InvalidInput(format!("failed to build block design operator: {e}"))
4185 })?;
4186 Ok(DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(
4187 Arc::new(block_op),
4188 )))
4189}
4190
4191pub fn select_columns(data: ArrayView2<'_, f64>, cols: &[usize]) -> Result<Array2<f64>, BasisError> {
4192 let n = data.nrows();
4193 let p = data.ncols();
4194 for &c in cols {
4195 if c >= p {
4196 crate::bail_dim_basis!("feature column {c} is out of bounds for data with {p} columns");
4197 }
4198 }
4199 let mut out = Array2::<f64>::zeros((n, cols.len()));
4200 for (j, &c) in cols.iter().enumerate() {
4201 out.column_mut(j).assign(&data.column(c));
4202 }
4203 Ok(out)
4204}
4205
4206pub fn nonfinite_value_label(value: f64) -> &'static str {
4207 if value.is_nan() {
4208 "NaN"
4209 } else if value.is_sign_positive() {
4210 "+Inf"
4211 } else {
4212 "-Inf"
4213 }
4214}
4215
4216pub fn validate_term_feature_column_finite(
4217 data: ArrayView2<'_, f64>,
4218 term_kind: &str,
4219 term_name: &str,
4220 feature_col: usize,
4221) -> Result<(), BasisError> {
4222 let p = data.ncols();
4223 if feature_col >= p {
4224 crate::bail_dim_basis!(
4225 "{term_kind} term '{term_name}' feature column {feature_col} out of bounds for {p} columns"
4226 );
4227 }
4228 for (row, &value) in data.column(feature_col).iter().enumerate() {
4229 if !value.is_finite() {
4230 crate::bail_invalid_basis!(
4231 "{term_kind} term '{term_name}' feature column {feature_col} row {row} contains non-finite value {}",
4232 nonfinite_value_label(value)
4233 );
4234 }
4235 }
4236 Ok(())
4237}
4238
4239pub fn validate_smooth_terms_finite_inputs(
4240 data: ArrayView2<'_, f64>,
4241 terms: &[SmoothTermSpec],
4242) -> Result<(), BasisError> {
4243 for term in terms {
4244 for feature_col in smooth_term_feature_cols(term) {
4245 validate_term_feature_column_finite(data, "smooth", &term.name, feature_col)?;
4246 }
4247 }
4248 Ok(())
4249}
4250
4251pub fn validate_term_collection_finite_inputs(
4252 data: ArrayView2<'_, f64>,
4253 spec: &TermCollectionSpec,
4254) -> Result<(), BasisError> {
4255 for term in &spec.linear_terms {
4256 validate_term_feature_column_finite(data, "linear", &term.name, term.feature_col)?;
4257 }
4258 for term in &spec.random_effect_terms {
4259 validate_term_feature_column_finite(data, "random-effect", &term.name, term.feature_col)?;
4260 }
4261 validate_smooth_terms_finite_inputs(data, &spec.smooth_terms)
4262}
4263
4264#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord)]
4265pub struct JointSpatialCenterGroupKey {
4266 feature_cols: Vec<usize>,
4267 strategy_kind: CenterStrategyKind,
4268 strategy_aux: usize,
4269 requested_num_centers: usize,
4270 input_scale_bits: Option<Vec<u64>>,
4271}
4272
4273pub fn spatial_term_min_center_count(term: &SmoothTermSpec) -> usize {
4274 match &term.basis {
4275 SmoothBasisSpec::ThinPlate { feature_cols, .. } => feature_cols.len() + 1,
4276 SmoothBasisSpec::Duchon {
4277 feature_cols, spec, ..
4278 } => match spec.nullspace_order {
4279 crate::basis::DuchonNullspaceOrder::Zero => 1,
4280 crate::basis::DuchonNullspaceOrder::Linear => feature_cols.len() + 1,
4281 crate::basis::DuchonNullspaceOrder::Degree(degree) => {
4282 crate::basis::duchon_nullspace_dimension(feature_cols.len(), degree)
4283 }
4284 },
4285 SmoothBasisSpec::Matern { .. } => 1,
4286 _ => 1,
4287 }
4288}
4289
4290pub fn spatial_term_group_key(term: &SmoothTermSpec) -> Option<JointSpatialCenterGroupKey> {
4291 let (feature_cols, strategy, input_scales) = match &term.basis {
4292 SmoothBasisSpec::ThinPlate {
4293 feature_cols,
4294 spec,
4295 input_scales,
4296 } => (feature_cols, &spec.center_strategy, input_scales.as_ref()),
4297 SmoothBasisSpec::Matern {
4298 feature_cols,
4299 spec,
4300 input_scales,
4301 } => (feature_cols, &spec.center_strategy, input_scales.as_ref()),
4302 SmoothBasisSpec::Duchon {
4303 feature_cols,
4304 spec,
4305 input_scales,
4306 } => (feature_cols, &spec.center_strategy, input_scales.as_ref()),
4307 _ => return None,
4308 };
4309 let strategy_kind = center_strategy_kind(strategy);
4310 let strategy_aux = match strategy {
4311 CenterStrategy::Auto(inner) => match inner.as_ref() {
4312 CenterStrategy::KMeans { max_iter, .. } => *max_iter,
4313 CenterStrategy::UniformGrid { points_per_dim } => *points_per_dim,
4314 _ => 0,
4315 },
4316 CenterStrategy::KMeans { max_iter, .. } => *max_iter,
4317 CenterStrategy::UniformGrid { points_per_dim } => *points_per_dim,
4318 _ => 0,
4319 };
4320 Some(JointSpatialCenterGroupKey {
4321 feature_cols: feature_cols.clone(),
4322 strategy_kind,
4323 strategy_aux,
4324 requested_num_centers: center_strategy_num_centers(strategy)?,
4325 input_scale_bits: input_scales
4326 .map(|values| values.iter().map(|value| value.to_bits()).collect()),
4327 })
4328}
4329
4330pub fn spatial_term_center_strategy(term: &SmoothTermSpec) -> Option<&CenterStrategy> {
4331 match &term.basis {
4332 SmoothBasisSpec::ThinPlate { spec, .. } => Some(&spec.center_strategy),
4333 SmoothBasisSpec::Matern { spec, .. } => Some(&spec.center_strategy),
4334 SmoothBasisSpec::Duchon { spec, .. } => Some(&spec.center_strategy),
4335 _ => None,
4336 }
4337}
4338
4339pub fn set_spatial_term_centers(
4340 term: &mut SmoothTermSpec,
4341 centers: Array2<f64>,
4342) -> Result<(), BasisError> {
4343 match &mut term.basis {
4344 SmoothBasisSpec::ThinPlate { spec, .. } => {
4345 spec.center_strategy = CenterStrategy::UserProvided(centers);
4346 Ok(())
4347 }
4348 SmoothBasisSpec::Matern { spec, .. } => {
4349 spec.center_strategy = CenterStrategy::UserProvided(centers);
4350 Ok(())
4351 }
4352 SmoothBasisSpec::Duchon { spec, .. } => {
4353 spec.center_strategy = CenterStrategy::UserProvided(centers);
4354 Ok(())
4355 }
4356 _ => Err(BasisError::InvalidInput(format!(
4357 "term '{}' does not support spatial center planning",
4358 term.name
4359 ))),
4360 }
4361}
4362
4363pub fn standardized_spatial_term_data(
4364 data: ArrayView2<'_, f64>,
4365 term: &SmoothTermSpec,
4366) -> Result<Array2<f64>, BasisError> {
4367 let (feature_cols, input_scales) = match &term.basis {
4368 SmoothBasisSpec::ThinPlate {
4369 feature_cols,
4370 input_scales,
4371 ..
4372 }
4373 | SmoothBasisSpec::Matern {
4374 feature_cols,
4375 input_scales,
4376 ..
4377 }
4378 | SmoothBasisSpec::Duchon {
4379 feature_cols,
4380 input_scales,
4381 ..
4382 } => (feature_cols, input_scales.as_ref()),
4383 _ => {
4384 crate::bail_invalid_basis!("term '{}' is not a spatial smooth", term.name);
4385 }
4386 };
4387 let mut x = select_columns(data, feature_cols)?;
4388 if let Some(scales) = input_scales {
4389 apply_input_standardization(&mut x, scales);
4390 } else if let Some(scales) = compute_spatial_input_scales(x.view()) {
4391 apply_input_standardization(&mut x, &scales);
4392 }
4393 Ok(x)
4394}
4395
4396pub fn plan_joint_spatial_centers_for_term_blocks(
4397 data: ArrayView2<'_, f64>,
4398 term_blocks: &[Vec<SmoothTermSpec>],
4399) -> Result<Vec<Vec<SmoothTermSpec>>, BasisError> {
4400 let mut planned_blocks = term_blocks.to_vec();
4401 let n = data.nrows();
4402 let mut groups: BTreeMap<JointSpatialCenterGroupKey, Vec<(usize, usize)>> = BTreeMap::new();
4403
4404 for (block_idx, terms) in planned_blocks.iter().enumerate() {
4405 for (term_idx, term) in terms.iter().enumerate() {
4406 let Some(strategy) = spatial_term_center_strategy(term) else {
4407 continue;
4408 };
4409 if !center_strategy_is_auto(strategy) {
4410 continue;
4411 }
4412 let Some(group_key) = spatial_term_group_key(term) else {
4413 continue;
4414 };
4415 if !matches!(
4416 group_key.strategy_kind,
4417 CenterStrategyKind::EqualMass
4418 | CenterStrategyKind::EqualMassCovarRepresentative
4419 | CenterStrategyKind::FarthestPoint
4420 | CenterStrategyKind::KMeans
4421 ) {
4422 continue;
4423 }
4424 if center_strategy_num_centers(strategy).is_none() {
4425 continue;
4426 }
4427 groups
4428 .entry(group_key)
4429 .or_default()
4430 .push((block_idx, term_idx));
4431 }
4432 }
4433
4434 for (group_key, members) in groups {
4435 if members.len() < 2 {
4436 continue;
4437 }
4438 let min_required = members
4439 .iter()
4440 .map(|&(block_idx, term_idx)| {
4441 spatial_term_min_center_count(&planned_blocks[block_idx][term_idx])
4442 })
4443 .max()
4444 .unwrap_or(1);
4445 let joint_centers = group_key
4446 .requested_num_centers
4447 .max(min_required)
4448 .min(n.max(1));
4449 let (first_block_idx, first_term_idx) = members[0];
4450 let prototype = &planned_blocks[first_block_idx][first_term_idx];
4451 let standardized = standardized_spatial_term_data(data, prototype)?;
4452 let strategy = spatial_term_center_strategy(prototype).ok_or_else(|| {
4453 BasisError::InvalidInput(format!(
4454 "term '{}' lost its spatial center strategy during joint planning",
4455 prototype.name
4456 ))
4457 })?;
4458 let joint_strategy = center_strategy_with_num_centers(strategy, joint_centers)?;
4459 let shared_centers = select_centers_by_strategy(standardized.view(), &joint_strategy)?;
4460 log::info!(
4461 "sharing {} spatial centers across {} smooth terms over columns {:?} (requested {} centers)",
4462 shared_centers.nrows(),
4463 members.len(),
4464 group_key.feature_cols,
4465 group_key.requested_num_centers,
4466 );
4467 for (block_idx, term_idx) in members {
4468 set_spatial_term_centers(
4469 &mut planned_blocks[block_idx][term_idx],
4470 shared_centers.clone(),
4471 )?;
4472 }
4473 }
4474
4475 for block in planned_blocks.iter_mut() {
4482 for term in block.iter_mut() {
4483 auto_init_length_scale_in_place(data, term);
4484 }
4485 }
4486
4487 Ok(planned_blocks)
4488}
4489
4490const AUTO_LENGTH_SCALE_FLOOR: f64 = 1e-6;
4493
4494fn feature_columns_max_range(data: ArrayView2<'_, f64>, feature_cols: &[usize]) -> Option<f64> {
4497 let mut max_range = 0.0_f64;
4498 for &c in feature_cols {
4499 if c >= data.ncols() {
4500 continue;
4501 }
4502 let col = data.column(c);
4503 let mut lo = f64::INFINITY;
4504 let mut hi = f64::NEG_INFINITY;
4505 for &v in col.iter() {
4506 if v.is_finite() {
4507 if v < lo {
4508 lo = v;
4509 }
4510 if v > hi {
4511 hi = v;
4512 }
4513 }
4514 }
4515 if hi > lo {
4516 let r = hi - lo;
4517 if r > max_range {
4518 max_range = r;
4519 }
4520 }
4521 }
4522 if max_range.is_finite() && max_range > 0.0 {
4523 Some(max_range)
4524 } else {
4525 None
4526 }
4527}
4528
4529pub fn auto_initial_length_scale(data: ArrayView2<'_, f64>, feature_cols: &[usize]) -> f64 {
4536 let n = data.nrows();
4537 if n == 0 || feature_cols.is_empty() {
4538 return 1.0;
4539 }
4540 let Some(max_range) = feature_columns_max_range(data, feature_cols) else {
4541 return 1.0;
4542 };
4543 let init = max_range / (n as f64).sqrt();
4544 init.max(AUTO_LENGTH_SCALE_FLOOR).min(max_range)
4545}
4546
4547pub fn auto_initial_length_scale_for_centers(
4570 data: ArrayView2<'_, f64>,
4571 feature_cols: &[usize],
4572 num_centers: usize,
4573) -> f64 {
4574 let n = data.nrows();
4575 if n == 0 || feature_cols.is_empty() {
4576 return 1.0;
4577 }
4578 let Some(max_range) = feature_columns_max_range(data, feature_cols) else {
4579 return 1.0;
4580 };
4581 let resolution_points = n.max(num_centers).max(1) as f64;
4587 let spacing = max_range / resolution_points.sqrt();
4588 spacing.max(AUTO_LENGTH_SCALE_FLOOR).min(max_range)
4589}
4590
4591pub fn auto_initial_length_scale_for_low_rank_centers(
4601 data: ArrayView2<'_, f64>,
4602 feature_cols: &[usize],
4603 num_centers: usize,
4604) -> f64 {
4605 if data.nrows() == 0 || feature_cols.is_empty() {
4606 return 1.0;
4607 }
4608 let Some(max_range) = feature_columns_max_range(data, feature_cols) else {
4609 return 1.0;
4610 };
4611 let resolution_points = num_centers.max(1) as f64;
4612 let spacing = max_range / resolution_points.sqrt();
4613 spacing.max(AUTO_LENGTH_SCALE_FLOOR).min(max_range)
4614}
4615
4616fn center_strategy_requested_count(strategy: &CenterStrategy) -> Option<usize> {
4619 match strategy {
4620 CenterStrategy::Auto(inner) => center_strategy_requested_count(inner),
4621 CenterStrategy::UserProvided(centers) => Some(centers.nrows()),
4622 CenterStrategy::EqualMass { num_centers }
4623 | CenterStrategy::EqualMassCovarRepresentative { num_centers }
4624 | CenterStrategy::FarthestPoint { num_centers }
4625 | CenterStrategy::KMeans { num_centers, .. } => Some(*num_centers),
4626 CenterStrategy::UniformGrid { .. } => None,
4627 }
4628}
4629
4630pub fn auto_init_length_scale_in_place(data: ArrayView2<'_, f64>, term: &mut SmoothTermSpec) {
4634 auto_init_length_scale_in_basis(data, &mut term.basis);
4635}
4636
4637pub fn auto_init_length_scale_in_basis(data: ArrayView2<'_, f64>, basis: &mut SmoothBasisSpec) {
4650 match basis {
4651 SmoothBasisSpec::Matern {
4652 feature_cols, spec, ..
4653 } => {
4654 if spec.length_scale == 0.0 {
4655 spec.length_scale = match center_strategy_requested_count(&spec.center_strategy) {
4664 Some(k) => auto_initial_length_scale_for_centers(data, feature_cols, k),
4665 None => auto_initial_length_scale(data, feature_cols),
4666 };
4667 }
4668 }
4669 SmoothBasisSpec::ThinPlate {
4670 feature_cols, spec, ..
4671 } => {
4672 if spec.length_scale == 0.0 {
4673 spec.length_scale = match center_strategy_requested_count(&spec.center_strategy) {
4674 Some(k) => {
4675 auto_initial_length_scale_for_low_rank_centers(data, feature_cols, k)
4676 }
4677 None => auto_initial_length_scale(data, feature_cols),
4678 };
4679 }
4680 }
4681 SmoothBasisSpec::ByVariable { inner, .. }
4682 | SmoothBasisSpec::FactorSumToZero { inner, .. } => {
4683 auto_init_length_scale_in_basis(data, inner);
4684 }
4685 SmoothBasisSpec::BySmooth { smooth, .. } => {
4686 auto_init_length_scale_in_basis(data, smooth);
4687 }
4688 _ => {}
4689 }
4690}
4691
4692impl LinearFitConditioning {
4693 pub fn from_columns(design: &TermCollectionDesign, selected_cols: &[usize]) -> Self {
4694 const SCALE_EPS: f64 = 1e-12;
4695 let n = design.design.nrows();
4696 let p = design.design.ncols();
4697 let mut columns = Vec::with_capacity(selected_cols.len());
4698 if n == 0 || selected_cols.is_empty() {
4699 return Self {
4700 intercept_idx: design.intercept_range.start,
4701 columns,
4702 };
4703 }
4704 let chunk_rows = gam_linalg::utils::row_chunk_for_byte_budget(n, p);
4705 let mut sums = vec![0.0_f64; selected_cols.len()];
4711 for start in (0..n).step_by(chunk_rows) {
4712 let end = (start + chunk_rows).min(n);
4713 let chunk = design
4714 .design
4715 .try_row_chunk(start..end)
4716 .expect("LinearFitConditioning::from_columns row chunk failed");
4717 for (k, &col_idx) in selected_cols.iter().enumerate() {
4718 let column = chunk.column(col_idx);
4719 for &v in column.iter() {
4720 sums[k] += v;
4721 }
4722 }
4723 }
4724 let inv_n = 1.0_f64 / n as f64;
4725 let means: Vec<f64> = sums.iter().map(|&s| s * inv_n).collect();
4726 let mut sq_devs = vec![0.0_f64; selected_cols.len()];
4727 for start in (0..n).step_by(chunk_rows) {
4728 let end = (start + chunk_rows).min(n);
4729 let chunk = design
4730 .design
4731 .try_row_chunk(start..end)
4732 .expect("LinearFitConditioning::from_columns row chunk failed");
4733 for (k, &col_idx) in selected_cols.iter().enumerate() {
4734 let mean_k = means[k];
4735 let column = chunk.column(col_idx);
4736 for &v in column.iter() {
4737 let d = v - mean_k;
4738 sq_devs[k] += d * d;
4739 }
4740 }
4741 }
4742 for (k, &col_idx) in selected_cols.iter().enumerate() {
4743 let mean = means[k];
4744 let var = sq_devs[k] * inv_n;
4745 let (mean, scale) = if var.is_finite() && var > SCALE_EPS * SCALE_EPS {
4746 (mean, var.sqrt())
4747 } else {
4748 (0.0, 1.0)
4751 };
4752 columns.push(LinearColumnConditioning {
4753 col_idx,
4754 mean,
4755 scale,
4756 });
4757 }
4758 Self {
4759 intercept_idx: design.intercept_range.start,
4760 columns,
4761 }
4762 }
4763
4764 pub fn apply_to_design(&self, design: &Array2<f64>) -> Array2<f64> {
4765 let mut out = design.clone();
4766 for col in &self.columns {
4767 {
4768 let mut dst = out.column_mut(col.col_idx);
4769 dst -= col.mean;
4770 }
4771 if col.scale != 1.0 {
4772 out.column_mut(col.col_idx).mapv_inplace(|v| v / col.scale);
4773 }
4774 }
4775 out
4776 }
4777
4778 fn transform_matrix_columnswith_a(&self, mat: &Array2<f64>) -> Array2<f64> {
4779 let mut out = mat.clone();
4780 let intercept = self.intercept_idx;
4781 for col in &self.columns {
4782 let intercept_col = out.column(intercept).to_owned();
4783 let mut target = out.column_mut(col.col_idx);
4784 target -= &(intercept_col * col.mean);
4785 if col.scale != 1.0 {
4786 target.mapv_inplace(|v| v / col.scale);
4787 }
4788 }
4789 out
4790 }
4791
4792 fn transform_matrixrowswith_a_transpose(&self, mat: &Array2<f64>) -> Array2<f64> {
4793 let mut out = mat.clone();
4794 let intercept = self.intercept_idx;
4795 for col in &self.columns {
4796 let interceptrow = out.row(intercept).to_owned();
4797 let mut target = out.row_mut(col.col_idx);
4798 target -= &(interceptrow * col.mean);
4799 if col.scale != 1.0 {
4800 target.mapv_inplace(|v| v / col.scale);
4801 }
4802 }
4803 out
4804 }
4805
4806 fn left_multiply_by_m_inv_transpose(&self, mat_internal: &Array2<f64>) -> Array2<f64> {
4811 let mut out = mat_internal.clone();
4812 let intercept = self.intercept_idx;
4813 let interceptrow_snapshot = mat_internal.row(intercept).to_owned();
4814 for col in &self.columns {
4815 if col.scale != 1.0 {
4816 out.row_mut(col.col_idx).mapv_inplace(|v| v * col.scale);
4817 }
4818 if col.mean != 0.0 {
4819 let mut target = out.row_mut(col.col_idx);
4820 target += &(&interceptrow_snapshot * col.mean);
4821 }
4822 }
4823 out
4824 }
4825
4826 fn right_multiply_by_m_inv(&self, mat_internal: &Array2<f64>) -> Array2<f64> {
4829 let mut out = mat_internal.clone();
4830 let intercept = self.intercept_idx;
4831 let intercept_col_snapshot = mat_internal.column(intercept).to_owned();
4832 for col in &self.columns {
4833 if col.scale != 1.0 {
4834 out.column_mut(col.col_idx).mapv_inplace(|v| v * col.scale);
4835 }
4836 if col.mean != 0.0 {
4837 let mut target = out.column_mut(col.col_idx);
4838 target += &(&intercept_col_snapshot * col.mean);
4839 }
4840 }
4841 out
4842 }
4843
4844 pub fn transform_blockwise_penalties_to_internal(
4851 &self,
4852 penalties: &[BlockwisePenalty],
4853 p: usize,
4854 ) -> Vec<crate::penalty_spec::PenaltySpec> {
4855 let conditioning_cols: std::collections::HashSet<usize> =
4856 self.columns.iter().map(|c| c.col_idx).collect();
4857 penalties
4858 .iter()
4859 .map(|bp| {
4860 let overlaps =
4861 (bp.col_range.start..bp.col_range.end).any(|j| conditioning_cols.contains(&j));
4862 if overlaps {
4863 let global = bp.to_global(p);
4866 let right = self.transform_matrix_columnswith_a(&global);
4867 let transformed = self.transform_matrixrowswith_a_transpose(&right);
4868 crate::penalty_spec::PenaltySpec::Dense(transformed)
4869 } else {
4870 crate::penalty_spec::PenaltySpec::from_blockwise(bp.clone())
4873 }
4874 })
4875 .collect()
4876 }
4877
4878 pub fn backtransform_beta(&self, beta_internal: &Array1<f64>) -> Array1<f64> {
4879 let mut beta = beta_internal.clone();
4880 let intercept = self.intercept_idx;
4881 for col in &self.columns {
4882 beta[intercept] -= beta_internal[col.col_idx] * col.mean / col.scale;
4883 beta[col.col_idx] = beta_internal[col.col_idx] / col.scale;
4884 }
4885 beta
4886 }
4887
4888 pub fn transform_penalized_hessian_to_original(&self, h_internal: &Array2<f64>) -> Array2<f64> {
4891 let right = self.right_multiply_by_m_inv(h_internal);
4892 self.left_multiply_by_m_inv_transpose(&right)
4893 }
4894
4895 pub fn internal_bounds_for(&self, col_idx: usize, min: f64, max: f64) -> (f64, f64) {
4896 if let Some(col) = self.columns.iter().find(|c| c.col_idx == col_idx) {
4897 (min * col.scale, max * col.scale)
4898 } else {
4899 (min, max)
4900 }
4901 }
4902}
4903
4904pub fn freeze_raw_spatial_metadata(metadata: BasisMetadata, raw_cols: usize) -> BasisMetadata {
4905 match metadata {
4906 BasisMetadata::ThinPlate {
4907 centers,
4908 length_scale,
4909 periodic,
4910 identifiability_transform: None,
4911 input_scales,
4912 radial_reparam,
4913 } => BasisMetadata::ThinPlate {
4914 centers,
4915 length_scale,
4916 periodic,
4917 identifiability_transform: Some(Array2::eye(raw_cols)),
4918 input_scales,
4919 radial_reparam,
4920 },
4921 BasisMetadata::Duchon {
4922 centers,
4923 length_scale,
4924 periodic,
4925 power,
4926 nullspace_order,
4927 identifiability_transform: None,
4928 input_scales,
4929 aniso_log_scales,
4930 operator_collocation_points,
4931 radial_reparam,
4932 } => BasisMetadata::Duchon {
4933 centers,
4934 length_scale,
4935 periodic,
4936 power,
4937 nullspace_order,
4938 identifiability_transform: Some(Array2::eye(raw_cols)),
4939 input_scales,
4940 aniso_log_scales,
4941 operator_collocation_points,
4942 radial_reparam,
4943 },
4944 other => other,
4945 }
4946}
4947
4948pub fn matern_operator_penalty_triplet_from_metadata(
4949 metadata: &BasisMetadata,
4950) -> Result<(Vec<Array2<f64>>, Vec<usize>, Vec<PenaltyInfo>), BasisError> {
4951 let BasisMetadata::Matern {
4952 centers,
4953 length_scale,
4954 periodic,
4955 nu,
4956 include_intercept,
4957 identifiability_transform,
4958 aniso_log_scales,
4959 input_scales,
4960 ..
4961 } = metadata
4962 else {
4963 crate::bail_invalid_basis!("Matérn operator penalties require Matérn metadata");
4964 };
4965 let penalty_length_scale = match input_scales.as_deref() {
4977 Some(scales) => compensate_length_scale_for_standardization(*length_scale, scales),
4978 None => *length_scale,
4979 };
4980 matern_operator_penalty_triplet_at_length_scale(
4981 centers.view(),
4982 periodic.as_deref(),
4983 identifiability_transform.as_ref(),
4984 *nu,
4985 *include_intercept,
4986 aniso_log_scales.as_deref(),
4987 penalty_length_scale,
4988 )
4989}
4990
4991pub fn matern_operator_penalty_triplet_at_length_scale(
5009 centers: ArrayView2<'_, f64>,
5010 periodic: Option<&[Option<f64>]>,
5011 identifiability_transform: Option<&Array2<f64>>,
5012 nu: crate::basis::MaternNu,
5013 include_intercept: bool,
5014 aniso_log_scales: Option<&[f64]>,
5015 effective_length_scale: f64,
5016) -> Result<(Vec<Array2<f64>>, Vec<usize>, Vec<PenaltyInfo>), BasisError> {
5017 let penalty_centers = crate::basis::expand_periodic_centers(¢ers.to_owned(), periodic)?;
5018 let ops = build_matern_collocation_operator_matrices(
5019 penalty_centers.view(),
5020 None,
5021 effective_length_scale,
5022 nu,
5023 include_intercept,
5024 identifiability_transform.map(|z| z.view()),
5025 aniso_log_scales,
5026 )?;
5027 const ORDER_EPS: f64 = 1e-9;
5035 let d = penalty_centers.ncols();
5036 let m = nu.half_integer_value() + 0.5 * d as f64;
5037 let mut candidates = Vec::with_capacity(3);
5038 for (raw, source, min_order) in [
5039 (ops.d0.t().dot(&ops.d0), PenaltySource::OperatorMass, 0.0),
5040 (ops.d1.t().dot(&ops.d1), PenaltySource::OperatorTension, 1.0),
5041 (
5042 ops.d2.t().dot(&ops.d2),
5043 PenaltySource::OperatorStiffness,
5044 2.0,
5045 ),
5046 ] {
5047 if min_order > 0.0 && m <= min_order + ORDER_EPS {
5048 continue;
5049 }
5050 let sym = (&raw + &raw.t()) * 0.5;
5051 let (matrix, normalization_scale) = normalize_penalty_in_constrained_space(&sym);
5052 candidates.push(PenaltyCandidate {
5053 matrix,
5054 nullspace_dim_hint: 0,
5055 source,
5056 normalization_scale,
5057 kronecker_factors: None,
5058 op: None,
5059 });
5060 }
5061 filter_active_penalty_candidates(candidates)
5062}
5063
5064pub fn normalize_penalty_in_constrained_space(matrix: &Array2<f64>) -> (Array2<f64>, f64) {
5065 let matrix = (matrix + &matrix.t().to_owned()) * 0.5;
5070 let matrix = crate::basis::project_penalty_to_psd_cone(&matrix);
5072 let c = matrix.iter().map(|v| v * v).sum::<f64>().sqrt();
5073 if c.is_finite() && c > 0.0 {
5074 (matrix.mapv(|v| v / c), c)
5075 } else {
5076 (matrix, 1.0)
5077 }
5078}
5079
5080pub fn tensor_product_design_from_sparse_marginals(
5081 marginal_sparse: &[&SparseColMat<usize, f64>],
5082) -> Result<SparseColMat<usize, f64>, BasisError> {
5083 if marginal_sparse.is_empty() {
5084 crate::bail_invalid_basis!("TensorBSpline requires at least one marginal basis");
5085 }
5086 let n = marginal_sparse[0].nrows();
5087 for (i, m) in marginal_sparse.iter().enumerate().skip(1) {
5088 if m.nrows() != n {
5089 crate::bail_dim_basis!(
5090 "tensor sparse marginal row mismatch at dim {i}: expected {n}, got {}",
5091 m.nrows()
5092 );
5093 }
5094 }
5095 let dims: Vec<usize> = marginal_sparse.iter().map(|m| m.ncols()).collect();
5096 let total_cols = dims.iter().try_fold(1usize, |acc, &q| {
5097 acc.checked_mul(q)
5098 .ok_or_else(|| BasisError::DimensionMismatch("tensor basis too large".to_string()))
5099 })?;
5100 let mut strides = vec![1usize; dims.len()];
5101 for d in (0..dims.len().saturating_sub(1)).rev() {
5102 strides[d] = strides[d + 1]
5103 .checked_mul(dims[d + 1])
5104 .ok_or_else(|| BasisError::DimensionMismatch("tensor basis too large".to_string()))?;
5105 }
5106
5107 use faer::sparse::SparseRowMat;
5108 let csrs: Vec<SparseRowMat<usize, f64>> = marginal_sparse
5109 .iter()
5110 .enumerate()
5111 .map(|(d, m)| {
5112 m.as_ref().to_row_major().map_err(|e| {
5113 BasisError::SparseCreation(format!(
5114 "tensor sparse marginal {d} CSR conversion failed: {e:?}"
5115 ))
5116 })
5117 })
5118 .collect::<Result<Vec<_>, _>>()?;
5119 let row_ptrs: Vec<&[usize]> = csrs.iter().map(|c| c.symbolic().row_ptr()).collect();
5120 let col_idxs: Vec<&[usize]> = csrs.iter().map(|c| c.symbolic().col_idx()).collect();
5121 let vals: Vec<&[f64]> = csrs.iter().map(|c| c.val()).collect();
5122
5123 use rayon::prelude::*;
5124 const CHUNK: usize = 1024;
5125 let num_chunks = n.div_ceil(CHUNK);
5126 let per_chunk: Vec<Vec<Triplet<usize, usize, f64>>> = (0..num_chunks)
5127 .into_par_iter()
5128 .map(|chunk_idx| {
5129 let row_start = chunk_idx * CHUNK;
5130 let row_end = (row_start + CHUNK).min(n);
5131 let mut chunk_triplets = Vec::<Triplet<usize, usize, f64>>::new();
5132 let mut cur_cols = Vec::<usize>::with_capacity(64);
5133 let mut cur_vals = Vec::<f64>::with_capacity(64);
5134 let mut next_cols = Vec::<usize>::with_capacity(64);
5135 let mut next_vals = Vec::<f64>::with_capacity(64);
5136 for i in row_start..row_end {
5137 cur_cols.clear();
5138 cur_vals.clear();
5139 cur_cols.push(0);
5140 cur_vals.push(1.0);
5141 let mut row_is_zero = false;
5142 for d in 0..dims.len() {
5143 let row_start_d = row_ptrs[d][i];
5144 let row_end_d = row_ptrs[d][i + 1];
5145 if row_start_d == row_end_d {
5146 row_is_zero = true;
5147 break;
5148 }
5149 let stride = strides[d];
5150 next_cols.clear();
5151 next_vals.clear();
5152 next_cols.reserve(cur_cols.len() * (row_end_d - row_start_d));
5153 next_vals.reserve(cur_vals.len() * (row_end_d - row_start_d));
5154 for (&prev_col, &prev_val) in cur_cols.iter().zip(cur_vals.iter()) {
5155 for ptr in row_start_d..row_end_d {
5156 let cj = col_idxs[d][ptr];
5157 let vj = vals[d][ptr];
5158 next_cols.push(prev_col + cj * stride);
5159 next_vals.push(prev_val * vj);
5160 }
5161 }
5162 std::mem::swap(&mut cur_cols, &mut next_cols);
5163 std::mem::swap(&mut cur_vals, &mut next_vals);
5164 }
5165 if row_is_zero {
5166 continue;
5167 }
5168 for (&col, &val) in cur_cols.iter().zip(cur_vals.iter()) {
5169 chunk_triplets.push(Triplet::new(i, col, val));
5170 }
5171 }
5172 chunk_triplets
5173 })
5174 .collect();
5175 let total_nnz: usize = per_chunk.iter().map(Vec::len).sum();
5176 let mut triplets = Vec::<Triplet<usize, usize, f64>>::with_capacity(total_nnz);
5177 for chunk in per_chunk {
5178 triplets.extend(chunk);
5179 }
5180 SparseColMat::try_new_from_triplets(n, total_cols, &triplets).map_err(|e| {
5181 BasisError::SparseCreation(format!(
5182 "failed to assemble sparse tensor product design: {e:?}"
5183 ))
5184 })
5185}
5186
5187pub fn dense_local_margin_to_sparse(
5188 dense: &Array2<f64>,
5189) -> Result<SparseColMat<usize, f64>, BasisError> {
5190 let expected_row_nnz = dense.ncols().min(4);
5191 let mut triplets =
5192 Vec::<Triplet<usize, usize, f64>>::with_capacity(dense.nrows() * expected_row_nnz);
5193 for ((row, col), &value) in dense.indexed_iter() {
5194 if value != 0.0 {
5195 triplets.push(Triplet::new(row, col, value));
5196 }
5197 }
5198 SparseColMat::try_new_from_triplets(dense.nrows(), dense.ncols(), &triplets).map_err(|e| {
5199 BasisError::SparseCreation(format!(
5200 "failed to convert tensor marginal design to sparse form: {e:?}"
5201 ))
5202 })
5203}
5204
5205pub struct TensorMarginRangeNullProjectors {
5206 range: Array2<f64>,
5207 null: Array2<f64>,
5208}
5209
5210pub fn projector_from_columns(columns: &Array2<f64>, indices: &[usize]) -> Array2<f64> {
5211 if indices.is_empty() {
5212 return Array2::<f64>::zeros((columns.nrows(), columns.nrows()));
5213 }
5214 let basis = columns.select(Axis(1), indices);
5215 basis.dot(&basis.t())
5216}
5217
5218pub fn tensor_margin_range_null_projectors(
5219 normalized_marginal_penalties: &[(Array2<f64>, f64)],
5220) -> Result<Vec<TensorMarginRangeNullProjectors>, BasisError> {
5221 normalized_marginal_penalties
5222 .iter()
5223 .enumerate()
5224 .map(|(dim, (penalty, _))| {
5225 let analysis = crate::basis::analyze_penalty_block(penalty)?;
5226 if analysis.rank == 0 {
5227 crate::bail_invalid_basis!(
5228 "t2 separable tensor penalty margin {dim} has rank-zero penalty; \
5229 cannot split penalized and null subspaces"
5230 );
5231 }
5232 let mut range_idx = Vec::<usize>::new();
5233 let mut null_idx = Vec::<usize>::new();
5234 for (idx, &ev) in analysis.eigenvalues.iter().enumerate() {
5235 if ev > analysis.tol {
5236 range_idx.push(idx);
5237 } else {
5238 null_idx.push(idx);
5239 }
5240 }
5241 Ok(TensorMarginRangeNullProjectors {
5242 range: projector_from_columns(&analysis.eigenvectors, &range_idx),
5243 null: projector_from_columns(&analysis.eigenvectors, &null_idx),
5244 })
5245 })
5246 .collect()
5247}
5248
5249pub fn build_tensor_bspline_basis(
5250 data: ArrayView2<'_, f64>,
5251 feature_cols: &[usize],
5252 spec: &TensorBSplineSpec,
5253) -> Result<BasisBuildResult, BasisError> {
5254 if feature_cols.is_empty() {
5255 crate::bail_invalid_basis!("TensorBSpline requires at least one feature column");
5256 }
5257 if feature_cols.len() != spec.marginalspecs.len() {
5258 crate::bail_dim_basis!(
5259 "TensorBSpline feature/spec mismatch: feature_cols={}, marginalspecs={}",
5260 feature_cols.len(),
5261 spec.marginalspecs.len()
5262 );
5263 }
5264 if !spec.periods.is_empty() && spec.periods.len() != feature_cols.len() {
5265 crate::bail_dim_basis!(
5266 "TensorBSpline periods length {} does not match feature count {}",
5267 spec.periods.len(),
5268 feature_cols.len()
5269 );
5270 }
5271 let p = data.ncols();
5272 for &c in feature_cols {
5273 if c >= p {
5274 crate::bail_dim_basis!(
5275 "tensor feature column {c} is out of bounds for data with {p} columns"
5276 );
5277 }
5278 }
5279
5280 let mut marginal_knots = Vec::<Array1<f64>>::with_capacity(feature_cols.len());
5281 let mut marginal_is_cr_flags = Vec::<bool>::with_capacity(feature_cols.len());
5284 let mut marginal_degrees = Vec::<usize>::with_capacity(feature_cols.len());
5285 let mut marginalnum_basis = Vec::<usize>::with_capacity(feature_cols.len());
5286 let mut marginal_penalties = Vec::<Array2<f64>>::with_capacity(feature_cols.len());
5287 let mut marginal_designs = Vec::<Array2<f64>>::with_capacity(feature_cols.len());
5288 let mut marginal_effective_periods = Vec::<Option<f64>>::with_capacity(feature_cols.len());
5296 let mut marginal_sparse =
5304 Vec::<Option<SparseColMat<usize, f64>>>::with_capacity(feature_cols.len());
5305
5306 for (dim, (&col, marginalspec)) in feature_cols
5309 .iter()
5310 .zip(spec.marginalspecs.iter())
5311 .enumerate()
5312 {
5313 let mut marginal_unconstrained = marginalspec.clone();
5318 marginal_unconstrained.identifiability = BSplineIdentifiability::None;
5319 let built = build_bspline_basis_1d(data.column(col), &marginal_unconstrained)?;
5320 let (knots, marginal_is_cr) = match built.metadata {
5325 BasisMetadata::BSpline1D { knots, .. } => (knots, false),
5326 BasisMetadata::CubicRegression1D { knots, .. } => (knots, true),
5327 _ => {
5328 crate::bail_invalid_basis!(
5329 "internal TensorBSpline error at dim {dim}: expected BSpline1D or CubicRegression1D metadata"
5330 );
5331 }
5332 };
5333 let metadata_knots = match marginalspec.knotspec {
5334 BSplineKnotSpec::PeriodicUniform {
5335 data_range,
5336 num_basis,
5337 } => Array1::linspace(data_range.0, data_range.1, num_basis),
5338 _ => knots,
5339 };
5340 marginal_knots.push(metadata_knots);
5341 marginal_is_cr_flags.push(marginal_is_cr);
5342 marginal_degrees.push(marginalspec.degree);
5343 marginalnum_basis.push(built.design.ncols());
5344 let dense_marginal = built.design.to_dense();
5349 let sparse_view: Option<SparseColMat<usize, f64>> = match built.design.as_sparse() {
5350 Some(sd) => {
5351 let inner: &SparseColMat<usize, f64> = sd;
5352 Some(inner.clone())
5353 }
5354 None => match marginalspec.knotspec {
5355 BSplineKnotSpec::PeriodicUniform { .. } => {
5356 Some(dense_local_margin_to_sparse(&dense_marginal)?)
5357 }
5358 _ => None,
5359 },
5360 };
5361 marginal_sparse.push(sparse_view);
5362 marginal_designs.push(dense_marginal);
5363 marginal_penalties.push(
5364 built
5365 .penalties
5366 .first()
5367 .ok_or_else(|| {
5368 BasisError::InvalidInput(format!(
5369 "internal TensorBSpline error at dim {dim}: missing marginal penalty"
5370 ))
5371 })?
5372 .clone(),
5373 );
5374 built.nullspace_dims.first().ok_or_else(|| {
5375 BasisError::InvalidInput(format!(
5376 "internal TensorBSpline error at dim {dim}: missing marginal nullspace dim"
5377 ))
5378 })?;
5379 let implied_period = match marginalspec.knotspec {
5387 BSplineKnotSpec::PeriodicUniform { data_range, .. } => {
5388 Some(data_range.1 - data_range.0)
5389 }
5390 _ => spec.periods.get(dim).and_then(|p| *p),
5391 };
5392 marginal_effective_periods.push(implied_period);
5393 }
5394
5395 let total_cols: usize = marginalnum_basis.iter().product();
5396 let mut dense_design = (!matches!(spec.identifiability, TensorBSplineIdentifiability::None))
5397 .then(|| tensor_product_design_from_marginals(&marginal_designs))
5398 .transpose()?;
5399 let mut candidates = Vec::<PenaltyCandidate>::with_capacity(
5400 match spec.penalty_decomposition {
5401 TensorBSplinePenaltyDecomposition::MarginalKroneckerSum => marginal_penalties.len(),
5402 TensorBSplinePenaltyDecomposition::Separable => marginal_penalties.len() * 2,
5403 } + if spec.double_penalty { 1 } else { 0 },
5404 );
5405
5406 let normalized_marginal_penalties: Vec<(Array2<f64>, f64)> = marginal_penalties
5414 .iter()
5415 .map(normalize_penalty_in_constrained_space)
5416 .collect();
5417 let mut kronecker_marginal_penalties =
5418 Vec::<Array2<f64>>::with_capacity(normalized_marginal_penalties.len());
5419
5420 match spec.penalty_decomposition {
5421 TensorBSplinePenaltyDecomposition::MarginalKroneckerSum => {
5422 let mut marginal_kron_sum = Array2::<f64>::zeros((total_cols, total_cols));
5428
5429 for dim in 0..normalized_marginal_penalties.len() {
5430 let mut s_dim = Array2::<f64>::eye(1);
5431 let mut factors = Vec::<Array2<f64>>::with_capacity(marginalnum_basis.len());
5432 for (j, &qj) in marginalnum_basis.iter().enumerate() {
5433 let factor = if j == dim {
5434 normalized_marginal_penalties[j].0.clone()
5435 } else {
5436 Array2::<f64>::eye(qj)
5437 };
5438 factors.push(factor.clone());
5439 s_dim = kronecker_product(&s_dim, &factor);
5440 }
5441 if dim == kronecker_marginal_penalties.len() {
5442 kronecker_marginal_penalties.push(normalized_marginal_penalties[dim].0.clone());
5443 }
5444 marginal_kron_sum += &s_dim;
5445
5446 candidates.push(PenaltyCandidate {
5447 matrix: s_dim,
5448 nullspace_dim_hint: 0,
5449 source: PenaltySource::TensorMarginal { dim },
5450 normalization_scale: normalized_marginal_penalties[dim].1,
5451 kronecker_factors: Some(factors),
5452 op: None,
5453 });
5454 }
5455
5456 if spec.double_penalty
5457 && let Some(shrink) =
5458 crate::basis::build_nullspace_shrinkage_penalty(&marginal_kron_sum)?
5459 {
5460 let (matrix, normalization_scale) =
5461 normalize_penalty_in_constrained_space(&shrink.sym_penalty);
5462 candidates.push(PenaltyCandidate {
5463 matrix,
5464 nullspace_dim_hint: 0,
5465 source: PenaltySource::TensorGlobalRidge,
5466 normalization_scale,
5467 kronecker_factors: None,
5468 op: None,
5469 });
5470 }
5471 }
5472 TensorBSplinePenaltyDecomposition::Separable => {
5473 let projectors = tensor_margin_range_null_projectors(&normalized_marginal_penalties)?;
5474 let n_masks = 1usize.checked_shl(projectors.len() as u32).ok_or_else(|| {
5475 BasisError::InvalidInput(format!(
5476 "t2 separable tensor penalty supports at most {} margins, got {}",
5477 usize::BITS - 1,
5478 projectors.len()
5479 ))
5480 })?;
5481 for mask in 1..n_masks {
5482 let mut matrix = Array2::<f64>::eye(1);
5483 let mut factors = Vec::<Array2<f64>>::with_capacity(projectors.len());
5484 let mut penalized_margins = Vec::<usize>::new();
5485 for (dim, projector) in projectors.iter().enumerate() {
5486 let use_range = ((mask >> dim) & 1) == 1;
5487 let factor = if use_range {
5488 penalized_margins.push(dim);
5489 projector.range.clone()
5490 } else {
5491 projector.null.clone()
5492 };
5493 matrix = kronecker_product(&matrix, &factor);
5494 factors.push(factor);
5495 }
5496 let (matrix, normalization_scale) = normalize_penalty_in_constrained_space(&matrix);
5497 candidates.push(PenaltyCandidate {
5498 matrix,
5499 nullspace_dim_hint: 0,
5500 source: PenaltySource::TensorSeparable { penalized_margins },
5501 normalization_scale,
5502 kronecker_factors: Some(factors),
5503 op: None,
5504 });
5505 }
5506
5507 if spec.double_penalty {
5508 let mut matrix = Array2::<f64>::eye(1);
5509 let mut factors = Vec::<Array2<f64>>::with_capacity(projectors.len());
5510 for projector in &projectors {
5511 matrix = kronecker_product(&matrix, &projector.null);
5512 factors.push(projector.null.clone());
5513 }
5514 let (matrix, normalization_scale) = normalize_penalty_in_constrained_space(&matrix);
5515 candidates.push(PenaltyCandidate {
5516 matrix,
5517 nullspace_dim_hint: 0,
5518 source: PenaltySource::TensorGlobalRidge,
5519 normalization_scale,
5520 kronecker_factors: Some(factors),
5521 op: None,
5522 });
5523 }
5524 }
5525 }
5526
5527 let z_opt = match &spec.identifiability {
5528 TensorBSplineIdentifiability::None => None,
5529 TensorBSplineIdentifiability::SumToZero => {
5530 if total_cols < 2 {
5531 crate::bail_invalid_basis!(
5532 "TensorBSpline requires at least 2 basis coefficients to enforce sum-to-zero identifiability"
5533 );
5534 }
5535 let dense_design_ref = dense_design.as_ref().ok_or_else(|| {
5536 BasisError::InvalidInput(
5537 "tensor sum-to-zero identifiability requires a realized basis".to_string(),
5538 )
5539 })?;
5540 let (_, z) = apply_sum_to_zero_constraint(dense_design_ref.view(), None)?;
5541 let gauge = gam_problem::Gauge::sum_to_zero(z);
5542 Some(gauge.block_transform(0))
5543 }
5544 TensorBSplineIdentifiability::MarginalSumToZero => {
5545 if marginal_designs.len() < 2 {
5556 crate::bail_invalid_basis!(
5557 "tensor interaction (ti) identifiability requires at least 2 margins"
5558 );
5559 }
5560 let mut z = Array2::<f64>::eye(1);
5561 for (dim, marginal) in marginal_designs.iter().enumerate() {
5562 if marginal.ncols() < 2 {
5563 crate::bail_invalid_basis!(
5564 "tensor interaction (ti) margin {dim} has fewer than 2 basis functions; \
5565 cannot remove its marginal main effect"
5566 );
5567 }
5568 let (_, z_dim) = apply_sum_to_zero_constraint(marginal.view(), None)?;
5569 let gauge_dim = gam_problem::Gauge::sum_to_zero(z_dim);
5570 let z_dim = gauge_dim.block_transform(0);
5571 z = kronecker_product(&z, &z_dim);
5572 }
5573 Some(z)
5574 }
5575 TensorBSplineIdentifiability::FrozenTransform { transform } => {
5576 if transform.nrows() != total_cols {
5577 crate::bail_dim_basis!(
5578 "frozen tensor identifiability transform mismatch: design has {} columns but transform has {} rows",
5579 total_cols,
5580 transform.nrows()
5581 );
5582 }
5583 Some(transform.clone())
5584 }
5585 };
5586
5587 if let Some(z) = z_opt.as_ref() {
5588 let gauge = gam_problem::Gauge::from_block_transforms(&[z.clone()]);
5589 let dense = dense_design.as_mut().ok_or_else(|| {
5590 BasisError::InvalidInput(
5591 "tensor identifiability transform requires a realized basis".to_string(),
5592 )
5593 })?;
5594 let restricted_design = gauge.restrict_design(dense);
5595 *dense = restricted_design;
5596 candidates = candidates
5597 .into_iter()
5598 .map(|candidate| -> Result<PenaltyCandidate, BasisError> {
5599 let matrix = gauge.restrict_penalty(&candidate.matrix);
5600 let (matrix, c_new) = normalize_penalty_in_constrained_space(&matrix);
5608 Ok(PenaltyCandidate {
5609 nullspace_dim_hint: candidate.nullspace_dim_hint,
5610 matrix,
5611 source: candidate.source,
5612 normalization_scale: candidate.normalization_scale * c_new,
5613 kronecker_factors: None,
5619 op: candidate.op.clone(),
5620 })
5621 })
5622 .collect::<Result<Vec<_>, _>>()?;
5623 }
5624
5625 let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
5626 filter_active_penalty_candidates_with_ops(candidates)?;
5627 let identifiability_is_none =
5628 matches!(spec.identifiability, TensorBSplineIdentifiability::None);
5629 let all_marginals_sparse = marginal_sparse.iter().all(Option::is_some);
5637 let design = if let Some(dense_design) = dense_design {
5638 DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(dense_design))
5639 } else if identifiability_is_none && all_marginals_sparse {
5640 let sparse_marginals: Vec<&SparseColMat<usize, f64>> = marginal_sparse
5646 .iter()
5647 .map(|m| m.as_ref().expect("all_marginals_sparse just verified"))
5648 .collect();
5649 let sparse_design = tensor_product_design_from_sparse_marginals(&sparse_marginals)?;
5650 DesignMatrix::Sparse(gam_linalg::matrix::SparseDesignMatrix::new(sparse_design))
5651 } else {
5652 let marginals: Vec<Arc<Array2<f64>>> = marginal_designs
5653 .iter()
5654 .map(|m| Arc::new(m.clone()))
5655 .collect();
5656 let op = TensorProductDesignOperator::new(marginals).map_err(|e| {
5657 BasisError::InvalidInput(format!("TensorProductDesignOperator build failed: {e}"))
5658 })?;
5659 DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(Arc::new(op)))
5660 };
5661
5662 Ok(BasisBuildResult {
5663 design,
5664 penalties,
5665 nullspace_dims,
5666 penaltyinfo,
5667 ops,
5668 null_eigenvectors,
5669 joint_null_rotation: None,
5670 metadata: BasisMetadata::TensorBSpline {
5671 feature_cols: feature_cols.to_vec(),
5672 knots: marginal_knots,
5673 degrees: marginal_degrees,
5674 periods: marginal_effective_periods,
5681 is_cr: marginal_is_cr_flags,
5682 identifiability_transform: z_opt,
5683 },
5684 kronecker_factored: if matches!(spec.identifiability, TensorBSplineIdentifiability::None)
5685 && matches!(
5686 spec.penalty_decomposition,
5687 TensorBSplinePenaltyDecomposition::MarginalKroneckerSum
5688 ) {
5689 Some(KroneckerFactoredBasis::new(
5690 marginal_designs,
5691 kronecker_marginal_penalties,
5692 marginalnum_basis.clone(),
5693 spec.double_penalty,
5694 ))
5695 } else {
5696 None
5697 },
5698 })
5699}
5700
5701pub fn tensor_product_design_from_marginals(
5702 marginal_designs: &[Array2<f64>],
5703) -> Result<Array2<f64>, BasisError> {
5704 if marginal_designs.is_empty() {
5705 crate::bail_invalid_basis!("TensorBSpline requires at least one marginal basis");
5706 }
5707 let n = marginal_designs[0].nrows();
5708 for (i, b) in marginal_designs.iter().enumerate().skip(1) {
5709 if b.nrows() != n {
5710 crate::bail_dim_basis!(
5711 "tensor marginal row mismatch at dim {i}: expected {n}, got {}",
5712 b.nrows()
5713 );
5714 }
5715 }
5716 let total_cols = marginal_designs.iter().try_fold(1usize, |acc, b| {
5717 acc.checked_mul(b.ncols())
5718 .ok_or_else(|| BasisError::DimensionMismatch("tensor basis too large".to_string()))
5719 })?;
5720 use ndarray::parallel::prelude::*;
5726 use rayon::iter::{IntoParallelIterator, ParallelIterator};
5727 let mut design = Array2::<f64>::zeros((n, total_cols));
5728 design
5729 .axis_chunks_iter_mut(ndarray::Axis(0), 1024)
5730 .into_par_iter()
5731 .enumerate()
5732 .for_each(|(chunk_idx, mut block)| {
5733 let row_offset = chunk_idx * 1024;
5734 let mut cur = Vec::<f64>::with_capacity(total_cols);
5736 let mut next = Vec::<f64>::with_capacity(total_cols);
5737 for (local_i, mut out_row) in block.outer_iter_mut().enumerate() {
5738 let i = row_offset + local_i;
5739 cur.clear();
5740 cur.push(1.0);
5741 for b in marginal_designs {
5742 let q = b.ncols();
5743 next.clear();
5744 next.resize(cur.len() * q, 0.0);
5745 let b_row = b.row(i);
5749 let b_slice = b_row
5750 .as_slice()
5751 .expect("Array2 row from outer_iter is contiguous");
5752 for (a_idx, &aval) in cur.iter().enumerate() {
5753 let off = a_idx * q;
5754 let dst = &mut next[off..off + q];
5755 for col in 0..q {
5756 dst[col] = aval * b_slice[col];
5757 }
5758 }
5759 std::mem::swap(&mut cur, &mut next);
5760 }
5761 let out_slice = out_row
5766 .as_slice_mut()
5767 .expect("design row is contiguous in C-major Array2");
5768 out_slice.copy_from_slice(&cur);
5769 }
5770 });
5771 Ok(design)
5772}
5773
5774fn fmt_level_value(v: f64) -> String {
5778 if v.is_finite() && v.fract() == 0.0 && v.abs() < 1e15 {
5779 format!("{}", v as i64)
5780 } else {
5781 format!("{v}")
5782 }
5783}
5784
5785pub fn build_random_effect_block(
5786 data: ArrayView2<'_, f64>,
5787 spec: &RandomEffectTermSpec,
5788) -> Result<RandomEffectBlock, BasisError> {
5789 let n = data.nrows();
5790 let p = data.ncols();
5791 if spec.feature_col >= p {
5792 crate::bail_dim_basis!(
5793 "random-effect term '{}' feature column {} out of bounds for {} columns",
5794 spec.name,
5795 spec.feature_col,
5796 p
5797 );
5798 }
5799
5800 let col = data.column(spec.feature_col);
5801 if col.iter().any(|v| !v.is_finite()) {
5802 crate::bail_invalid_basis!(
5803 "random-effect term '{}' contains non-finite group values",
5804 spec.name
5805 );
5806 }
5807
5808 let kept_levels: Vec<u64> = if let Some(levels) = spec.frozen_levels.as_ref() {
5809 if levels.is_empty() {
5810 crate::bail_invalid_basis!(
5811 "random-effect term '{}' has empty frozen_levels",
5812 spec.name
5813 );
5814 }
5815 levels
5819 .iter()
5820 .map(|&b| gam_data::canonical_level_bits(f64::from_bits(b)))
5821 .collect()
5822 } else {
5823 let mut seen = BTreeSet::<u64>::new();
5824 let mut levels = Vec::<u64>::new();
5825 for &v in col {
5826 let bits = gam_data::canonical_level_bits(v);
5827 if seen.insert(bits) {
5828 levels.push(bits);
5829 }
5830 }
5831 if levels.is_empty() {
5832 crate::bail_invalid_basis!("random-effect term '{}' has no observed levels", spec.name);
5833 }
5834 let start_idx = if spec.drop_first_level && levels.len() > 1 {
5835 1usize
5836 } else {
5837 0usize
5838 };
5839 levels[start_idx..].to_vec()
5840 };
5841
5842 if kept_levels.is_empty() {
5843 crate::bail_invalid_basis!(
5844 "random-effect term '{}' drops all levels; keep at least one level",
5845 spec.name
5846 );
5847 }
5848
5849 let q = kept_levels.len();
5850 let mut level_to_col = BTreeMap::<u64, usize>::new();
5851 for (idx, &bits) in kept_levels.iter().enumerate() {
5852 if level_to_col.insert(bits, idx).is_some() {
5853 crate::bail_invalid_basis!(
5854 "random-effect term '{}' has duplicate frozen level bits {bits}",
5855 spec.name
5856 );
5857 }
5858 }
5859 let strict_unseen = !spec.lenient_unseen && !spec.drop_first_level && spec.frozen_levels.is_some();
5873 let mut group_ids = Vec::with_capacity(n);
5874 for (row, &v) in col.iter().enumerate() {
5875 let bits = gam_data::canonical_level_bits(v);
5876 let group_id = level_to_col.get(&bits).copied();
5877 if strict_unseen && group_id.is_none() {
5878 crate::bail_invalid_basis!(
5879 "unseen level '{}' in fixed factor column '{}' at row {}; the factor's levels \
5880 were fixed at fit time and an out-of-vocabulary level cannot be predicted \
5881 (use group({}) for a random effect that tolerates held-out levels)",
5882 fmt_level_value(v),
5883 spec.name,
5884 row,
5885 spec.name
5886 );
5887 }
5888 group_ids.push(group_id);
5889 }
5890
5891 Ok(RandomEffectBlock {
5892 name: spec.name.clone(),
5893 group_ids,
5894 num_groups: q,
5895 kept_levels,
5896 })
5897}
5898
5899#[cfg(test)]
5900mod random_effect_signed_zero_tests {
5901 use super::{RandomEffectTermSpec, build_random_effect_block};
5902 use ndarray::array;
5903
5904 fn spec() -> RandomEffectTermSpec {
5905 RandomEffectTermSpec {
5906 name: "g".to_string(),
5907 feature_col: 0,
5908 drop_first_level: false,
5909 penalized: true,
5910 frozen_levels: None,
5911 lenient_unseen: true,
5912 }
5913 }
5914
5915 #[test]
5916 fn signed_zero_rows_share_one_group() {
5917 let data = array![[-0.0_f64], [0.0], [1.0], [-0.0], [1.0]];
5921 let block = build_random_effect_block(data.view(), &spec()).unwrap();
5922 assert_eq!(block.num_groups, 2, "0.0/-0.0 must not split into two groups");
5923 assert_eq!(block.group_ids[0], block.group_ids[1]);
5925 assert_eq!(block.group_ids[0], block.group_ids[3]);
5926 assert_eq!(block.group_ids[2], block.group_ids[4]);
5927 assert_ne!(block.group_ids[0], block.group_ids[2]);
5928 }
5929
5930 #[test]
5931 fn frozen_positive_zero_matches_negative_zero_row() {
5932 let mut s = spec();
5935 s.frozen_levels = Some(vec![0.0_f64.to_bits(), 1.0_f64.to_bits()]);
5936 let data = array![[-0.0_f64], [1.0]];
5937 let block = build_random_effect_block(data.view(), &s).unwrap();
5938 assert_eq!(block.group_ids[0], Some(0), "-0.0 must match the +0.0 column");
5939 assert_eq!(block.group_ids[1], Some(1));
5940 }
5941
5942 #[test]
5943 fn frozen_negative_zero_matches_positive_zero_row() {
5944 let mut s = spec();
5947 s.frozen_levels = Some(vec![(-0.0_f64).to_bits(), 1.0_f64.to_bits()]);
5948 let data = array![[0.0_f64], [1.0]];
5949 let block = build_random_effect_block(data.view(), &s).unwrap();
5950 assert_eq!(block.group_ids[0], Some(0), "+0.0 must match the -0.0 column");
5951 }
5952
5953 fn fixed_factor_spec() -> RandomEffectTermSpec {
5956 let mut s = spec();
5959 s.name = "year".to_string();
5960 s.lenient_unseen = false;
5961 s
5962 }
5963
5964 #[test]
5965 fn fixed_factor_rejects_unseen_numeric_level_at_predict() {
5966 let mut s = fixed_factor_spec();
5971 s.frozen_levels = Some(vec![2000.0_f64.to_bits(), 2001.0_f64.to_bits()]);
5972 let data = array![[2000.0_f64], [1999.0]];
5973 let err = build_random_effect_block(data.view(), &s)
5974 .expect_err("an unseen fixed-factor level must be rejected");
5975 let msg = format!("{err}");
5976 assert!(msg.contains("unseen level"), "message must name the defect: {msg}");
5977 assert!(msg.contains("1999"), "message must name the integer level (not 1999.0): {msg}");
5978 assert!(msg.contains("year"), "message must name the column: {msg}");
5979 }
5980
5981 #[test]
5982 fn fixed_factor_accepts_seen_numeric_levels_at_predict() {
5983 let mut s = fixed_factor_spec();
5986 s.frozen_levels = Some(vec![2000.0_f64.to_bits(), 2001.0_f64.to_bits()]);
5987 let data = array![[2001.0_f64], [2000.0]];
5988 let block = build_random_effect_block(data.view(), &s).unwrap();
5989 assert_eq!(block.group_ids[0], Some(1));
5990 assert_eq!(block.group_ids[1], Some(0));
5991 }
5992
5993 #[test]
5994 fn fixed_factor_at_fit_time_derives_vocabulary_and_never_false_rejects() {
5995 let mut s = fixed_factor_spec();
5999 s.frozen_levels = None;
6000 let data = array![[2000.0_f64], [2001.0], [2002.0], [2000.0]];
6001 let block = build_random_effect_block(data.view(), &s)
6002 .expect("fit-time build must not reject its own levels");
6003 assert_eq!(block.num_groups, 3);
6004 }
6005
6006 #[test]
6007 fn random_effect_still_tolerates_unseen_numeric_level() {
6008 let mut s = spec(); s.frozen_levels = Some(vec![2000.0_f64.to_bits(), 2001.0_f64.to_bits()]);
6013 let data = array![[2000.0_f64], [1999.0]];
6014 let block = build_random_effect_block(data.view(), &s)
6015 .expect("a random effect tolerates unseen levels");
6016 assert_eq!(block.group_ids[0], Some(0));
6017 assert_eq!(block.group_ids[1], None, "unseen level → population mean, not a reject");
6018 }
6019}
6020
6021impl SmoothDesign {
6022 pub fn map_term_coefficients(
6025 unconstrained: &Array1<f64>,
6026 shape: ShapeConstraint,
6027 ) -> Result<Array1<f64>, BasisError> {
6028 if unconstrained.is_empty() {
6029 crate::bail_invalid_basis!("unconstrained coefficient vector cannot be empty");
6030 }
6031 let mapped = match shape {
6032 ShapeConstraint::None => unconstrained.clone(),
6033 ShapeConstraint::MonotoneIncreasing => cumulative_exp(unconstrained, 1.0),
6034 ShapeConstraint::MonotoneDecreasing => cumulative_exp(unconstrained, -1.0),
6035 ShapeConstraint::Convex => second_cumulative_exp(unconstrained, 1.0),
6036 ShapeConstraint::Concave => second_cumulative_exp(unconstrained, -1.0),
6037 };
6038 Ok(mapped)
6039 }
6040}
6041
6042pub struct LocalSmoothTermBuild {
6043 pub dim: usize,
6044 pub design: DesignMatrix,
6045 pub penalties: Vec<Array2<f64>>,
6046 pub ops: Vec<Option<std::sync::Arc<dyn crate::analytic_penalties::PenaltyOp>>>,
6047 pub nullspaces: Vec<usize>,
6048 pub null_eigenvectors: Vec<Option<Array2<f64>>>,
6056 pub joint_null_rotation: Option<crate::basis::JointNullRotation>,
6063 pub penaltyinfo: Vec<PenaltyInfo>,
6064 pub pre_dropped_penaltyinfo: Vec<PenaltyInfo>,
6065 pub metadata: BasisMetadata,
6066 pub linear_constraints: Option<LinearInequalityConstraints>,
6067 pub box_reparam: bool,
6068 pub kronecker_factored: Option<KroneckerFactoredBasis>,
6069}
6070
6071#[derive(Clone)]
6072pub struct PcaScoresMemmapDesignOperator {
6073 mmap: Arc<memmap2::Mmap>,
6074 data_offset: usize,
6075 nrows: usize,
6076 ncols: usize,
6077 chunk_size: usize,
6078}
6079
6080impl PcaScoresMemmapDesignOperator {
6081 fn open(path: PathBuf, chunk_size: usize) -> Result<Self, BasisError> {
6082 let file = File::open(&path).map_err(|err| {
6083 BasisError::InvalidInput(format!(
6084 "failed to open lazy Pca .npy scores '{}': {err}",
6085 path.display()
6086 ))
6087 })?;
6088 let mmap = unsafe {
6094 memmap2::Mmap::map(&file).map_err(|err| {
6095 BasisError::InvalidInput(format!(
6096 "failed to memmap lazy Pca .npy scores '{}': {err}",
6097 path.display()
6098 ))
6099 })?
6100 };
6101 let (data_offset, nrows, ncols) = parse_f64_2d_npy_header(&mmap, &path)?;
6102 let expected = data_offset
6103 .checked_add(nrows.saturating_mul(ncols).saturating_mul(8))
6104 .ok_or_else(|| {
6105 BasisError::InvalidInput(format!(
6106 "lazy Pca .npy scores '{}' shape is too large",
6107 path.display()
6108 ))
6109 })?;
6110 if mmap.len() < expected {
6111 crate::bail_invalid_basis!(
6112 "lazy Pca .npy scores '{}' is truncated: header expects {} bytes, file has {}",
6113 path.display(),
6114 expected,
6115 mmap.len()
6116 );
6117 }
6118 Ok(Self {
6119 mmap: Arc::new(mmap),
6120 data_offset,
6121 nrows,
6122 ncols,
6123 chunk_size: chunk_size.max(1),
6124 })
6125 }
6126
6127 fn value(&self, row: usize, col: usize) -> f64 {
6128 let offset = self.data_offset + (row * self.ncols + col) * 8;
6129 let mut bytes = [0_u8; 8];
6130 bytes.copy_from_slice(&self.mmap[offset..offset + 8]);
6131 f64::from_le_bytes(bytes)
6132 }
6133
6134 fn chunk_rows(&self) -> usize {
6135 self.chunk_size.min(self.nrows.max(1))
6136 }
6137}
6138
6139impl LinearOperator for PcaScoresMemmapDesignOperator {
6140 fn nrows(&self) -> usize {
6141 self.nrows
6142 }
6143
6144 fn ncols(&self) -> usize {
6145 self.ncols
6146 }
6147
6148 fn apply(&self, vector: &Array1<f64>) -> Array1<f64> {
6149 assert_eq!(
6150 vector.len(),
6151 self.ncols,
6152 "lazy Pca apply vector length mismatch"
6153 );
6154 let mut out = Array1::<f64>::zeros(self.nrows);
6155 for start in (0..self.nrows).step_by(self.chunk_rows()) {
6156 let end = (start + self.chunk_rows()).min(self.nrows);
6157 for row in start..end {
6158 let mut acc = 0.0;
6159 for col in 0..self.ncols {
6160 acc += self.value(row, col) * vector[col];
6161 }
6162 out[row] = acc;
6163 }
6164 }
6165 out
6166 }
6167
6168 fn apply_transpose(&self, vector: &Array1<f64>) -> Array1<f64> {
6169 assert_eq!(
6170 vector.len(),
6171 self.nrows,
6172 "lazy Pca apply_transpose vector length mismatch"
6173 );
6174 let mut out = Array1::<f64>::zeros(self.ncols);
6175 for start in (0..self.nrows).step_by(self.chunk_rows()) {
6176 let end = (start + self.chunk_rows()).min(self.nrows);
6177 for row in start..end {
6178 let scale = vector[row];
6179 if scale == 0.0 {
6180 continue;
6181 }
6182 for col in 0..self.ncols {
6183 out[col] += scale * self.value(row, col);
6184 }
6185 }
6186 }
6187 out
6188 }
6189
6190 fn diag_xtw_x(&self, weights: &Array1<f64>) -> Result<Array2<f64>, String> {
6191 if weights.len() != self.nrows {
6192 return Err(format!(
6193 "lazy Pca diag_xtw_x weight length mismatch: weights={}, nrows={}",
6194 weights.len(),
6195 self.nrows
6196 ));
6197 }
6198 let mut gram = Array2::<f64>::zeros((self.ncols, self.ncols));
6199 for start in (0..self.nrows).step_by(self.chunk_rows()) {
6200 let end = (start + self.chunk_rows()).min(self.nrows);
6201 for row in start..end {
6202 let w = weights[row];
6203 if w == 0.0 {
6204 continue;
6205 }
6206 for a in 0..self.ncols {
6207 let xa = self.value(row, a);
6208 if xa == 0.0 {
6209 continue;
6210 }
6211 for b in a..self.ncols {
6212 gram[[a, b]] += w * xa * self.value(row, b);
6213 }
6214 }
6215 }
6216 }
6217 for a in 0..self.ncols {
6218 for b in 0..a {
6219 gram[[a, b]] = gram[[b, a]];
6220 }
6221 }
6222 Ok(gram)
6223 }
6224
6225 fn apply_weighted_normal(
6226 &self,
6227 weights: &Array1<f64>,
6228 vector: &Array1<f64>,
6229 penalty: Option<&Array2<f64>>,
6230 ridge: f64,
6231 ) -> Array1<f64> {
6232 assert_eq!(
6233 weights.len(),
6234 self.nrows,
6235 "lazy Pca weighted-normal weight mismatch"
6236 );
6237 assert_eq!(
6238 vector.len(),
6239 self.ncols,
6240 "lazy Pca weighted-normal vector mismatch"
6241 );
6242 let mut out = Array1::<f64>::zeros(self.ncols);
6243 for start in (0..self.nrows).step_by(self.chunk_rows()) {
6244 let end = (start + self.chunk_rows()).min(self.nrows);
6245 for row in start..end {
6246 let w = weights[row].max(0.0);
6247 if w == 0.0 {
6248 continue;
6249 }
6250 let mut row_dot = 0.0;
6251 for col in 0..self.ncols {
6252 row_dot += self.value(row, col) * vector[col];
6253 }
6254 if row_dot == 0.0 {
6255 continue;
6256 }
6257 let scaled = w * row_dot;
6258 for col in 0..self.ncols {
6259 out[col] += scaled * self.value(row, col);
6260 }
6261 }
6262 }
6263 if let Some(pen) = penalty {
6264 out += &pen.dot(vector);
6265 }
6266 if ridge > 0.0 {
6267 out += &vector.mapv(|x| ridge * x);
6268 }
6269 out
6270 }
6271}
6272
6273impl DenseDesignOperator for PcaScoresMemmapDesignOperator {
6274 fn compute_xtwy(&self, weights: &Array1<f64>, y: &Array1<f64>) -> Result<Array1<f64>, String> {
6275 if weights.len() != self.nrows || y.len() != self.nrows {
6276 return Err(format!(
6277 "lazy Pca compute_xtwy dimension mismatch: weights={}, y={}, nrows={}",
6278 weights.len(),
6279 y.len(),
6280 self.nrows
6281 ));
6282 }
6283 let mut out = Array1::<f64>::zeros(self.ncols);
6284 for start in (0..self.nrows).step_by(self.chunk_rows()) {
6285 let end = (start + self.chunk_rows()).min(self.nrows);
6286 for row in start..end {
6287 let scale = weights[row] * y[row];
6288 if scale == 0.0 {
6289 continue;
6290 }
6291 for col in 0..self.ncols {
6292 out[col] += scale * self.value(row, col);
6293 }
6294 }
6295 }
6296 Ok(out)
6297 }
6298
6299 fn row_chunk_into(
6300 &self,
6301 rows: Range<usize>,
6302 mut out: ArrayViewMut2<'_, f64>,
6303 ) -> Result<(), MatrixMaterializationError> {
6304 if rows.end > self.nrows || rows.start > rows.end {
6305 return Err(MatrixMaterializationError::MissingRowChunk {
6306 context: "lazy Pca row range out of bounds",
6307 });
6308 }
6309 if out.nrows() != rows.end - rows.start || out.ncols() != self.ncols {
6310 return Err(MatrixMaterializationError::MissingRowChunk {
6311 context: "lazy Pca row_chunk_into shape mismatch",
6312 });
6313 }
6314 for (local, row) in (rows.start..rows.end).enumerate() {
6315 for col in 0..self.ncols {
6316 out[[local, col]] = self.value(row, col);
6317 }
6318 }
6319 Ok(())
6320 }
6321
6322 fn to_dense(&self) -> Array2<f64> {
6323 let mut out = Array2::<f64>::zeros((self.nrows, self.ncols));
6324 self.row_chunk_into(0..self.nrows, out.view_mut())
6325 .expect("lazy Pca full materialization failed");
6326 out
6327 }
6328}
6329
6330pub fn parse_f64_2d_npy_header(
6331 bytes: &[u8],
6332 path: &PathBuf,
6333) -> Result<(usize, usize, usize), BasisError> {
6334 if bytes.len() < 10 || &bytes[0..6] != b"\x93NUMPY" {
6335 crate::bail_invalid_basis!("lazy Pca scores '{}' is not a .npy file", path.display());
6336 }
6337 let major = bytes[6];
6338 let header_len = match major {
6339 1 => u16::from_le_bytes([bytes[8], bytes[9]]) as usize,
6340 2 | 3 => {
6341 if bytes.len() < 12 {
6342 crate::bail_invalid_basis!(
6343 "lazy Pca scores '{}' has a truncated .npy header",
6344 path.display()
6345 );
6346 }
6347 u32::from_le_bytes([bytes[8], bytes[9], bytes[10], bytes[11]]) as usize
6348 }
6349 other => {
6350 crate::bail_invalid_basis!(
6351 "lazy Pca scores '{}' uses unsupported .npy version {}",
6352 path.display(),
6353 other
6354 );
6355 }
6356 };
6357 let header_start = if major == 1 { 10 } else { 12 };
6358 let data_offset = header_start + header_len;
6359 if bytes.len() < data_offset {
6360 crate::bail_invalid_basis!(
6361 "lazy Pca scores '{}' has a truncated .npy header",
6362 path.display()
6363 );
6364 }
6365 let header = std::str::from_utf8(&bytes[header_start..data_offset]).map_err(|err| {
6366 BasisError::InvalidInput(format!(
6367 "lazy Pca scores '{}' has a non-UTF8 .npy header: {err}",
6368 path.display()
6369 ))
6370 })?;
6371 if !(header.contains("'descr': '<f8'")
6372 || header.contains("\"descr\": \"<f8\"")
6373 || header.contains("'descr': '|f8'")
6374 || header.contains("\"descr\": \"|f8\""))
6375 {
6376 crate::bail_invalid_basis!(
6377 "lazy Pca scores '{}' must be float64 little-endian .npy",
6378 path.display()
6379 );
6380 }
6381 if header.contains("True") {
6382 crate::bail_invalid_basis!(
6383 "lazy Pca scores '{}' must be C-contiguous, not Fortran-ordered",
6384 path.display()
6385 );
6386 }
6387 let shape_pos = header.find("shape").ok_or_else(|| {
6388 BasisError::InvalidInput(format!(
6389 "lazy Pca scores '{}' .npy header is missing shape",
6390 path.display()
6391 ))
6392 })?;
6393 let open = header[shape_pos..].find('(').ok_or_else(|| {
6394 BasisError::InvalidInput(format!(
6395 "lazy Pca scores '{}' .npy header has malformed shape",
6396 path.display()
6397 ))
6398 })? + shape_pos;
6399 let close = header[open..].find(')').ok_or_else(|| {
6400 BasisError::InvalidInput(format!(
6401 "lazy Pca scores '{}' .npy header has malformed shape",
6402 path.display()
6403 ))
6404 })? + open;
6405 let dims = header[open + 1..close]
6406 .split(',')
6407 .map(str::trim)
6408 .filter(|part| !part.is_empty())
6409 .map(|part| part.parse::<usize>())
6410 .collect::<Result<Vec<_>, _>>()
6411 .map_err(|err| {
6412 BasisError::InvalidInput(format!(
6413 "lazy Pca scores '{}' .npy shape is not integral: {err}",
6414 path.display()
6415 ))
6416 })?;
6417 if dims.len() != 2 {
6418 crate::bail_invalid_basis!(
6419 "lazy Pca scores '{}' must have shape (N, K), got {:?}",
6420 path.display(),
6421 dims
6422 );
6423 }
6424 Ok((data_offset, dims[0], dims[1]))
6425}
6426
6427pub fn pca_center_mean(x: ArrayView2<'_, f64>) -> Result<Array1<f64>, BasisError> {
6428 if x.nrows() == 0 {
6429 crate::bail_invalid_basis!("Pca basis requires at least one row to compute center mean");
6430 }
6431 let mut mean = Array1::<f64>::zeros(x.ncols());
6432 for row in x.rows() {
6433 mean += &row;
6434 }
6435 mean.mapv_inplace(|v| v / x.nrows() as f64);
6436 Ok(mean)
6437}
6438
6439pub fn build_pca_smooth_basis(
6440 data: ArrayView2<'_, f64>,
6441 feature_cols: &[usize],
6442 basis_matrix: &Array2<f64>,
6443 centered: bool,
6444 smooth_penalty: f64,
6445 center_mean: Option<&Array1<f64>>,
6446 pca_basis_path: Option<&PathBuf>,
6447 chunk_size: usize,
6448) -> Result<BasisBuildResult, BasisError> {
6449 if let Some(path) = pca_basis_path {
6450 let op = PcaScoresMemmapDesignOperator::open(path.clone(), chunk_size)?;
6451 if op.nrows != data.nrows() {
6452 crate::bail_dim_basis!(
6453 "lazy Pca scores row mismatch: .npy has {}, data has {}",
6454 op.nrows,
6455 data.nrows()
6456 );
6457 }
6458 let k = op.ncols;
6459 let mut penalty = Array2::<f64>::eye(k);
6460 penalty.mapv_inplace(|v| v * smooth_penalty);
6461 let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
6462 filter_active_penalty_candidates_with_ops(vec![PenaltyCandidate {
6463 matrix: penalty,
6464 nullspace_dim_hint: 0,
6465 source: PenaltySource::Other("PcaRidge".to_string()),
6466 normalization_scale: 1.0,
6467 kronecker_factors: None,
6468 op: None,
6469 }])?;
6470 return Ok(BasisBuildResult {
6471 design: DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(Arc::new(op))),
6472 penalties,
6473 nullspace_dims,
6474 penaltyinfo,
6475 ops,
6476 null_eigenvectors,
6477 joint_null_rotation: None,
6478 metadata: BasisMetadata::Pca {
6479 feature_cols: feature_cols.to_vec(),
6480 basis_matrix: basis_matrix.clone(),
6481 centered,
6482 smooth_penalty,
6483 center_mean: center_mean.cloned(),
6484 pca_basis_path: Some(path.clone()),
6485 chunk_size: chunk_size.max(1),
6486 },
6487 kronecker_factored: None,
6488 });
6489 }
6490 if basis_matrix.nrows() != feature_cols.len() {
6491 crate::bail_dim_basis!(
6492 "Pca basis row mismatch: basis rows={}, feature columns={}",
6493 basis_matrix.nrows(),
6494 feature_cols.len()
6495 );
6496 }
6497 let mut x = select_columns(data, feature_cols)?;
6498 let mean = if centered {
6499 match center_mean {
6500 Some(mean) => mean.clone(),
6501 None => pca_center_mean(x.view())?,
6502 }
6503 } else {
6504 Array1::<f64>::zeros(feature_cols.len())
6505 };
6506 if centered {
6507 for mut row in x.rows_mut() {
6508 row -= &mean;
6509 }
6510 }
6511 let design = fast_ab(&x, basis_matrix);
6512 let k = basis_matrix.ncols();
6513 let mut penalty = Array2::<f64>::eye(k);
6514 penalty.mapv_inplace(|v| v * smooth_penalty);
6515 let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
6516 filter_active_penalty_candidates_with_ops(vec![PenaltyCandidate {
6517 matrix: penalty,
6518 nullspace_dim_hint: 0,
6519 source: PenaltySource::Other("PcaRidge".to_string()),
6520 normalization_scale: 1.0,
6521 kronecker_factors: None,
6522 op: None,
6523 }])?;
6524 Ok(BasisBuildResult {
6525 design: DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(design)),
6526 penalties,
6527 nullspace_dims,
6528 penaltyinfo,
6529 ops,
6530 null_eigenvectors,
6531 joint_null_rotation: None,
6532 metadata: BasisMetadata::Pca {
6533 feature_cols: feature_cols.to_vec(),
6534 basis_matrix: basis_matrix.clone(),
6535 centered,
6536 smooth_penalty,
6537 center_mean: centered.then_some(mean),
6538 pca_basis_path: None,
6539 chunk_size: chunk_size.max(1),
6540 },
6541 kronecker_factored: None,
6542 })
6543}
6544
6545pub fn defer_inner_model_centering_to_factor_level_wrapper(basis: &mut SmoothBasisSpec) {
6561 if let SmoothBasisSpec::BSpline1D { spec, .. } = basis
6562 && matches!(
6563 spec.identifiability,
6564 BSplineIdentifiability::WeightedSumToZero { .. }
6565 )
6566 {
6567 spec.identifiability = BSplineIdentifiability::None;
6568 }
6569}
6570
6571pub fn apply_by_variable_to_local_build(
6572 mut built: LocalSmoothTermBuild,
6573 data: ArrayView2<'_, f64>,
6574 by_col: usize,
6575 by: &ByVariableSpec,
6576 term_name: &str,
6577) -> Result<LocalSmoothTermBuild, BasisError> {
6578 if by_col >= data.ncols() {
6579 crate::bail_dim_basis!(
6580 "by-variable smooth term '{term_name}' references column {by_col}, but data has {} columns",
6581 data.ncols()
6582 );
6583 }
6584 let weights = match by {
6585 ByVariableSpec::Numeric => data.column(by_col).to_owned(),
6586 ByVariableSpec::Level { value_bits, .. } => {
6587 let value_bits = gam_data::canonical_level_bits(f64::from_bits(*value_bits));
6588 data.column(by_col).mapv(|value| {
6589 if gam_data::canonical_level_bits(value) == value_bits {
6590 1.0
6591 } else {
6592 0.0
6593 }
6594 })
6595 }
6596 };
6597 if weights.iter().any(|value| !value.is_finite()) {
6598 crate::bail_invalid_basis!(
6599 "by-variable smooth term '{term_name}' has non-finite by-column values"
6600 );
6601 }
6602
6603 let mut dense = built
6604 .design
6605 .try_to_dense_by_chunks("by-variable smooth row gating")
6606 .map_err(BasisError::InvalidInput)?;
6607 for (mut row, &weight) in dense.rows_mut().into_iter().zip(weights.iter()) {
6608 row.mapv_inplace(|value| value * weight);
6609 }
6610 built.design = DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(dense));
6611 built.kronecker_factored = None;
6612 Ok(built)
6613}
6614
6615pub fn build_by_smooth_local(
6626 data: ArrayView2<'_, f64>,
6627 term: &SmoothTermSpec,
6628 smooth: &SmoothBasisSpec,
6629 by_kind: &ByVarKind,
6630 workspace: &mut crate::basis::BasisWorkspace,
6631) -> Result<LocalSmoothTermBuild, BasisError> {
6632 let inner_term = SmoothTermSpec {
6633 name: term.name.clone(),
6634 basis: (*smooth).clone(),
6635 shape: term.shape,
6636 joint_null_rotation: None,
6637 };
6638 let inner = build_single_local_smooth_term(data, &inner_term, workspace)?;
6639
6640 match by_kind {
6641 ByVarKind::Numeric { feature_col } => {
6642 let inner_meta = inner.metadata.clone();
6643 let mut built = apply_by_variable_to_local_build(
6644 inner,
6645 data,
6646 *feature_col,
6647 &ByVariableSpec::Numeric,
6648 &term.name,
6649 )?;
6650 built.metadata = BasisMetadata::BySmooth {
6651 inner: Box::new(inner_meta),
6652 by_col: *feature_col,
6653 levels: None,
6654 ordered: false,
6655 };
6656 Ok(built)
6657 }
6658 ByVarKind::Factor {
6659 feature_col,
6660 frozen_levels,
6661 ordered,
6662 } => {
6663 let level_bits: Vec<u64> = if let Some(fl) = frozen_levels {
6666 fl.iter()
6667 .map(|&b| gam_data::canonical_level_bits(f64::from_bits(b)))
6668 .collect()
6669 } else {
6670 let col = data.column(*feature_col);
6671 let mut seen = BTreeSet::<u64>::new();
6672 for &v in col.iter() {
6673 if v.is_finite() {
6674 seen.insert(gam_data::canonical_level_bits(v));
6675 }
6676 }
6677 seen.into_iter().collect()
6678 };
6679 let n_levels = level_bits.len();
6680 if n_levels == 0 {
6681 crate::bail_invalid_basis!(
6682 "by-factor smooth term '{}': factor column {} has no observed levels",
6683 term.name,
6684 feature_col
6685 );
6686 }
6687 let p = inner.dim;
6688 let q = n_levels * p;
6689 let n = data.nrows();
6690
6691 let inner_dense = inner
6692 .design
6693 .try_to_dense_by_chunks("by-factor smooth design gating")
6694 .map_err(BasisError::InvalidInput)?;
6695
6696 let mut combined = Array2::<f64>::zeros((n, q));
6698 for (lvl_idx, &bits) in level_bits.iter().enumerate() {
6699 let col_start = lvl_idx * p;
6700 for row in 0..n {
6701 if gam_data::canonical_level_bits(data[[row, *feature_col]]) == bits {
6702 combined
6703 .slice_mut(s![row, col_start..col_start + p])
6704 .assign(&inner_dense.row(row));
6705 }
6706 }
6707 }
6708
6709 let inner_meta = inner.metadata.clone();
6721 let n_penalties = inner.penalties.len();
6722 let n_blocks = n_penalties.saturating_mul(n_levels);
6723 let mut penalties = Vec::<Array2<f64>>::with_capacity(n_blocks);
6724 let mut penaltyinfo = Vec::<PenaltyInfo>::with_capacity(n_blocks);
6725 let mut nullspaces = Vec::<usize>::with_capacity(n_blocks);
6726 for (pen_pos, s_inner) in inner.penalties.iter().enumerate() {
6727 for lvl in 0..n_levels {
6728 let off = lvl * p;
6729 let mut s_big = Array2::<f64>::zeros((q, q));
6730 s_big
6731 .slice_mut(s![off..off + p, off..off + p])
6732 .assign(s_inner);
6733 let (s_big, scale) = normalize_penalty_in_constrained_space(&s_big);
6734 let mut info = inner.penaltyinfo[pen_pos].clone();
6735 info.original_index = pen_pos * n_levels + lvl;
6738 info.normalization_scale *= scale;
6739 info.kronecker_factors = None;
6742 penalties.push(s_big);
6743 penaltyinfo.push(info);
6744 nullspaces.push(inner.nullspaces[pen_pos]);
6745 }
6746 }
6747
6748 let null_eigenvectors = vec![None; penalties.len()];
6749 let ops = vec![None; penalties.len()];
6750
6751 Ok(LocalSmoothTermBuild {
6752 dim: q,
6753 design: DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(combined)),
6754 penalties,
6755 ops,
6756 nullspaces,
6757 null_eigenvectors,
6758 joint_null_rotation: None,
6759 penaltyinfo,
6760 pre_dropped_penaltyinfo: inner.pre_dropped_penaltyinfo,
6761 metadata: BasisMetadata::BySmooth {
6762 inner: Box::new(inner_meta),
6763 by_col: *feature_col,
6764 levels: Some(level_bits),
6765 ordered: *ordered,
6766 },
6767 linear_constraints: None,
6768 box_reparam: false,
6769 kronecker_factored: None,
6770 })
6771 }
6772 }
6773}
6774
6775pub fn ensure_by_variable_specs_match(
6776 kind: &BySmoothKind,
6777 by: &ByVariableSpec,
6778 term_name: &str,
6779) -> Result<(), BasisError> {
6780 match (kind, by) {
6781 (BySmoothKind::Numeric, ByVariableSpec::Numeric) => Ok(()),
6782 (BySmoothKind::Level { level_bits }, ByVariableSpec::Level { value_bits, .. })
6783 if level_bits == value_bits =>
6784 {
6785 Ok(())
6786 }
6787 _ => Err(BasisError::InvalidInput(format!(
6788 "by-variable smooth term '{term_name}' has inconsistent by-variable specifications"
6789 ))),
6790 }
6791}
6792
6793pub fn build_factor_smooth(
6821 data: ArrayView2<'_, f64>,
6822 spec: &FactorSmoothSpec,
6823 term_name: &str,
6824 workspace: &mut crate::basis::BasisWorkspace,
6825) -> Result<LocalSmoothTermBuild, BasisError> {
6826 if spec.continuous_cols.len() != 1 {
6827 crate::bail_invalid_basis!(
6828 "factor smooth term '{}' currently supports exactly one continuous covariate; found {}",
6829 term_name,
6830 spec.continuous_cols.len()
6831 );
6832 }
6833 let feature_col = spec.continuous_cols[0];
6834 let group_col = spec.group_col;
6835 if feature_col >= data.ncols() || group_col >= data.ncols() {
6836 crate::bail_dim_basis!(
6837 "factor smooth term '{}' references columns ({}, {}) out of bounds for {} columns",
6838 term_name,
6839 feature_col,
6840 group_col,
6841 data.ncols()
6842 );
6843 }
6844
6845 if matches!(spec.flavour, FactorSmoothFlavour::Sz) {
6848 let levels = resolve_factor_smooth_levels(data, group_col, spec, term_name)?;
6849 let inner = SmoothBasisSpec::BSpline1D {
6850 feature_col,
6851 spec: factor_smooth_marginal_for_replay(&spec.marginal),
6852 };
6853 let sz_term = SmoothTermSpec {
6854 name: term_name.to_string(),
6855 basis: SmoothBasisSpec::FactorSumToZero {
6856 inner: Box::new(inner),
6857 by_col: group_col,
6858 levels: levels.clone(),
6859 frozen_global_orthogonality: None,
6860 },
6861 shape: ShapeConstraint::None,
6862 joint_null_rotation: None,
6863 };
6864 let mut built = build_single_local_smooth_term(data, &sz_term, workspace)?;
6865 let (knots, degree, periodic, marginal_is_cr) = match &built.metadata {
6886 BasisMetadata::BSpline1D {
6887 knots,
6888 periodic,
6889 degree,
6890 ..
6891 } => (
6892 knots.clone(),
6893 degree.unwrap_or(spec.marginal.degree),
6894 *periodic,
6895 false,
6896 ),
6897 BasisMetadata::CubicRegression1D { knots, .. } => {
6898 (knots.clone(), spec.marginal.degree, None, true)
6899 }
6900 other => {
6901 crate::bail_invalid_basis!(
6902 "sz factor smooth term '{}' produced an unexpected marginal metadata variant {:?}",
6903 term_name,
6904 other
6905 );
6906 }
6907 };
6908 built.metadata = BasisMetadata::FactorSmooth {
6909 continuous_cols: spec.continuous_cols.clone(),
6910 group_col,
6911 knots,
6912 degree,
6913 periodic,
6914 group_levels: levels,
6915 flavour: "sz".to_string(),
6916 marginal_is_cr,
6917 };
6918 return Ok(built);
6919 }
6920
6921 let levels = resolve_factor_smooth_levels(data, group_col, spec, term_name)?;
6922 let n_levels = levels.len();
6923 if n_levels < 2 {
6924 crate::bail_invalid_basis!(
6925 "factor smooth term '{}' requires at least two grouping levels; found {}",
6926 term_name,
6927 n_levels
6928 );
6929 }
6930
6931 let use_per_dim_null = matches!(
6939 &spec.flavour,
6940 FactorSmoothFlavour::Fs { m_null_penalty_orders }
6941 if m_null_penalty_orders.iter().copied().max().unwrap_or(0) >= 1
6942 );
6943
6944 let mut marginal_spec = factor_smooth_marginal_for_replay(&spec.marginal);
6950 if use_per_dim_null {
6951 marginal_spec.double_penalty = false;
6952 }
6953 let inner_term = SmoothTermSpec {
6954 name: format!("{term_name}::marginal"),
6955 basis: SmoothBasisSpec::BSpline1D {
6956 feature_col,
6957 spec: marginal_spec,
6958 },
6959 shape: ShapeConstraint::None,
6960 joint_null_rotation: None,
6961 };
6962 let inner = build_single_local_smooth_term(data, &inner_term, workspace)?;
6963 let mut base = inner
6964 .design
6965 .try_to_dense_by_chunks("factor smooth marginal")
6966 .map_err(BasisError::InvalidInput)?;
6967 if matches!(spec.flavour, FactorSmoothFlavour::Re) {
6968 let center = match &inner.metadata {
6978 BasisMetadata::BSpline1D { knots, .. } if !knots.is_empty() => {
6979 0.5 * (knots[0] + knots[knots.len() - 1])
6980 }
6981 _ => 0.0,
6982 };
6983 let mut linear = Array2::<f64>::ones((data.nrows(), 2));
6984 linear
6985 .column_mut(1)
6986 .assign(&data.column(feature_col).mapv(|x| x - center));
6987 base = linear;
6988 }
6989 let n = base.nrows();
6990 let p = base.ncols();
6991 let q = p * n_levels;
6992
6993 let mut dense = Array2::<f64>::zeros((n, q));
6996 for i in 0..n {
6997 let bits = gam_data::canonical_level_bits(data[[i, group_col]]);
6998 let level_idx = levels.iter().position(|b| *b == bits).ok_or_else(|| {
6999 BasisError::InvalidInput(format!(
7000 "factor smooth term '{term_name}' saw an unseen grouping level at row {}",
7001 i + 1
7002 ))
7003 })?;
7004 let start = level_idx * p;
7005 dense
7006 .slice_mut(s![i, start..start + p])
7007 .assign(&base.row(i));
7008 }
7009
7010 let marginal_penalties: Vec<Array2<f64>> = if matches!(spec.flavour, FactorSmoothFlavour::Re) {
7016 (0..p)
7017 .map(|j| {
7018 let mut s = Array2::<f64>::zeros((p, p));
7019 s[[j, j]] = 1.0;
7020 s
7021 })
7022 .collect()
7023 } else {
7024 inner.penalties.clone()
7025 };
7026 let marginal_penaltyinfo: Vec<PenaltyInfo> = if matches!(spec.flavour, FactorSmoothFlavour::Re)
7027 {
7028 (0..p)
7029 .map(|j| PenaltyInfo {
7030 source: PenaltySource::Primary,
7031 original_index: j,
7032 active: true,
7033 effective_rank: 1,
7034 dropped_reason: None,
7035 nullspace_dim_hint: p.saturating_sub(1),
7036 normalization_scale: 1.0,
7037 kronecker_factors: None,
7038 })
7039 .collect()
7040 } else {
7041 inner.penaltyinfo.clone()
7042 };
7043 if marginal_penalties.len() != marginal_penaltyinfo.len() {
7044 crate::bail_invalid_basis!(
7045 "internal factor-smooth penalty metadata mismatch for term '{}': penalties={}, infos={}",
7046 term_name,
7047 marginal_penalties.len(),
7048 marginal_penaltyinfo.len()
7049 );
7050 }
7051
7052 let mut penalties = Vec::<Array2<f64>>::with_capacity(marginal_penalties.len());
7053 let mut penaltyinfo = Vec::<PenaltyInfo>::with_capacity(marginal_penalties.len());
7054 for (penalty_pos, s_inner) in marginal_penalties.iter().enumerate() {
7055 let mut s_big = Array2::<f64>::zeros((q, q));
7056 for level in 0..n_levels {
7057 let start = level * p;
7058 s_big
7059 .slice_mut(s![start..start + p, start..start + p])
7060 .assign(s_inner);
7061 }
7062 let (s_big, factor_smooth_scale) = normalize_penalty_in_constrained_space(&s_big);
7063 let mut info = marginal_penaltyinfo[penalty_pos].clone();
7064 info.original_index = penalty_pos;
7065 info.normalization_scale *= factor_smooth_scale;
7066 info.nullspace_dim_hint = info.nullspace_dim_hint.saturating_mul(n_levels);
7067 info.kronecker_factors = None;
7068 penalties.push(s_big);
7069 penaltyinfo.push(info);
7070 }
7071
7072 let mut nullspaces: Vec<usize> = if matches!(spec.flavour, FactorSmoothFlavour::Re) {
7073 vec![q.saturating_sub(n_levels); p]
7074 } else {
7075 inner
7076 .nullspaces
7077 .iter()
7078 .map(|ns| ns.saturating_mul(n_levels))
7079 .collect()
7080 };
7081
7082 if use_per_dim_null
7112 && let Some(Some(z)) = inner.null_eigenvectors.first()
7113 && z.nrows() == p
7114 {
7115 for k in 0..z.ncols() {
7116 let zk = z.column(k);
7121 let mut p_k = Array2::<f64>::zeros((p, p));
7122 for a in 0..p {
7123 for b in 0..p {
7124 p_k[[a, b]] = zk[a] * zk[b];
7125 }
7126 }
7127 let mut s_null = Array2::<f64>::zeros((q, q));
7128 for level in 0..n_levels {
7129 let start = level * p;
7130 s_null
7131 .slice_mut(s![start..start + p, start..start + p])
7132 .assign(&p_k);
7133 }
7134 let (s_null, null_scale) = normalize_penalty_in_constrained_space(&s_null);
7135 let null_block = crate::basis::analyze_penalty_block_with_op(&s_null, None)?;
7136 if null_block.rank > 0 {
7137 let original_index = penalties.len();
7138 penalties.push(null_block.sym_penalty);
7139 nullspaces.push(null_block.nullity);
7140 penaltyinfo.push(PenaltyInfo {
7141 source: PenaltySource::Primary,
7142 original_index,
7143 active: true,
7144 effective_rank: null_block.rank,
7145 dropped_reason: None,
7146 nullspace_dim_hint: null_block.nullity,
7147 normalization_scale: null_scale,
7148 kronecker_factors: None,
7149 });
7150 }
7151 }
7152 }
7153 let null_eigenvectors = crate::basis::recompute_null_eigenvectors(&penalties)?;
7154 let joint_null_rotation = crate::basis::compute_joint_null_rotation(&penalties)?;
7155
7156 let (knots, degree, periodic) = match &inner.metadata {
7159 BasisMetadata::BSpline1D {
7160 knots,
7161 periodic,
7162 degree,
7163 ..
7164 } => (
7165 knots.clone(),
7166 degree.unwrap_or(spec.marginal.degree),
7167 *periodic,
7168 ),
7169 other => {
7170 crate::bail_invalid_basis!(
7171 "factor smooth term '{}' produced an unexpected marginal metadata variant {:?}",
7172 term_name,
7173 other
7174 );
7175 }
7176 };
7177 let flavour_tag = match &spec.flavour {
7178 FactorSmoothFlavour::Fs { .. } => "fs",
7179 FactorSmoothFlavour::Sz => "sz",
7180 FactorSmoothFlavour::Re => "re",
7181 }
7182 .to_string();
7183 let metadata = BasisMetadata::FactorSmooth {
7184 continuous_cols: spec.continuous_cols.clone(),
7185 group_col,
7186 knots,
7187 degree,
7188 periodic,
7189 group_levels: levels,
7190 flavour: flavour_tag,
7191 marginal_is_cr: false,
7194 };
7195
7196 let ops = vec![None; penalties.len()];
7197 Ok(LocalSmoothTermBuild {
7198 dim: q,
7199 design: DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(dense)),
7200 penalties,
7201 ops,
7202 nullspaces,
7203 null_eigenvectors,
7204 joint_null_rotation,
7205 penaltyinfo,
7206 pre_dropped_penaltyinfo: Vec::new(),
7207 metadata,
7208 linear_constraints: None,
7209 box_reparam: false,
7210 kronecker_factored: None,
7211 })
7212}
7213
7214pub fn resolve_factor_smooth_levels(
7218 data: ArrayView2<'_, f64>,
7219 group_col: usize,
7220 spec: &FactorSmoothSpec,
7221 term_name: &str,
7222) -> Result<Vec<u64>, BasisError> {
7223 if let Some(frozen) = &spec.group_frozen_levels {
7224 if frozen.is_empty() {
7225 crate::bail_invalid_basis!(
7226 "factor smooth term '{}' has an empty frozen level list",
7227 term_name
7228 );
7229 }
7230 return Ok(frozen
7231 .iter()
7232 .map(|&b| gam_data::canonical_level_bits(f64::from_bits(b)))
7233 .collect());
7234 }
7235 let mut bits: Vec<u64> = data
7236 .column(group_col)
7237 .iter()
7238 .map(|v| gam_data::canonical_level_bits(*v))
7239 .collect();
7240 bits.sort_by(|a, b| {
7241 f64::from_bits(*a)
7242 .partial_cmp(&f64::from_bits(*b))
7243 .unwrap_or(std::cmp::Ordering::Equal)
7244 });
7245 bits.dedup();
7246 Ok(bits)
7247}
7248
7249pub fn factor_smooth_marginal_for_replay(marginal: &BSplineBasisSpec) -> BSplineBasisSpec {
7256 let mut m = marginal.clone();
7257 m.identifiability = BSplineIdentifiability::None;
7258 m
7259}
7260
7261pub fn build_single_local_smooth_term(
7262 data: ArrayView2<'_, f64>,
7263 term: &SmoothTermSpec,
7264 workspace: &mut crate::basis::BasisWorkspace,
7265) -> Result<LocalSmoothTermBuild, BasisError> {
7266 if term.shape != ShapeConstraint::None && !shape_supports_basis(term) {
7267 crate::bail_invalid_basis!(
7268 "ShapeConstraint::{:?} is unsupported for term '{}'",
7269 term.shape,
7270 term.name
7271 );
7272 }
7273 if let SmoothBasisSpec::ByVariable {
7274 inner,
7275 by_col,
7276 kind,
7277 by,
7278 } = &term.basis
7279 {
7280 ensure_by_variable_specs_match(kind, by, &term.name)?;
7281 let mut inner_basis = (**inner).clone();
7282 if matches!(by, ByVariableSpec::Level { .. }) {
7289 defer_inner_model_centering_to_factor_level_wrapper(&mut inner_basis);
7290 }
7291 let inner_term = SmoothTermSpec {
7292 name: term.name.clone(),
7293 basis: inner_basis,
7294 shape: term.shape,
7295 joint_null_rotation: None,
7296 };
7297 let built = build_single_local_smooth_term(data, &inner_term, workspace)?;
7298 return apply_by_variable_to_local_build(built, data, *by_col, by, &term.name);
7299 }
7300
7301 if let SmoothBasisSpec::BySmooth { smooth, by_kind } = &term.basis {
7304 return build_by_smooth_local(data, term, smooth, by_kind, workspace);
7305 }
7306
7307 let mut shape_axis_col: Option<usize> = None;
7308 let mut built: BasisBuildResult = match &term.basis {
7309 SmoothBasisSpec::FactorSumToZero {
7310 inner,
7311 by_col,
7312 levels,
7313 ..
7314 } => {
7315 if *by_col >= data.ncols() {
7316 crate::bail_dim_basis!(
7317 "term '{}' by column {} out of bounds for {} columns",
7318 term.name,
7319 by_col,
7320 data.ncols()
7321 );
7322 }
7323 if levels.len() < 2 {
7324 crate::bail_invalid_basis!(
7325 "sum-to-zero factor smooth term '{}' requires at least two levels",
7326 term.name
7327 );
7328 }
7329 if term.shape != ShapeConstraint::None {
7330 crate::bail_invalid_basis!(
7331 "ShapeConstraint::{:?} is unsupported for sum-to-zero factor smooth term '{}'",
7332 term.shape,
7333 term.name
7334 );
7335 }
7336 let inner_term = SmoothTermSpec {
7337 name: format!("{}::inner", term.name),
7338 basis: (**inner).clone(),
7339 shape: ShapeConstraint::None,
7340 joint_null_rotation: None,
7341 };
7342 let mut inner_built = build_single_local_smooth_term(data, &inner_term, workspace)?;
7343 let inner_null_eigenvectors = inner_built.null_eigenvectors.clone();
7347 let base = inner_built
7348 .design
7349 .try_to_dense_by_chunks("sum-to-zero factor smooth")
7350 .map_err(BasisError::InvalidInput)?;
7351 let n = base.nrows();
7352 let p = base.ncols();
7353 let l_minus_one = levels.len() - 1;
7354 let canon_levels: Vec<u64> = levels
7357 .iter()
7358 .map(|&b| gam_data::canonical_level_bits(f64::from_bits(b)))
7359 .collect();
7360 let mut dense = Array2::<f64>::zeros((n, p * l_minus_one));
7361 for i in 0..n {
7362 let bits = gam_data::canonical_level_bits(data[[i, *by_col]]);
7363 let level_idx = canon_levels.iter().position(|b| *b == bits).ok_or_else(|| {
7364 BasisError::InvalidInput(format!(
7365 "sum-to-zero factor smooth term '{}' saw an unseen level at row {}",
7366 term.name,
7367 i + 1
7368 ))
7369 })?;
7370 if level_idx < l_minus_one {
7371 let start = level_idx * p;
7372 dense
7373 .slice_mut(s![i, start..start + p])
7374 .assign(&base.row(i));
7375 } else {
7376 for level in 0..l_minus_one {
7377 let start = level * p;
7378 dense
7379 .slice_mut(s![i, start..start + p])
7380 .assign(&base.row(i).mapv(|v| -v));
7381 }
7382 }
7383 }
7384 let mut penalties = Vec::<Array2<f64>>::with_capacity(inner_built.penalties.len());
7385 let active_penalty_indices = inner_built
7386 .penaltyinfo
7387 .iter()
7388 .enumerate()
7389 .filter_map(|(idx, info)| info.active.then_some(idx))
7390 .collect::<Vec<_>>();
7391 if active_penalty_indices.len() != inner_built.penalties.len() {
7392 crate::bail_invalid_basis!(
7393 "internal sz penalty metadata mismatch: activeinfos={}, penalties={}",
7394 active_penalty_indices.len(),
7395 inner_built.penalties.len()
7396 );
7397 }
7398 let stz_per_group_penalty = |s_inner: &Array2<f64>, which_level: usize| -> Array2<f64> {
7433 let mut s_big = Array2::<f64>::zeros((p * l_minus_one, p * l_minus_one));
7434 if which_level < l_minus_one {
7435 let k = which_level;
7437 let mut block = s_big.slice_mut(s![k * p..(k + 1) * p, k * p..(k + 1) * p]);
7438 block.assign(s_inner);
7439 } else {
7440 for a in 0..l_minus_one {
7442 for b in 0..l_minus_one {
7443 let mut block =
7444 s_big.slice_mut(s![a * p..(a + 1) * p, b * p..(b + 1) * p]);
7445 block.assign(s_inner);
7446 }
7447 }
7448 }
7449 s_big
7450 };
7451 let mut nullspaces = Vec::<usize>::with_capacity(penalties.capacity());
7457 for (penalty_pos, s_inner) in inner_built.penalties.iter().enumerate() {
7458 let info_idx = active_penalty_indices[penalty_pos];
7459 let base_info = inner_built.penaltyinfo[info_idx].clone();
7460 let marginal_nullity = inner_built.nullspaces.get(penalty_pos).copied().unwrap_or(0);
7461 for which_level in 0..=l_minus_one {
7463 let raw = stz_per_group_penalty(s_inner, which_level);
7464 let (s_big, group_scale) = normalize_penalty_in_constrained_space(&raw);
7465 let block = crate::basis::analyze_penalty_block_with_op(&s_big, None)?;
7466 if block.rank == 0 {
7467 continue;
7468 }
7469 if which_level == 0 {
7470 inner_built.penaltyinfo[info_idx].normalization_scale *= group_scale;
7473 inner_built.penaltyinfo[info_idx].original_index = penalties.len();
7474 inner_built.penaltyinfo[info_idx].effective_rank = block.rank;
7475 inner_built.penaltyinfo[info_idx].nullspace_dim_hint = block.nullity;
7476 } else {
7477 let mut info = base_info.clone();
7478 info.original_index = penalties.len();
7479 info.normalization_scale = base_info.normalization_scale * group_scale;
7480 info.effective_rank = block.rank;
7481 info.nullspace_dim_hint = block.nullity;
7482 info.kronecker_factors = None;
7483 inner_built.penaltyinfo.push(info);
7484 }
7485 penalties.push(block.sym_penalty);
7486 nullspaces.push(marginal_nullity);
7492 }
7493 }
7494
7495 if let Some(Some(z)) = inner_null_eigenvectors.first()
7513 && z.nrows() == p
7514 {
7515 for k in 0..z.ncols() {
7516 let zk = z.column(k);
7517 let mut p_k = Array2::<f64>::zeros((p, p));
7518 for a in 0..p {
7519 for b in 0..p {
7520 p_k[[a, b]] = zk[a] * zk[b];
7521 }
7522 }
7523 let stz_pooled_null = {
7528 let mut s_big = Array2::<f64>::zeros((p * l_minus_one, p * l_minus_one));
7529 for a in 0..l_minus_one {
7530 for b in 0..l_minus_one {
7531 let factor = if a == b { 2.0 } else { 1.0 };
7532 let mut block =
7533 s_big.slice_mut(s![a * p..(a + 1) * p, b * p..(b + 1) * p]);
7534 block.assign(&p_k.mapv(|v| v * factor));
7535 }
7536 }
7537 s_big
7538 };
7539 let (s_null, null_scale) =
7540 normalize_penalty_in_constrained_space(&stz_pooled_null);
7541 let null_block = crate::basis::analyze_penalty_block_with_op(&s_null, None)?;
7542 if null_block.rank > 0 {
7543 let original_index = penalties.len();
7544 penalties.push(null_block.sym_penalty);
7545 nullspaces.push(null_block.nullity);
7546 inner_built.penaltyinfo.push(PenaltyInfo {
7547 source: PenaltySource::Primary,
7548 original_index,
7549 active: true,
7550 effective_rank: null_block.rank,
7551 dropped_reason: None,
7552 nullspace_dim_hint: null_block.nullity,
7553 normalization_scale: null_scale,
7554 kronecker_factors: None,
7555 });
7556 }
7557 }
7558 }
7559 inner_built.dim = p * l_minus_one;
7560 inner_built.design = DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(dense));
7561 inner_built.penalties = penalties;
7562 inner_built.ops = vec![None; inner_built.penalties.len()];
7563 inner_built.nullspaces = nullspaces;
7564 inner_built.null_eigenvectors =
7571 crate::basis::recompute_null_eigenvectors(&inner_built.penalties)?;
7572 inner_built.joint_null_rotation =
7573 crate::basis::compute_joint_null_rotation(&inner_built.penalties)?;
7574 inner_built.kronecker_factored = None;
7575 return Ok(inner_built);
7576 }
7577 SmoothBasisSpec::BSpline1D { feature_col, spec } => {
7578 if *feature_col >= data.ncols() {
7579 crate::bail_dim_basis!(
7580 "term '{}' feature column {} out of bounds for {} columns",
7581 term.name,
7582 feature_col,
7583 data.ncols()
7584 );
7585 }
7586 let mut spec_local = spec.clone();
7587 if term.shape != ShapeConstraint::None {
7588 spec_local.identifiability = BSplineIdentifiability::None;
7591 }
7592 build_bspline_basis_1d(data.column(*feature_col), &spec_local)?
7596 }
7597 SmoothBasisSpec::ThinPlate {
7598 feature_cols,
7599 spec,
7600 input_scales,
7601 } => {
7602 if term.shape != ShapeConstraint::None {
7603 if feature_cols.len() != 1 {
7604 crate::bail_invalid_basis!(
7605 "ShapeConstraint::{:?} for term '{}' on ThinPlate basis requires exactly 1 feature axis; found {}",
7606 term.shape,
7607 term.name,
7608 feature_cols.len()
7609 );
7610 }
7611 shape_axis_col = Some(feature_cols[0]);
7612 }
7613 let mut x = select_columns(data, feature_cols)?;
7614 let (scales, length_scale_eff) = if let Some(s) = input_scales {
7620 apply_input_standardization(&mut x, s);
7621 (
7622 Some(s.clone()),
7623 compensate_length_scale_for_standardization(spec.length_scale, s),
7624 )
7625 } else if let Some(s) = compute_spatial_input_scales(x.view()) {
7626 apply_input_standardization(&mut x, &s);
7627 let l_eff = compensate_length_scale_for_standardization(spec.length_scale, &s);
7628 (Some(s), l_eff)
7629 } else {
7630 (None, spec.length_scale)
7631 };
7632 let mut spec_local = spec.clone();
7633 spec_local.length_scale = length_scale_eff;
7634 if matches!(
7635 spec_local.identifiability,
7636 SpatialIdentifiability::OrthogonalToParametric
7637 ) {
7638 spec_local.identifiability = SpatialIdentifiability::None;
7639 }
7640 let mut result = build_thin_plate_basis(x.view(), &spec_local).map_err(|err| {
7641 rewrite_thin_plate_knots_error(err, &term.name, feature_cols.len(), spec)
7642 })?;
7643 match &mut result.metadata {
7651 BasisMetadata::ThinPlate {
7652 input_scales: ms,
7653 length_scale,
7654 ..
7655 } => {
7656 *ms = scales;
7657 *length_scale = spec.length_scale;
7658 }
7659 BasisMetadata::Duchon {
7660 input_scales: ms,
7661 length_scale,
7662 ..
7663 } => {
7664 if let (Some(s), Some(realized)) = (scales.as_ref(), *length_scale) {
7689 let inv_sigma_geom =
7690 compensate_length_scale_for_standardization(1.0, s);
7691 if inv_sigma_geom.is_finite() && inv_sigma_geom > 0.0 {
7692 *length_scale = Some(realized / inv_sigma_geom);
7693 }
7694 }
7695 *ms = scales;
7696 }
7697 _ => {}
7698 }
7699 result
7700 }
7701 SmoothBasisSpec::Sphere { feature_cols, spec } => {
7702 if term.shape != ShapeConstraint::None {
7703 crate::bail_invalid_basis!(
7704 "ShapeConstraint::{:?} for term '{}' is not supported on spherical splines",
7705 term.shape,
7706 term.name
7707 );
7708 }
7709 let x = select_columns(data, feature_cols)?;
7710 build_spherical_spline_basis(x.view(), spec)?
7711 }
7712 SmoothBasisSpec::ConstantCurvature { feature_cols, spec } => {
7713 if term.shape != ShapeConstraint::None {
7714 crate::bail_invalid_basis!(
7715 "ShapeConstraint::{:?} for term '{}' is not supported on constant-curvature smooths",
7716 term.shape,
7717 term.name
7718 );
7719 }
7720 let x = select_columns(data, feature_cols)?;
7727 build_constant_curvature_basis(x.view(), spec)?
7728 }
7729 SmoothBasisSpec::MeasureJet {
7730 feature_cols,
7731 spec,
7732 input_scales,
7733 } => {
7734 if term.shape != ShapeConstraint::None {
7735 crate::bail_invalid_basis!(
7736 "ShapeConstraint::{:?} for term '{}' is not supported on measure-jet smooths",
7737 term.shape,
7738 term.name
7739 );
7740 }
7741 let mut x = select_columns(data, feature_cols)?;
7742 let (scales, length_scale_eff) = if let Some(s) = input_scales {
7754 apply_input_standardization(&mut x, s);
7755 (Some(s.clone()), spec.length_scale)
7756 } else if let Some(s) = compute_spatial_input_scales(x.view()) {
7757 apply_input_standardization(&mut x, &s);
7758 let l_eff = if spec.length_scale > 0.0 {
7759 compensate_length_scale_for_standardization(spec.length_scale, &s)
7760 } else {
7761 spec.length_scale
7762 };
7763 (Some(s), l_eff)
7764 } else {
7765 (None, spec.length_scale)
7766 };
7767 let mut spec_local = spec.clone();
7768 spec_local.length_scale = length_scale_eff;
7769 let mut result = build_measure_jet_basis(x.view(), &spec_local)?;
7770 if let BasisMetadata::MeasureJet {
7771 input_scales: ms, ..
7772 } = &mut result.metadata
7773 {
7774 *ms = scales;
7775 }
7776 result
7777 }
7778 SmoothBasisSpec::Matern {
7779 feature_cols,
7780 spec,
7781 input_scales,
7782 } => {
7783 if term.shape != ShapeConstraint::None {
7784 if feature_cols.len() != 1 {
7785 crate::bail_invalid_basis!(
7786 "ShapeConstraint::{:?} for term '{}' on Matern basis requires exactly 1 feature axis; found {}",
7787 term.shape,
7788 term.name,
7789 feature_cols.len()
7790 );
7791 }
7792 shape_axis_col = Some(feature_cols[0]);
7793 }
7794 let mut x = select_columns(data, feature_cols)?;
7795 let (scales, length_scale_eff) = if let Some(s) = input_scales {
7810 apply_input_standardization(&mut x, s);
7811 (
7812 Some(s.clone()),
7813 compensate_length_scale_for_standardization(spec.length_scale, s),
7814 )
7815 } else if let Some(s) = compute_spatial_input_scales(x.view()) {
7816 apply_input_standardization(&mut x, &s);
7817 let l_eff = compensate_length_scale_for_standardization(spec.length_scale, &s);
7818 (Some(s), l_eff)
7819 } else {
7820 (None, spec.length_scale)
7821 };
7822 let mut spec_local = spec.clone();
7823 spec_local.length_scale = length_scale_eff;
7824 let mut result = build_matern_basiswithworkspace(x.view(), &spec_local, workspace)?;
7825 if let BasisMetadata::Matern {
7826 input_scales,
7827 length_scale,
7828 ..
7829 } = &mut result.metadata
7830 {
7831 *input_scales = scales;
7832 *length_scale = spec.length_scale;
7833 }
7834 result
7835 }
7836 SmoothBasisSpec::Duchon {
7837 feature_cols,
7838 spec,
7839 input_scales,
7840 } => {
7841 if term.shape != ShapeConstraint::None {
7842 if feature_cols.len() != 1 {
7843 crate::bail_invalid_basis!(
7844 "ShapeConstraint::{:?} for term '{}' on Duchon basis requires exactly 1 feature axis; found {}",
7845 term.shape,
7846 term.name,
7847 feature_cols.len()
7848 );
7849 }
7850 shape_axis_col = Some(feature_cols[0]);
7851 }
7852 let mut x = select_columns(data, feature_cols)?;
7853 let (scales, length_scale_eff) = if let Some(s) = input_scales {
7864 apply_input_standardization(&mut x, s);
7865 (
7866 Some(s.clone()),
7867 compensate_optional_length_scale_for_standardization(spec.length_scale, s),
7868 )
7869 } else if let Some(s) = compute_spatial_input_scales(x.view()) {
7870 apply_input_standardization(&mut x, &s);
7871 let l_eff =
7872 compensate_optional_length_scale_for_standardization(spec.length_scale, &s);
7873 (Some(s), l_eff)
7874 } else {
7875 (None, spec.length_scale)
7876 };
7877 let mut spec_local = spec.clone();
7878 spec_local.length_scale = length_scale_eff;
7879 if let (Some(s), crate::basis::OneDimensionalBoundary::Cyclic { start, end }) =
7890 (scales.as_ref(), spec_local.boundary.clone())
7891 && s.len() == 1
7892 && s[0] > 0.0
7893 {
7894 spec_local.boundary = crate::basis::OneDimensionalBoundary::Cyclic {
7895 start: start / s[0],
7896 end: end / s[0],
7897 };
7898 }
7899 if matches!(
7900 spec_local.identifiability,
7901 SpatialIdentifiability::OrthogonalToParametric
7902 ) {
7903 spec_local.identifiability = SpatialIdentifiability::None;
7904 }
7905 let mut result = build_duchon_basiswithworkspace(x.view(), &spec_local, workspace)?;
7906 if let BasisMetadata::Duchon {
7907 input_scales,
7908 length_scale,
7909 ..
7910 } = &mut result.metadata
7911 {
7912 *input_scales = scales;
7913 *length_scale = spec.length_scale;
7914 }
7915 result
7916 }
7917 SmoothBasisSpec::Pca {
7918 feature_cols,
7919 basis_matrix,
7920 centered,
7921 smooth_penalty,
7922 center_mean,
7923 pca_basis_path,
7924 chunk_size,
7925 } => {
7926 if term.shape != ShapeConstraint::None {
7927 crate::bail_invalid_basis!(
7928 "ShapeConstraint::{:?} for term '{}' is not supported on Pca basis",
7929 term.shape,
7930 term.name
7931 );
7932 }
7933 build_pca_smooth_basis(
7934 data,
7935 feature_cols,
7936 basis_matrix,
7937 *centered,
7938 *smooth_penalty,
7939 center_mean.as_ref(),
7940 pca_basis_path.as_ref(),
7941 *chunk_size,
7942 )?
7943 }
7944 SmoothBasisSpec::TensorBSpline { feature_cols, spec } => {
7945 build_tensor_bspline_basis(data, feature_cols, spec)?
7946 }
7947 SmoothBasisSpec::ByVariable { .. } => {
7948 crate::bail_invalid_basis!(
7949 "internal: ByVariable smooths must return before inner basis dispatch"
7950 );
7951 }
7952 SmoothBasisSpec::BySmooth { .. } => {
7953 crate::bail_invalid_basis!("internal: BySmooth smooths must be lowered to ByVariable before inner basis dispatch"
7954 .to_string(),);
7955 }
7956 SmoothBasisSpec::FactorSmooth { spec } => {
7957 if term.shape != ShapeConstraint::None {
7958 crate::bail_invalid_basis!(
7959 "ShapeConstraint::{:?} is unsupported for factor smooth term '{}'",
7960 term.shape,
7961 term.name
7962 );
7963 }
7964 return build_factor_smooth(data, spec, &term.name, workspace);
7965 }
7966 };
7967
7968 if let SmoothBasisSpec::Matern { .. } = &term.basis {
7984 let (penalties, nullspace_dims, penaltyinfo) =
7985 matern_operator_penalty_triplet_from_metadata(&built.metadata)?;
7986 built.penalties = penalties;
7987 built.nullspace_dims = nullspace_dims;
7988 built.penaltyinfo = penaltyinfo;
7989 }
7990
7991 let p_local = built.design.ncols();
7992 let mut metadata = built.metadata.clone();
7993 let kron_factored = if term.shape == ShapeConstraint::None {
7996 built.kronecker_factored
7997 } else {
7998 None
7999 };
8000 let mut design_t = built.design;
8001 let mut penalties_t: Vec<Array2<f64>> = built.penalties;
8002 let mut ops_t: Vec<Option<std::sync::Arc<dyn crate::analytic_penalties::PenaltyOp>>> =
8007 built.ops;
8008 if matches!(
8009 spatial_identifiability_policy(term),
8010 Some(SpatialIdentifiability::OrthogonalToParametric)
8011 ) {
8012 metadata = freeze_raw_spatial_metadata(metadata, design_t.ncols());
8013 }
8014
8015 let active_penaltyinfo_t = built
8016 .penaltyinfo
8017 .iter()
8018 .filter(|info| info.active)
8019 .cloned()
8020 .collect::<Vec<_>>();
8021 let pre_dropped_penaltyinfo_t = built
8022 .penaltyinfo
8023 .iter()
8024 .filter(|info| !info.active)
8025 .cloned()
8026 .collect::<Vec<_>>();
8027 let use_box_reparam =
8028 term.shape != ShapeConstraint::None && shape_uses_box_reparameterization(&term.basis);
8029 if let Some((order, sign)) = shape_order_and_sign(term.shape)
8030 && use_box_reparam
8031 {
8032 let t = if order == 2 {
8047 let bspline_meta = match &metadata {
8048 BasisMetadata::BSpline1D {
8049 knots,
8050 degree,
8051 periodic,
8052 ..
8053 } if periodic.is_none() => Some((knots.clone(), degree.unwrap_or(0))),
8054 _ => None,
8055 };
8056 match bspline_meta {
8057 Some((knots, degree)) if degree >= 1 => {
8058 let greville = crate::basis::compute_greville_abscissae(&knots, degree)?;
8059 if greville.len() != p_local {
8060 crate::bail_invalid_basis!(
8061 "shape-constraint Greville abscissae count {} does not match basis dim {} for term '{}'",
8062 greville.len(),
8063 p_local,
8064 term.name
8065 );
8066 }
8067 convex_divided_difference_transform_matrix(&greville, sign)?
8068 }
8069 _ => cumulative_sum_transform_matrix(p_local, order, sign),
8070 }
8071 } else {
8072 cumulative_sum_transform_matrix(p_local, order, sign)
8073 };
8074 let inner_dense = match design_t {
8078 DesignMatrix::Dense(d) => d,
8079 DesignMatrix::Sparse(sp) => gam_linalg::matrix::DenseDesignMatrix::from(
8080 sp.try_to_dense_arc("shape-constrained coefficient transform")
8081 .map_err(BasisError::InvalidInput)?,
8082 ),
8083 };
8084 let coeff_op = gam_linalg::matrix::CoefficientTransformOperator::new(inner_dense, t.clone())
8085 .map_err(|e| BasisError::InvalidInput(format!("CoefficientTransformOperator: {e}")))?;
8086 design_t = DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(Arc::new(coeff_op)));
8087 if penalties_t.len() != active_penaltyinfo_t.len() {
8088 crate::bail_invalid_basis!(
8089 "internal box-reparam penalty/info mismatch for term '{}': penalties={}, infos={}",
8090 term.name,
8091 penalties_t.len(),
8092 active_penaltyinfo_t.len()
8093 );
8094 }
8095 let transformed_wiggliness = penalties_t
8111 .iter()
8112 .zip(active_penaltyinfo_t.iter())
8113 .find(|(_, info)| !matches!(info.source, PenaltySource::DoublePenaltyNullspace))
8114 .map(|(s_local, _)| {
8115 let tt_s = fast_atb(&t, s_local);
8116 fast_ab(&tt_s, &t)
8117 });
8118 let mut rebuilt = Vec::with_capacity(penalties_t.len());
8119 for (s_local, info) in penalties_t.iter().zip(active_penaltyinfo_t.iter()) {
8120 if matches!(info.source, PenaltySource::DoublePenaltyNullspace) {
8121 if order == 2 {
8156 let tt_s = fast_atb(&t, s_local);
8157 rebuilt.push(fast_ab(&tt_s, &t));
8158 } else {
8159 let s_wiggle_t = transformed_wiggliness.as_ref().ok_or_else(|| {
8160 BasisError::InvalidInput(format!(
8161 "box-reparam term '{}' has a double-penalty ridge but no primary wiggliness penalty to derive its nullspace from",
8162 term.name
8163 ))
8164 })?;
8165 let ridge = crate::basis::build_nullspace_shrinkage_penalty(s_wiggle_t)?
8166 .map(|shrink| shrink.sym_penalty)
8167 .unwrap_or_else(|| Array2::<f64>::zeros((p_local, p_local)));
8168 rebuilt.push(ridge);
8169 }
8170 } else {
8171 let tt_s = fast_atb(&t, s_local);
8172 rebuilt.push(fast_ab(&tt_s, &t));
8173 }
8174 }
8175 penalties_t = rebuilt;
8176 ops_t = vec![None; penalties_t.len()];
8179 }
8180 if penalties_t.len() != active_penaltyinfo_t.len() {
8181 crate::bail_invalid_basis!(
8182 "internal penalty metadata mismatch for term '{}': active penalties={}, active infos={}",
8183 term.name,
8184 penalties_t.len(),
8185 active_penaltyinfo_t.len()
8186 );
8187 }
8188 if ops_t.len() != penalties_t.len() {
8189 ops_t = vec![None; penalties_t.len()];
8190 }
8191 let penalty_candidates = penalties_t
8192 .into_iter()
8193 .zip(active_penaltyinfo_t.into_iter())
8194 .zip(ops_t.into_iter())
8195 .map(
8196 |((matrix, info), op_in)| -> Result<PenaltyCandidate, BasisError> {
8197 let (matrix, c_new) = normalize_penalty_in_constrained_space(&matrix);
8198 let normalization_scale = info.normalization_scale * c_new;
8199 let op_scale = 1.0 / c_new;
8200 let kronecker_scale = 1.0 / c_new;
8201 let scaled_op = if op_scale > 0.0 && op_scale.is_finite() {
8204 op_in.map(|op| {
8205 std::sync::Arc::new(crate::analytic_penalties::ScaledPenaltyOp::new(
8206 op, op_scale,
8207 ))
8208 as std::sync::Arc<dyn crate::analytic_penalties::PenaltyOp>
8209 })
8210 } else {
8211 None
8212 };
8213 let kronecker_factors = info.kronecker_factors.map(|mut factors| {
8214 if let Some(first) = factors.first_mut() {
8215 first.mapv_inplace(|v| v * kronecker_scale);
8216 }
8217 factors
8218 });
8219 Ok(PenaltyCandidate {
8220 nullspace_dim_hint: info.nullspace_dim_hint,
8221 matrix,
8222 source: info.source,
8223 normalization_scale,
8224 kronecker_factors,
8225 op: scaled_op,
8226 })
8227 },
8228 )
8229 .collect::<Result<Vec<_>, _>>()?;
8230 let (penalties_t, nullspaces_t, penaltyinfo_t, null_eigenvectors_t, ops_t) =
8231 crate::basis::filter_active_penalty_candidates_with_ops(penalty_candidates)?;
8232 let shape_linear_constraints = if term.shape != ShapeConstraint::None && !use_box_reparam {
8233 let axis = shape_axis_col.ok_or_else(|| {
8234 BasisError::InvalidInput(format!(
8235 "internal shape-constraint axis missing for term '{}'",
8236 term.name
8237 ))
8238 })?;
8239 let (x_shape_eval, design_shape_eval) =
8240 build_shape_constraint_design_1d(data, term, &metadata, axis)?;
8241 build_shape_linear_constraints_1d(
8242 x_shape_eval.view(),
8243 design_shape_eval.view(),
8244 term.shape,
8245 )?
8246 } else {
8247 None
8248 };
8249 let linear_constraints_local = merge_linear_constraints_global(shape_linear_constraints, None);
8250
8251 let joint_null_rotation = match term.joint_null_rotation.clone() {
8270 Some(persisted) => Some(persisted),
8271 None if smooth_has_frozen_identifiability(term) => None,
8272 None if kron_factored.is_some() => None,
8273 None => crate::basis::compute_joint_null_rotation(&penalties_t)?,
8274 };
8275
8276 Ok(LocalSmoothTermBuild {
8277 dim: p_local,
8278 design: design_t,
8279 penalties: penalties_t,
8280 ops: ops_t,
8281 nullspaces: nullspaces_t,
8282 null_eigenvectors: null_eigenvectors_t,
8283 joint_null_rotation,
8284 penaltyinfo: penaltyinfo_t,
8285 pre_dropped_penaltyinfo: pre_dropped_penaltyinfo_t,
8286 metadata,
8287 linear_constraints: linear_constraints_local,
8288 box_reparam: use_box_reparam,
8289 kronecker_factored: kron_factored,
8290 })
8291}
8292
8293pub fn build_smooth_design(
8294 data: ArrayView2<'_, f64>,
8295 terms: &[SmoothTermSpec],
8296) -> Result<RawSmoothDesign, BasisError> {
8297 let mut ws = crate::basis::BasisWorkspace::new();
8298 build_smooth_design_withworkspace(data, terms, &mut ws)
8299}
8300
8301pub fn build_smooth_design_withworkspace(
8308 data: ArrayView2<'_, f64>,
8309 terms: &[SmoothTermSpec],
8310 workspace: &mut crate::basis::BasisWorkspace,
8311) -> Result<RawSmoothDesign, BasisError> {
8312 validate_smooth_terms_finite_inputs(data, terms)?;
8313 build_smooth_design_withworkspace_unvalidated(data, terms, workspace)
8314}
8315
8316pub fn build_smooth_design_withworkspace_unvalidated(
8317 data: ArrayView2<'_, f64>,
8318 terms: &[SmoothTermSpec],
8319 workspace: &mut crate::basis::BasisWorkspace,
8320) -> Result<RawSmoothDesign, BasisError> {
8321 let mut planned_blocks = plan_joint_spatial_centers_for_term_blocks(data, &[terms.to_vec()])?;
8322 let planned_terms = planned_blocks.pop().ok_or_else(|| {
8323 BasisError::InvalidInput(
8324 "joint spatial center planner returned no smooth blocks".to_string(),
8325 )
8326 })?;
8327 let policy = workspace.policy().clone();
8328 let local_builds: Vec<LocalSmoothTermBuild> = {
8329 use rayon::iter::{IntoParallelIterator, ParallelIterator};
8330 planned_terms
8331 .into_par_iter()
8332 .map(|term| {
8333 let mut term_workspace = crate::basis::BasisWorkspace::with_policy(policy.clone());
8334 build_single_local_smooth_term(data, &term, &mut term_workspace)
8335 })
8336 .collect::<Result<Vec<_>, _>>()?
8337 };
8338
8339 let total_p: usize = local_builds.iter().map(|built| built.dim).sum();
8340
8341 let mut local_designs: Vec<DesignMatrix> = Vec::with_capacity(local_builds.len());
8342 let mut terms_out = Vec::<SmoothTerm>::with_capacity(terms.len());
8343 let mut penalties_global = Vec::<BlockwisePenalty>::new();
8344 let mut nullspace_dims_global = Vec::<usize>::new();
8345 let mut penaltyinfo_global = Vec::<PenaltyBlockInfo>::new();
8346 let mut dropped_penaltyinfo_global = Vec::<DroppedPenaltyBlockInfo>::new();
8347 let mut coefficient_lower_bounds = Array1::<f64>::from_elem(total_p, f64::NEG_INFINITY);
8348 let mut any_bounds = false;
8349 let mut linear_constraintsrows: Vec<(usize, usize, Array1<f64>)> = Vec::new();
8354 let mut linear_constraints_b: Vec<f64> = Vec::new();
8355
8356 let mut col_start = 0usize;
8357 for (term, mut built) in terms.iter().zip(local_builds.into_iter()) {
8358 let p_local = built.dim;
8359 let col_end = col_start + p_local;
8360 let lb_local = if built.box_reparam {
8361 shape_lower_bounds_local(term.shape, p_local)
8362 } else {
8363 None
8364 };
8365
8366 let applied_rotation: Option<crate::basis::JointNullRotation> = match (
8398 built.joint_null_rotation.take(),
8399 lb_local.is_some(),
8400 built.linear_constraints.is_some(),
8401 ) {
8402 (Some(rot), false, false) => {
8403 let q = &rot.rotation;
8404 let dense = built
8405 .design
8406 .try_to_dense_by_chunks("joint-null absorption rotation")
8407 .map_err(BasisError::InvalidInput)?;
8408 let rotated = gam_linalg::faer_ndarray::fast_ab(&dense, q);
8409 built.design = DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(rotated));
8410 built.penalties = built
8411 .penalties
8412 .into_iter()
8413 .map(|s_local| {
8414 let qt_s = gam_linalg::faer_ndarray::fast_atb(q, &s_local);
8415 gam_linalg::faer_ndarray::fast_ab(&qt_s, q)
8416 })
8417 .collect();
8418 built.ops = vec![None; built.penalties.len()];
8419 built.kronecker_factored = None;
8420 Some(rot)
8421 }
8422 (Some(_), _, _) => None,
8423 (None, _, _) => None,
8424 };
8425
8426 let activeinfos = built
8427 .penaltyinfo
8428 .iter()
8429 .filter(|info| info.active)
8430 .collect::<Vec<_>>();
8431 if activeinfos.len() != built.penalties.len() {
8432 crate::bail_invalid_basis!(
8433 "internal penalty info mismatch for term '{}': activeinfos={}, penalties={}",
8434 term.name,
8435 activeinfos.len(),
8436 built.penalties.len()
8437 );
8438 }
8439 for (((s_local, &ns), info), op_local) in built
8440 .penalties
8441 .iter()
8442 .zip(built.nullspaces.iter())
8443 .zip(activeinfos.into_iter())
8444 .zip(built.ops.iter())
8445 {
8446 let global_index = penalties_global.len();
8447 penalties_global.push(
8448 BlockwisePenalty::new(col_start..col_end, s_local.clone())
8449 .with_op(op_local.clone()),
8450 );
8451 nullspace_dims_global.push(ns);
8452 let mut penalty = info.clone();
8453 penalty.nullspace_dim_hint = ns;
8454 penaltyinfo_global.push(PenaltyBlockInfo {
8455 global_index,
8456 termname: Some(term.name.clone()),
8457 penalty,
8458 });
8459 }
8460 for info in built.penaltyinfo.iter().filter(|info| !info.active) {
8461 dropped_penaltyinfo_global.push(DroppedPenaltyBlockInfo {
8462 termname: Some(term.name.clone()),
8463 penalty: info.clone(),
8464 });
8465 }
8466 for info in &built.pre_dropped_penaltyinfo {
8467 dropped_penaltyinfo_global.push(DroppedPenaltyBlockInfo {
8468 termname: Some(term.name.clone()),
8469 penalty: info.clone(),
8470 });
8471 }
8472
8473 if let Some(lin_local) = &built.linear_constraints {
8474 for r in 0..lin_local.a.nrows() {
8475 linear_constraintsrows.push((col_start, col_end, lin_local.a.row(r).to_owned()));
8476 linear_constraints_b.push(lin_local.b[r]);
8477 }
8478 }
8479 if let Some(lb_local) = &lb_local {
8480 coefficient_lower_bounds
8481 .slice_mut(s![col_start..col_end])
8482 .assign(lb_local);
8483 any_bounds = true;
8484 }
8485
8486 local_designs.push(built.design);
8488
8489 terms_out.push(SmoothTerm {
8490 name: term.name.clone(),
8491 coeff_range: col_start..col_end,
8492 shape: term.shape,
8493 penalties_local: built.penalties,
8494 nullspace_dims: built.nullspaces,
8495 penaltyinfo_local: built.penaltyinfo,
8496 metadata: built.metadata,
8497 lower_bounds_local: lb_local,
8498 linear_constraints_local: built.linear_constraints,
8499 kronecker_factored: built.kronecker_factored.take(),
8500 joint_null_rotation: applied_rotation,
8501 unabsorbed_global_orthogonality: None,
8502 });
8503
8504 col_start = col_end;
8505 }
8506
8507 assert_eq!(
8508 penalties_global.len(),
8509 nullspace_dims_global.len(),
8510 "global smooth penalty/nullspace bookkeeping diverged"
8511 );
8512 assert_eq!(
8513 penalties_global.len(),
8514 penaltyinfo_global.len(),
8515 "global smooth penalty metadata bookkeeping diverged"
8516 );
8517
8518 Ok(RawSmoothDesign {
8519 term_designs: local_designs,
8520 penalties: penalties_global,
8521 nullspace_dims: nullspace_dims_global,
8522 penaltyinfo: penaltyinfo_global,
8523 dropped_penaltyinfo: dropped_penaltyinfo_global,
8524 terms: terms_out,
8525 coefficient_lower_bounds: if any_bounds {
8526 Some(coefficient_lower_bounds)
8527 } else {
8528 None
8529 },
8530 linear_constraints: if linear_constraintsrows.is_empty() {
8531 None
8532 } else {
8533 let mut a = Array2::<f64>::zeros((linear_constraintsrows.len(), total_p));
8534 for (i, (cs, ce, values)) in linear_constraintsrows.iter().enumerate() {
8535 a.row_mut(i).slice_mut(s![*cs..*ce]).assign(values);
8536 }
8537 Some(LinearInequalityConstraints {
8538 a,
8539 b: Array1::from_vec(linear_constraints_b),
8540 })
8541 },
8542 })
8543}