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)]
1091 pub categorical_levels: Vec<(usize, u64)>,
1092 #[serde(default = "default_linear_term_double_penalty")]
1099 pub double_penalty: bool,
1100 #[serde(default)]
1101 pub coefficient_geometry: LinearCoefficientGeometry,
1102 #[serde(default)]
1103 pub coefficient_min: Option<f64>,
1104 #[serde(default)]
1105 pub coefficient_max: Option<f64>,
1106}
1107
1108impl LinearTermSpec {
1109 pub fn effective_feature_cols(&self) -> Vec<usize> {
1112 if self.feature_cols.is_empty() {
1113 vec![self.feature_col]
1114 } else {
1115 self.feature_cols.clone()
1116 }
1117 }
1118
1119 pub fn is_interaction(&self) -> bool {
1121 self.feature_cols.len() > 1 || !self.categorical_levels.is_empty()
1122 }
1123
1124 pub fn realized_design_column(&self, data: ArrayView2<'_, f64>) -> Result<Array1<f64>, String> {
1136 let n = data.nrows();
1137 let p = data.ncols();
1138 let bounds = |col: usize| -> Result<(), String> {
1139 if col >= p {
1140 Err(format!(
1141 "linear term '{}' feature column {} out of bounds for {} columns",
1142 self.name, col, p
1143 ))
1144 } else {
1145 Ok(())
1146 }
1147 };
1148
1149 let mut column = if self.categorical_levels.is_empty() {
1154 let cols = self.effective_feature_cols();
1155 for &c in &cols {
1156 bounds(c)?;
1157 }
1158 let mut acc = data.column(cols[0]).to_owned();
1159 for &c in cols.iter().skip(1) {
1160 acc *= &data.column(c);
1161 }
1162 acc
1163 } else {
1164 let mut acc = Array1::<f64>::ones(n);
1165 for &c in &self.feature_cols {
1166 bounds(c)?;
1167 acc *= &data.column(c);
1168 }
1169 acc
1170 };
1171
1172 for &(col, level_bits) in &self.categorical_levels {
1173 bounds(col)?;
1174 let gate = data.column(col);
1175 for (out, &v) in column.iter_mut().zip(gate.iter()) {
1176 if v.to_bits() != level_bits {
1177 *out = 0.0;
1178 }
1179 }
1180 }
1181
1182 Ok(column)
1183 }
1184}
1185
1186pub const fn default_linear_term_double_penalty() -> bool {
1187 false
1195}
1196
1197pub const fn default_pca_smooth_penalty() -> f64 {
1198 1.0
1199}
1200
1201pub const fn default_pca_chunk_size() -> usize {
1202 4096
1203}
1204
1205#[derive(Debug, Clone, Serialize, Deserialize)]
1211pub struct RandomEffectTermSpec {
1212 pub name: String,
1213 pub feature_col: usize,
1214 pub drop_first_level: bool,
1217 #[serde(default = "default_random_effect_penalized")]
1221 pub penalized: bool,
1222 #[serde(default)]
1225 pub frozen_levels: Option<Vec<u64>>,
1226 #[serde(default = "default_random_effect_lenient_unseen")]
1239 pub lenient_unseen: bool,
1240}
1241
1242pub fn default_random_effect_penalized() -> bool {
1243 true
1244}
1245
1246pub fn default_random_effect_lenient_unseen() -> bool {
1247 true
1248}
1249
1250pub fn validate_measure_jet_positive_vec_len(
1251 label: &str,
1252 term_name: &str,
1253 field: &str,
1254 values: &[f64],
1255 expected: usize,
1256) -> Result<(), String> {
1257 if values.len() != expected {
1258 return Err(SmoothError::invalid_config(format!(
1259 "{label} term '{term_name}' frozen MeasureJet {field} has length {}, expected {expected}",
1260 values.len()
1261 ))
1262 .into());
1263 }
1264 if values
1265 .iter()
1266 .any(|value| !(value.is_finite() && *value > 0.0))
1267 {
1268 return Err(SmoothError::invalid_config(format!(
1269 "{label} term '{term_name}' frozen MeasureJet {field} values must be positive and finite"
1270 ))
1271 .into());
1272 }
1273 Ok(())
1274}
1275
1276#[derive(Debug, Clone, Serialize, Deserialize)]
1277pub struct TermCollectionSpec {
1278 pub linear_terms: Vec<LinearTermSpec>,
1279 pub random_effect_terms: Vec<RandomEffectTermSpec>,
1280 pub smooth_terms: Vec<SmoothTermSpec>,
1281}
1282
1283pub fn validate_smooth_basis_frozen(
1284 basis: &SmoothBasisSpec,
1285 label: &str,
1286 term_name: &str,
1287) -> Result<(), String> {
1288 match basis {
1289 SmoothBasisSpec::ByVariable { inner, .. }
1290 | SmoothBasisSpec::FactorSumToZero { inner, .. } => {
1291 validate_smooth_basis_frozen(inner, label, term_name)
1292 }
1293 SmoothBasisSpec::BSpline1D { spec, .. } => {
1294 if !matches!(
1295 spec.knotspec,
1296 BSplineKnotSpec::Provided(_)
1297 | BSplineKnotSpec::PeriodicUniform { .. }
1298 | BSplineKnotSpec::NaturalCubicRegression { .. }
1299 ) {
1300 return Err(format!(
1301 "{label} term '{term_name}' is not frozen: BSpline knotspec must be Provided, PeriodicUniform, or NaturalCubicRegression"
1302 ));
1303 }
1304 Ok(())
1305 }
1306 SmoothBasisSpec::ThinPlate { spec, .. } => {
1307 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1308 return Err(format!(
1309 "{label} term '{term_name}' is not frozen: ThinPlate centers must be UserProvided"
1310 ));
1311 }
1312 if matches!(
1313 spec.identifiability,
1314 SpatialIdentifiability::OrthogonalToParametric
1315 ) {
1316 return Err(format!(
1317 "{label} term '{term_name}' is not frozen: ThinPlate identifiability must be FrozenTransform or None"
1318 ));
1319 }
1320 Ok(())
1321 }
1322 _ => Ok(()),
1323 }
1324}
1325
1326impl TermCollectionSpec {
1327 pub fn write_structural_shape_hash(&self, h: &mut gam_runtime::warm_start::Fingerprinter) {
1341 h.write_str("term-collection");
1342 h.write_usize(self.linear_terms.len());
1343 for linear in &self.linear_terms {
1344 h.write_str(&linear.name);
1345 }
1346 h.write_usize(self.random_effect_terms.len());
1347 h.write_usize(self.smooth_terms.len());
1348 for smooth in &self.smooth_terms {
1349 h.write_str(&smooth.name);
1350 h.write_str(smooth.basis.structural_kind());
1351 for col in smooth.basis.structural_feature_cols() {
1352 h.write_usize(col);
1353 }
1354 }
1355 }
1356
1357 pub fn validate_frozen(&self, label: &str) -> Result<(), String> {
1361 for linear in &self.linear_terms {
1362 if let (Some(min), Some(max)) = (linear.coefficient_min, linear.coefficient_max)
1363 && (!min.is_finite() || !max.is_finite() || min > max)
1364 {
1365 return Err(SmoothError::invalid_config(format!(
1366 "{label} linear term '{}' has invalid coefficient constraint [{min}, {max}]",
1367 linear.name
1368 ))
1369 .into());
1370 }
1371 if let Some(min) = linear.coefficient_min
1372 && !min.is_finite()
1373 {
1374 return Err(SmoothError::invalid_config(format!(
1375 "{label} linear term '{}' has non-finite coefficient minimum {min}",
1376 linear.name
1377 ))
1378 .into());
1379 }
1380 if let Some(max) = linear.coefficient_max
1381 && !max.is_finite()
1382 {
1383 return Err(SmoothError::invalid_config(format!(
1384 "{label} linear term '{}' has non-finite coefficient maximum {max}",
1385 linear.name
1386 ))
1387 .into());
1388 }
1389 if let LinearCoefficientGeometry::Bounded { min, max, prior } =
1390 &linear.coefficient_geometry
1391 {
1392 if !min.is_finite() || !max.is_finite() || min >= max {
1393 return Err(SmoothError::invalid_config(format!(
1394 "{label} bounded term '{}' has invalid bounds [{min}, {max}]",
1395 linear.name
1396 ))
1397 .into());
1398 }
1399 match prior {
1400 BoundedCoefficientPriorSpec::None | BoundedCoefficientPriorSpec::Uniform => {}
1401 BoundedCoefficientPriorSpec::Beta { a, b } => {
1402 if !a.is_finite() || !b.is_finite() || *a < 1.0 || *b < 1.0 {
1403 return Err(SmoothError::invalid_config(format!(
1404 "{label} bounded term '{}' has invalid Beta prior ({a}, {b})",
1405 linear.name
1406 ))
1407 .into());
1408 }
1409 }
1410 }
1411 }
1412 }
1413 for st in &self.smooth_terms {
1414 match &st.basis {
1415 SmoothBasisSpec::ByVariable { inner, .. } => {
1416 validate_smooth_basis_frozen(inner, label, &st.name)?;
1417 let nested = SmoothTermSpec {
1418 name: st.name.clone(),
1419 basis: (**inner).clone(),
1420 shape: st.shape,
1421 joint_null_rotation: None,
1422 };
1423 TermCollectionSpec {
1424 linear_terms: Vec::new(),
1425 random_effect_terms: Vec::new(),
1426 smooth_terms: vec![nested],
1427 }
1428 .validate_frozen(label)?;
1429 }
1430 SmoothBasisSpec::FactorSumToZero { inner, levels, .. } => {
1431 if levels.len() < 2 {
1432 return Err(format!(
1433 "{label} term '{}' has invalid frozen sz levels",
1434 st.name
1435 ));
1436 }
1437 validate_smooth_basis_frozen(inner, label, &st.name)?;
1438 }
1439 SmoothBasisSpec::BSpline1D { spec, .. } => {
1440 if !matches!(
1441 spec.knotspec,
1442 BSplineKnotSpec::Provided(_)
1443 | BSplineKnotSpec::PeriodicUniform { .. }
1444 | BSplineKnotSpec::NaturalCubicRegression { .. }
1445 ) {
1446 return Err(SmoothError::invalid_config(format!(
1447 "{label} term '{}' is not frozen: BSpline knotspec must be Provided, PeriodicUniform, or NaturalCubicRegression",
1448 st.name
1449 ))
1450 .into());
1451 }
1452 }
1453 SmoothBasisSpec::ThinPlate { spec, .. } => {
1454 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1455 return Err(SmoothError::invalid_config(format!(
1456 "{label} term '{}' is not frozen: ThinPlate centers must be UserProvided",
1457 st.name
1458 ))
1459 .into());
1460 }
1461 if matches!(
1462 spec.identifiability,
1463 SpatialIdentifiability::OrthogonalToParametric
1464 ) {
1465 return Err(SmoothError::invalid_config(format!(
1466 "{label} term '{}' is not frozen: ThinPlate identifiability must be FrozenTransform or None",
1467 st.name
1468 ))
1469 .into());
1470 }
1471 }
1472 SmoothBasisSpec::Sphere { spec, .. } => {
1473 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1474 return Err(SmoothError::invalid_config(format!(
1475 "{label} term '{}' is not frozen: Sphere centers must be UserProvided",
1476 st.name
1477 ))
1478 .into());
1479 }
1480 if matches!(spec.method, crate::basis::SphereMethod::Harmonic)
1481 && spec.max_degree.is_none_or(|d| d == 0)
1482 {
1483 return Err(format!(
1484 "{label} term '{}' is not frozen: sphere max_degree must be positive",
1485 st.name
1486 ));
1487 }
1488 }
1489 SmoothBasisSpec::ConstantCurvature { spec, .. } => {
1490 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1491 return Err(SmoothError::invalid_config(format!(
1492 "{label} term '{}' is not frozen: ConstantCurvature centers must be UserProvided",
1493 st.name
1494 ))
1495 .into());
1496 }
1497 if !(spec.length_scale.is_finite() && spec.length_scale > 0.0) {
1498 return Err(SmoothError::invalid_config(format!(
1499 "{label} term '{}' is not frozen: ConstantCurvature length_scale must be the realized positive value",
1500 st.name
1501 ))
1502 .into());
1503 }
1504 }
1505 SmoothBasisSpec::MeasureJet { spec, .. } => {
1506 let centers = match &spec.center_strategy {
1507 CenterStrategy::UserProvided(centers) => centers,
1508 _ => {
1509 return Err(SmoothError::invalid_config(format!(
1510 "{label} term '{}' is not frozen: MeasureJet centers must be UserProvided",
1511 st.name
1512 ))
1513 .into());
1514 }
1515 };
1516 if centers.nrows() == 0 {
1517 return Err(SmoothError::invalid_config(format!(
1518 "{label} term '{}' is not frozen: MeasureJet centers are empty",
1519 st.name
1520 ))
1521 .into());
1522 }
1523 if !(spec.length_scale.is_finite() && spec.length_scale > 0.0) {
1524 return Err(SmoothError::invalid_config(format!(
1525 "{label} term '{}' is not frozen: MeasureJet length_scale must be the realized positive value",
1526 st.name
1527 ))
1528 .into());
1529 }
1530 let frozen = spec.frozen_quadrature.as_ref().ok_or_else(|| {
1533 SmoothError::invalid_config(format!(
1534 "{label} term '{}' is not frozen: MeasureJet frozen_quadrature payload is missing",
1535 st.name
1536 ))
1537 })?;
1538 if frozen.masses.len() != centers.nrows() {
1539 return Err(SmoothError::invalid_config(format!(
1540 "{label} term '{}' frozen MeasureJet has {} masses for {} centers",
1541 st.name,
1542 frozen.masses.len(),
1543 centers.nrows()
1544 ))
1545 .into());
1546 }
1547 let total_mass = frozen.masses.sum();
1548 if frozen
1549 .masses
1550 .iter()
1551 .any(|mass| !(mass.is_finite() && *mass >= 0.0))
1552 || !(total_mass.is_finite() && total_mass > 0.0)
1553 {
1554 return Err(SmoothError::invalid_config(format!(
1555 "{label} term '{}' frozen MeasureJet masses must be finite, nonnegative, and have positive total mass",
1556 st.name
1557 ))
1558 .into());
1559 }
1560 let n_levels = frozen.eps_band.len();
1561 if n_levels == 0
1562 || frozen
1563 .eps_band
1564 .iter()
1565 .any(|eps| !(eps.is_finite() && *eps > 0.0))
1566 {
1567 return Err(SmoothError::invalid_config(format!(
1568 "{label} term '{}' frozen MeasureJet eps_band must be nonempty, finite, and positive",
1569 st.name
1570 ))
1571 .into());
1572 }
1573 for (idx, pair) in frozen.eps_band.windows(2).enumerate() {
1574 if pair[1] <= pair[0] {
1575 return Err(SmoothError::invalid_config(format!(
1576 "{label} term '{}' frozen MeasureJet eps_band is not strictly ascending at {idx}: {} then {}",
1577 st.name,
1578 pair[0],
1579 pair[1]
1580 ))
1581 .into());
1582 }
1583 }
1584 validate_measure_jet_positive_vec_len(
1585 label,
1586 &st.name,
1587 "support_means",
1588 &frozen.support_means,
1589 n_levels,
1590 )?;
1591 let per_level = crate::basis::measure_jet_multiscale_mode(spec);
1599 if per_level {
1600 validate_measure_jet_positive_vec_len(
1601 label,
1602 &st.name,
1603 "penalty_normalization_scales",
1604 &frozen.penalty_normalization_scales,
1605 n_levels,
1606 )?;
1607 validate_measure_jet_positive_vec_len(
1608 label,
1609 &st.name,
1610 "raw_penalty_normalization_scales",
1611 &frozen.raw_penalty_normalization_scales,
1612 n_levels,
1613 )?;
1614 if frozen.fused_penalty_normalization_scale.is_some() {
1615 return Err(SmoothError::invalid_config(format!(
1616 "{label} term '{}' per-level MeasureJet must not carry a fused penalty normalization scale",
1617 st.name
1618 ))
1619 .into());
1620 }
1621 } else {
1622 if !frozen.penalty_normalization_scales.is_empty()
1623 || !frozen.raw_penalty_normalization_scales.is_empty()
1624 {
1625 return Err(SmoothError::invalid_config(format!(
1626 "{label} term '{}' fused MeasureJet must not carry per-level penalty normalization scales",
1627 st.name
1628 ))
1629 .into());
1630 }
1631 match frozen.fused_penalty_normalization_scale {
1632 Some(scale) if scale.is_finite() && scale > 0.0 => {}
1633 Some(scale) => {
1634 return Err(SmoothError::invalid_config(format!(
1635 "{label} term '{}' fused MeasureJet penalty normalization scale must be positive and finite, got {scale}",
1636 st.name
1637 ))
1638 .into());
1639 }
1640 None => {
1641 return Err(SmoothError::invalid_config(format!(
1642 "{label} term '{}' fused MeasureJet is missing its penalty normalization scale",
1643 st.name
1644 ))
1645 .into());
1646 }
1647 }
1648 }
1649 }
1650 SmoothBasisSpec::Matern { spec, .. } => {
1651 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1652 return Err(SmoothError::invalid_config(format!(
1653 "{label} term '{}' is not frozen: Matern centers must be UserProvided",
1654 st.name
1655 ))
1656 .into());
1657 }
1658 }
1659 SmoothBasisSpec::Duchon { spec, .. } => {
1660 if !matches!(spec.center_strategy, CenterStrategy::UserProvided(_)) {
1661 return Err(SmoothError::invalid_config(format!(
1662 "{label} term '{}' is not frozen: Duchon centers must be UserProvided",
1663 st.name
1664 ))
1665 .into());
1666 }
1667 if matches!(
1668 spec.identifiability,
1669 SpatialIdentifiability::OrthogonalToParametric
1670 ) {
1671 return Err(SmoothError::invalid_config(format!(
1672 "{label} term '{}' is not frozen: Duchon identifiability must be FrozenTransform or None",
1673 st.name
1674 ))
1675 .into());
1676 }
1677 }
1678 SmoothBasisSpec::Pca {
1679 centered,
1680 center_mean,
1681 pca_basis_path,
1682 ..
1683 } => {
1684 if *centered && center_mean.is_none() && pca_basis_path.is_none() {
1685 return Err(SmoothError::invalid_config(format!(
1686 "{label} term '{}' is not frozen: centered Pca missing center_mean",
1687 st.name
1688 ))
1689 .into());
1690 }
1691 }
1692 SmoothBasisSpec::BySmooth { smooth, by_kind } => {
1693 if let SmoothBasisSpec::BySmooth { .. } = smooth.as_ref() {
1694 return Err(format!("{label} term '{}' has nested by-smooths", st.name));
1695 }
1696 match by_kind {
1697 ByVarKind::Numeric { .. } => {}
1698 ByVarKind::Factor { frozen_levels, .. } if frozen_levels.is_none() => {
1699 return Err(format!(
1700 "{label} term '{}' is not frozen: by-factor levels missing",
1701 st.name
1702 ));
1703 }
1704 ByVarKind::Factor { .. } => {}
1705 }
1706 let nested = TermCollectionSpec {
1707 linear_terms: vec![],
1708 random_effect_terms: vec![],
1709 smooth_terms: vec![SmoothTermSpec {
1710 name: st.name.clone(),
1711 basis: (**smooth).clone(),
1712 shape: st.shape,
1713 joint_null_rotation: None,
1714 }],
1715 };
1716 nested.validate_frozen(label)?;
1717 }
1718 SmoothBasisSpec::FactorSmooth { spec } => {
1719 if spec.group_frozen_levels.is_none() {
1720 return Err(format!(
1721 "{label} term '{}' is not frozen: factor-smooth levels missing",
1722 st.name
1723 ));
1724 }
1725 if !matches!(
1726 spec.marginal.knotspec,
1727 BSplineKnotSpec::Provided(_)
1728 | BSplineKnotSpec::PeriodicUniform { .. }
1729 | BSplineKnotSpec::NaturalCubicRegression { .. }
1741 ) {
1742 return Err(format!(
1743 "{label} term '{}' is not frozen: factor-smooth marginal knots missing",
1744 st.name
1745 ));
1746 }
1747 }
1748 SmoothBasisSpec::TensorBSpline { spec, .. } => {
1749 for (dim, marginal) in spec.marginalspecs.iter().enumerate() {
1750 if !matches!(
1751 marginal.knotspec,
1752 BSplineKnotSpec::Provided(_)
1753 | BSplineKnotSpec::PeriodicUniform { .. }
1754 | BSplineKnotSpec::NaturalCubicRegression { .. }
1755 ) {
1756 return Err(SmoothError::invalid_config(format!(
1757 "{label} term '{}' dim {} is not frozen: tensor marginal knotspec must be Provided, PeriodicUniform, or NaturalCubicRegression",
1758 st.name, dim
1759 ))
1760 .into());
1761 }
1762 }
1763 if matches!(
1764 spec.identifiability,
1765 TensorBSplineIdentifiability::SumToZero
1766 | TensorBSplineIdentifiability::MarginalSumToZero
1767 ) {
1768 return Err(SmoothError::invalid_config(format!(
1769 "{label} term '{}' is not frozen: tensor identifiability must be FrozenTransform or None",
1770 st.name
1771 ))
1772 .into());
1773 }
1774 }
1775 }
1776 }
1777
1778 for rt in &self.random_effect_terms {
1779 if rt.frozen_levels.is_none() {
1780 return Err(SmoothError::invalid_config(format!(
1781 "{label} random-effect term '{}' is not frozen: missing frozen_levels",
1782 rt.name
1783 ))
1784 .into());
1785 }
1786 }
1787
1788 Ok(())
1789 }
1790
1791 pub fn remap_feature_columns<E, F>(&self, mut remap: F) -> Result<TermCollectionSpec, E>
1810 where
1811 F: FnMut(usize) -> Result<usize, E>,
1812 {
1813 let mut out = self.clone();
1814 for lt in &mut out.linear_terms {
1815 lt.feature_col = remap(lt.feature_col)?;
1816 for fc in lt.feature_cols.iter_mut() {
1826 *fc = remap(*fc)?;
1827 }
1828 for (col, _bits) in lt.categorical_levels.iter_mut() {
1833 *col = remap(*col)?;
1834 }
1835 }
1836 for rt in &mut out.random_effect_terms {
1837 rt.feature_col = remap(rt.feature_col)?;
1838 }
1839 for st in &mut out.smooth_terms {
1840 remap_smooth_basis_feature_columns(&mut st.basis, &mut remap)?;
1841 }
1842 Ok(out)
1843 }
1844}
1845
1846pub fn remap_smooth_basis_feature_columns<E, F>(
1851 basis: &mut SmoothBasisSpec,
1852 remap: &mut F,
1853) -> Result<(), E>
1854where
1855 F: FnMut(usize) -> Result<usize, E>,
1856{
1857 match basis {
1858 SmoothBasisSpec::ByVariable { inner, by_col, .. }
1859 | SmoothBasisSpec::FactorSumToZero { inner, by_col, .. } => {
1860 *by_col = remap(*by_col)?;
1861 remap_smooth_basis_feature_columns(inner, remap)?;
1862 }
1863 SmoothBasisSpec::BSpline1D { feature_col, .. } => {
1864 *feature_col = remap(*feature_col)?;
1865 }
1866 SmoothBasisSpec::BySmooth { smooth, by_kind } => {
1867 let by_feature_col = match by_kind {
1868 ByVarKind::Numeric { feature_col } | ByVarKind::Factor { feature_col, .. } => {
1869 feature_col
1870 }
1871 };
1872 *by_feature_col = remap(*by_feature_col)?;
1873 remap_smooth_basis_feature_columns(smooth, remap)?;
1874 }
1875 SmoothBasisSpec::FactorSmooth { spec } => {
1876 for fc in spec.continuous_cols.iter_mut() {
1877 *fc = remap(*fc)?;
1878 }
1879 spec.group_col = remap(spec.group_col)?;
1880 }
1881 SmoothBasisSpec::ThinPlate { feature_cols, .. }
1882 | SmoothBasisSpec::Sphere { feature_cols, .. }
1883 | SmoothBasisSpec::ConstantCurvature { feature_cols, .. }
1884 | SmoothBasisSpec::Matern { feature_cols, .. }
1885 | SmoothBasisSpec::MeasureJet { feature_cols, .. }
1886 | SmoothBasisSpec::Duchon { feature_cols, .. }
1887 | SmoothBasisSpec::Pca { feature_cols, .. }
1888 | SmoothBasisSpec::TensorBSpline { feature_cols, .. } => {
1889 for fc in feature_cols.iter_mut() {
1890 *fc = remap(*fc)?;
1891 }
1892 }
1893 }
1894 Ok(())
1895}
1896
1897#[derive(Debug, Clone)]
1898pub enum PenaltyStructureHint {
1899 Ridge(f64),
1900 Kronecker(Vec<Array2<f64>>),
1901}
1902
1903#[derive(Clone)]
1910pub struct BlockwisePenalty {
1911 pub col_range: Range<usize>,
1913 pub local: Array2<f64>,
1916 pub prior_mean: gam_problem::CoefficientPriorMean,
1918 pub structure_hint: Option<PenaltyStructureHint>,
1921 pub op: Option<std::sync::Arc<dyn crate::analytic_penalties::PenaltyOp>>,
1926}
1927
1928impl std::fmt::Debug for BlockwisePenalty {
1929 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
1930 f.debug_struct("BlockwisePenalty")
1931 .field("col_range", &self.col_range)
1932 .field(
1933 "local",
1934 &format_args!("{}×{}", self.local.nrows(), self.local.ncols()),
1935 )
1936 .field("prior_mean", &self.prior_mean)
1937 .field("structure_hint", &self.structure_hint)
1938 .field("op", &self.op.as_ref().map(|o| o.dim()))
1939 .finish()
1940 }
1941}
1942
1943impl BlockwisePenalty {
1944 pub fn new(col_range: Range<usize>, local: Array2<f64>) -> Self {
1946 assert_eq!(col_range.len(), local.nrows());
1947 assert_eq!(col_range.len(), local.ncols());
1948 Self {
1949 col_range,
1950 local,
1951 prior_mean: gam_problem::CoefficientPriorMean::Zero,
1952 structure_hint: None,
1953 op: None,
1954 }
1955 }
1956
1957 pub fn with_prior_mean(
1958 mut self,
1959 prior_mean: gam_problem::CoefficientPriorMean,
1960 ) -> Self {
1961 self.prior_mean = prior_mean;
1962 self
1963 }
1964
1965 pub fn with_op(
1967 mut self,
1968 op: Option<std::sync::Arc<dyn crate::analytic_penalties::PenaltyOp>>,
1969 ) -> Self {
1970 self.op = op;
1971 self
1972 }
1973
1974 pub fn ridge(col_range: Range<usize>, scale: f64) -> Self {
1975 let block_size = col_range.len();
1976 let mut local = Array2::<f64>::zeros((block_size, block_size));
1977 for i in 0..block_size {
1978 local[[i, i]] = scale;
1979 }
1980 Self {
1981 col_range,
1982 local,
1983 prior_mean: gam_problem::CoefficientPriorMean::Zero,
1984 structure_hint: Some(PenaltyStructureHint::Ridge(scale)),
1985 op: None,
1986 }
1987 }
1988
1989 pub fn kronecker(
1990 col_range: Range<usize>,
1991 local: Array2<f64>,
1992 factors: Vec<Array2<f64>>,
1993 ) -> Self {
1994 assert_eq!(col_range.len(), local.nrows());
1995 assert_eq!(col_range.len(), local.ncols());
1996 Self {
1997 col_range,
1998 local,
1999 prior_mean: gam_problem::CoefficientPriorMean::Zero,
2000 structure_hint: Some(PenaltyStructureHint::Kronecker(factors)),
2001 op: None,
2002 }
2003 }
2004
2005 pub fn to_global(&self, p_total: usize) -> Array2<f64> {
2009 let mut g = Array2::<f64>::zeros((p_total, p_total));
2010 let r = &self.col_range;
2011 assert!(
2012 r.end <= p_total && self.local.nrows() == r.len() && self.local.ncols() == r.len(),
2013 "BlockwisePenalty::to_global shape invariant violated: \
2014 col_range={}..{}, local={}x{}, p_total={}",
2015 r.start,
2016 r.end,
2017 self.local.nrows(),
2018 self.local.ncols(),
2019 p_total,
2020 );
2021 g.slice_mut(s![r.start..r.end, r.start..r.end])
2022 .assign(&self.local);
2023 g
2024 }
2025
2026 pub fn to_penalty_matrix(
2029 &self,
2030 total_dim: usize,
2031 ) -> gam_problem::PenaltyMatrix {
2032 gam_problem::PenaltyMatrix::Blockwise {
2033 local: self.local.clone(),
2034 col_range: self.col_range.clone(),
2035 total_dim,
2036 }
2037 }
2038
2039 #[inline]
2041 pub fn block_size(&self) -> usize {
2042 self.col_range.len()
2043 }
2044}
2045
2046pub fn weighted_blockwise_penalty_sum(
2050 penalties: &[BlockwisePenalty],
2051 lambdas: &[f64],
2052 p_total: usize,
2053) -> Array2<f64> {
2054 assert_eq!(penalties.len(), lambdas.len());
2055 for (idx, &lam) in lambdas.iter().enumerate() {
2062 assert!(
2063 lam.is_finite() && lam >= 0.0,
2064 "weighted_blockwise_penalty_sum: lambdas[{idx}] = {lam} is invalid (must be finite and non-negative; negative smoothing parameters violate S_λ ⪰ 0)",
2065 );
2066 }
2067 for (idx, bp) in penalties.iter().enumerate() {
2071 let r = &bp.col_range;
2072 assert!(
2073 r.end <= p_total,
2074 "weighted_blockwise_penalty_sum: penalties[{idx}] col_range {:?} exceeds p_total = {p_total}",
2075 r,
2076 );
2077 }
2078 let mut out = Array2::<f64>::zeros((p_total, p_total));
2079 for (bp, &lam) in penalties.iter().zip(lambdas.iter()) {
2080 let r = &bp.col_range;
2081 let mut slice = out.slice_mut(s![r.start..r.end, r.start..r.end]);
2082 slice.scaled_add(lam, &bp.local);
2083 }
2084 out
2085}
2086
2087#[derive(Debug, Clone)]
2094pub struct KroneckerPenaltySystem {
2095 pub marginal_penalties: Vec<Array2<f64>>,
2097 pub marginal_eigensystems: Vec<(Array1<f64>, Array2<f64>)>,
2099 pub marginal_dims: Vec<usize>,
2101 pub has_double_penalty: bool,
2103}
2104
2105impl KroneckerPenaltySystem {
2106 pub fn new(
2107 marginal_penalties: Vec<Array2<f64>>,
2108 marginal_dims: Vec<usize>,
2109 has_double_penalty: bool,
2110 ) -> Result<Self, BasisError> {
2111 if marginal_penalties.len() != marginal_dims.len() {
2112 crate::bail_dim_basis!(
2113 "KroneckerPenaltySystem: {} penalties vs {} dims",
2114 marginal_penalties.len(),
2115 marginal_dims.len()
2116 );
2117 }
2118 let eigensystems =
2119 kronecker_marginal_eigensystems(&marginal_penalties, "KroneckerPenaltySystem")
2120 .map_err(|e| BasisError::InvalidInput(e.to_string()))?;
2121 Ok(Self {
2122 marginal_penalties,
2123 marginal_eigensystems: eigensystems,
2124 marginal_dims,
2125 has_double_penalty,
2126 })
2127 }
2128
2129 pub fn p_total(&self) -> usize {
2130 self.marginal_dims.iter().copied().product()
2131 }
2132
2133 pub fn ndim(&self) -> usize {
2134 self.marginal_dims.len()
2135 }
2136
2137 pub fn num_penalties(&self) -> usize {
2138 self.marginal_dims.len() + if self.has_double_penalty { 1 } else { 0 }
2139 }
2140
2141 pub fn logdet_and_derivatives(
2145 &self,
2146 lambdas: &[f64],
2147 ridge: f64,
2148 ) -> (f64, Array1<f64>, Array2<f64>) {
2149 let n_pen = self.num_penalties();
2150 assert_eq!(lambdas.len(), n_pen, "lambda count mismatch");
2151 let marginal_evals: Vec<_> = self
2152 .marginal_eigensystems
2153 .iter()
2154 .map(|(evals, _)| evals.view())
2155 .collect();
2156 kronecker_logdet_and_derivatives(
2157 &marginal_evals,
2158 &self.marginal_dims,
2159 lambdas,
2160 self.has_double_penalty,
2161 ridge,
2162 )
2163 }
2164
2165 pub fn logdet_rank_and_derivatives(
2166 &self,
2167 lambdas: &[f64],
2168 ridge: f64,
2169 ) -> (f64, usize, Array1<f64>, Array2<f64>) {
2170 let n_pen = self.num_penalties();
2171 assert_eq!(lambdas.len(), n_pen, "lambda count mismatch");
2172 let d = self.marginal_dims.len();
2173 let mut logdet = 0.0;
2174 let mut rank = 0usize;
2175 let mut grad = Array1::<f64>::zeros(n_pen);
2176 let mut hess = Array2::<f64>::zeros((n_pen, n_pen));
2177 const EIGENVALUE_POSITIVITY_FLOOR: f64 = 1e-12;
2181 const STRUCTURAL_ZERO_FLOOR: f64 = 1e-12;
2185 let mut multi_idx = vec![0usize; d];
2186 loop {
2187 let mut sigma = 0.0;
2188 let mut structural_sigma = 0.0;
2189 for k in 0..d {
2190 let marginal_eigenvalue = self.marginal_eigensystems[k].0[multi_idx[k]];
2191 structural_sigma += marginal_eigenvalue;
2192 sigma += lambdas[k] * marginal_eigenvalue;
2193 }
2194 let joint_null = structural_sigma <= STRUCTURAL_ZERO_FLOOR;
2195 if self.has_double_penalty && joint_null {
2196 sigma += lambdas[d];
2197 }
2198 if structural_sigma > STRUCTURAL_ZERO_FLOOR {
2199 sigma += ridge;
2200 }
2201
2202 if sigma > EIGENVALUE_POSITIVITY_FLOOR {
2203 rank += 1;
2204 logdet += sigma.ln();
2205 let inv_sigma = 1.0 / sigma;
2206 let inv_sigma2 = inv_sigma * inv_sigma;
2207 for k in 0..n_pen {
2208 let ck = if k < d {
2209 lambdas[k] * self.marginal_eigensystems[k].0[multi_idx[k]]
2210 } else if joint_null {
2211 lambdas[d]
2212 } else {
2213 0.0
2214 };
2215 grad[k] += ck * inv_sigma;
2216 hess[[k, k]] += ck * inv_sigma - ck * ck * inv_sigma2;
2217 for l in (k + 1)..n_pen {
2218 let cl = if l < d {
2219 lambdas[l] * self.marginal_eigensystems[l].0[multi_idx[l]]
2220 } else if joint_null {
2221 lambdas[d]
2222 } else {
2223 0.0
2224 };
2225 let off = -ck * cl * inv_sigma2;
2226 hess[[k, l]] += off;
2227 hess[[l, k]] += off;
2228 }
2229 }
2230 }
2231
2232 let mut carry = true;
2233 for dim in (0..d).rev() {
2234 if carry {
2235 multi_idx[dim] += 1;
2236 if multi_idx[dim] < self.marginal_dims[dim] {
2237 carry = false;
2238 } else {
2239 multi_idx[dim] = 0;
2240 }
2241 }
2242 }
2243 if carry {
2244 break;
2245 }
2246 }
2247 (logdet, rank, grad, hess)
2248 }
2249}
2250
2251#[cfg(test)]
2252mod joint_unpenalized_dim_tests {
2253 use super::joint_unpenalized_dim;
2254 use ndarray::{Array2, array};
2255
2256 #[test]
2257 fn no_penalty_is_fully_unpenalized() {
2258 assert_eq!(joint_unpenalized_dim(4, &[], &[]), 4);
2259 }
2260
2261 #[test]
2262 fn single_penalty_returns_its_own_null_space() {
2263 let s = array![[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 5.0]];
2266 assert_eq!(joint_unpenalized_dim(3, std::slice::from_ref(&s), &[2]), 2);
2267 }
2268
2269 #[test]
2270 fn complementary_double_penalty_has_empty_joint_null_space() {
2271 let bending = array![[0.0, 0.0, 0.0], [0.0, 4.0, 0.0], [0.0, 0.0, 4.0]];
2278 let ridge = array![[2.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]];
2279 assert_eq!(joint_unpenalized_dim(3, &[bending, ridge], &[1, 2]), 0);
2280 }
2281
2282 #[test]
2283 fn partial_overlap_keeps_shared_null_direction() {
2284 let a = array![[0.0, 0.0, 0.0], [0.0, 3.0, 0.0], [0.0, 0.0, 0.0]];
2288 let b = array![[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 3.0]];
2289 assert_eq!(joint_unpenalized_dim(3, &[a, b], &[2, 2]), 1);
2290 }
2291
2292 #[test]
2293 fn non_materialized_penalty_falls_back_conservatively() {
2294 let full: Array2<f64> = array![[0.0, 0.0], [0.0, 1.0]];
2298 let factor: Array2<f64> = array![[1.0]]; assert_eq!(
2300 joint_unpenalized_dim(2, &[full, factor.clone()], &[1, 0]),
2301 0
2302 );
2303 assert_eq!(joint_unpenalized_dim(4, std::slice::from_ref(&factor), &[2]), 2);
2305 }
2306}
2307
2308#[cfg(test)]
2309mod kronecker_penalty_system_tests {
2310 use super::KroneckerPenaltySystem;
2311 use ndarray::array;
2312
2313 #[test]
2314 fn double_penalty_rank_derivatives_use_only_joint_null_space() {
2315 let penalties = vec![
2316 array![[0.0, 0.0], [0.0, 2.0]],
2317 array![[0.0, 0.0], [0.0, 3.0]],
2318 ];
2319 let system = KroneckerPenaltySystem::new(penalties, vec![2usize, 2usize], true).unwrap();
2320 let lambdas = vec![5.0, 7.0, 11.0];
2321
2322 let (logdet, rank, grad, hess) = system.logdet_rank_and_derivatives(&lambdas, 0.0);
2323
2324 let expected_diag = [11.0_f64, 21.0, 10.0, 31.0];
2325 let expected_logdet: f64 = expected_diag.iter().map(|v| v.ln()).sum();
2326 assert_eq!(rank, 4);
2327 assert!((logdet - expected_logdet).abs() <= 1e-12);
2328 assert!(
2329 (grad[2] - 1.0).abs() <= 1e-12,
2330 "double-penalty rank derivative must count only the joint null mode, got {}",
2331 grad[2]
2332 );
2333 assert!(hess[[2, 2]].abs() <= 1e-12);
2334 }
2335}
2336
2337#[derive(Clone, Debug)]
2338pub struct TermCollectionDesign {
2339 pub design: DesignMatrix,
2348 pub penalties: Vec<BlockwisePenalty>,
2349 pub nullspace_dims: Vec<usize>,
2350 pub penaltyinfo: Vec<PenaltyBlockInfo>,
2351 pub dropped_penaltyinfo: Vec<DroppedPenaltyBlockInfo>,
2352 pub coefficient_lower_bounds: Option<Array1<f64>>,
2355 pub linear_constraints: Option<LinearInequalityConstraints>,
2358 pub intercept_range: Range<usize>,
2359 pub linear_ranges: Vec<(String, Range<usize>)>,
2360 pub random_effect_ranges: Vec<(String, Range<usize>)>,
2361 pub random_effect_levels: Vec<(String, Vec<u64>)>,
2362 pub smooth: SmoothDesign,
2363}
2364
2365impl TermCollectionDesign {
2366 pub fn leading_penalty_blocks_before_smooth(&self) -> usize {
2374 self.penaltyinfo
2375 .iter()
2376 .take_while(|info| {
2377 matches!(
2378 &info.penalty.source,
2379 crate::basis::PenaltySource::Other(source)
2380 if source == "LinearTermRidge"
2381 || source.starts_with("RandomEffectRidge(")
2382 )
2383 })
2384 .count()
2385 }
2386
2387 pub fn penalties_as_penalty_matrix(&self) -> Vec<gam_problem::PenaltyMatrix> {
2391 let p = self.design.ncols();
2392 self.penalties
2393 .iter()
2394 .map(|bp| bp.to_penalty_matrix(p))
2395 .collect()
2396 }
2397
2398 #[inline]
2400 pub fn num_penalties(&self) -> usize {
2401 self.penalties.len()
2402 }
2403
2404 pub fn realize_coefficient_groups(
2407 &self,
2408 groups: &[CoefficientGroupSpec],
2409 base_prior: &gam_spec::RhoPrior,
2410 ) -> Result<RealizedCoefficientGroups, BasisError> {
2411 realize_coefficient_groups(self, groups, base_prior)
2412 }
2413
2414 pub fn kronecker_penalty_system(&self) -> Option<KroneckerPenaltySystem> {
2425 let [only_term] = self.smooth.terms.as_slice() else {
2426 return None;
2427 };
2428 let kron = only_term.kronecker_factored.as_ref()?;
2429 if kron.marginal_dims.len() < 2
2435 || kron.marginal_penalties.len() != kron.marginal_dims.len()
2436 || kron.marginal_designs.len() != kron.marginal_dims.len()
2437 {
2438 return None;
2439 }
2440 KroneckerPenaltySystem::new(
2441 kron.marginal_penalties.clone(),
2442 kron.marginal_dims.clone(),
2443 kron.has_double_penalty,
2444 )
2445 .ok()
2446 }
2447}
2448
2449#[derive(Clone)]
2455pub struct StandardLatentCoordConfig {
2456 pub values: std::sync::Arc<crate::latent::LatentCoordValues>,
2457 pub term_index: gam_problem::types::SmoothTermIdx,
2458 pub feature_cols: Vec<usize>,
2459 pub manifold: crate::latent::LatentManifold,
2460 pub manifold_auto: bool,
2461 pub retraction_registry: gam_problem::LatentRetractionRegistry,
2462 pub analytic_penalties: Option<std::sync::Arc<crate::AnalyticPenaltyRegistry>>,
2463}
2464
2465#[derive(Clone, Debug, Serialize, Deserialize)]
2466pub struct AdaptiveSpatialMap {
2467 pub termname: String,
2468 pub feature_cols: Vec<usize>,
2469 pub collocation_points: Array2<f64>,
2470 pub inv_magweight: Array1<f64>,
2471 pub invgradweight: Array1<f64>,
2472 pub inv_lapweight: Array1<f64>,
2473}
2474
2475#[derive(Clone, Debug, Serialize, Deserialize)]
2476pub struct AdaptiveRegularizationDiagnostics {
2477 pub epsilon_0: f64,
2478 pub epsilon_g: f64,
2479 pub epsilon_c: f64,
2480 pub epsilon_outer_iterations: usize,
2481 pub mm_iterations: usize,
2482 pub converged: bool,
2483 pub maps: Vec<AdaptiveSpatialMap>,
2484}
2485
2486#[derive(Debug, Clone)]
2487pub struct LinearColumnConditioning {
2488 col_idx: usize,
2489 mean: f64,
2490 scale: f64,
2491}
2492
2493#[derive(Debug, Clone, Default)]
2494pub struct LinearFitConditioning {
2495 pub intercept_idx: usize,
2496 pub columns: Vec<LinearColumnConditioning>,
2497}
2498
2499#[derive(Clone)]
2500pub struct SpatialPsiDerivative {
2501 pub penalty_index: usize,
2503 pub penalty_indices: Vec<usize>,
2504 pub global_range: Range<usize>,
2505 pub total_p: usize,
2506 pub x_psi_local: Array2<f64>,
2507 pub s_psi_components_local: Vec<Array2<f64>>,
2508 pub x_psi_psi_local: Array2<f64>,
2509 pub s_psi_psi_components_local: Vec<Array2<f64>>,
2510 pub aniso_group_id: Option<usize>,
2511 pub aniso_cross_designs: Option<Vec<(usize, Array2<f64>)>>,
2514 pub aniso_cross_penalty_provider: Option<
2518 std::sync::Arc<
2519 dyn Fn(usize) -> Result<Vec<Array2<f64>>, EstimationError> + Send + Sync + 'static,
2520 >,
2521 >,
2522 pub implicit_operator: Option<std::sync::Arc<crate::basis::ImplicitDesignPsiDerivative>>,
2527 pub implicit_axis: usize,
2529}
2530
2531#[derive(Debug, Clone)]
2532pub struct SpatialLogKappaCoords {
2533 pub values: Array1<f64>,
2536 pub dims_per_term: Vec<usize>,
2538}
2539
2540#[derive(Clone, Copy)]
2545pub enum AnisoBoundEnd {
2546 Lower,
2547 Upper,
2548}
2549
2550impl SpatialLogKappaCoords {
2551 pub fn new_with_dims(values: Array1<f64>, dims_per_term: Vec<usize>) -> Self {
2553 assert_eq!(
2554 values.len(),
2555 dims_per_term.iter().sum::<usize>(),
2556 "SpatialLogKappaCoords: values length {} != sum of dims_per_term {}",
2557 values.len(),
2558 dims_per_term.iter().sum::<usize>(),
2559 );
2560 Self {
2561 values,
2562 dims_per_term,
2563 }
2564 }
2565
2566 pub fn from_length_scales(
2568 spec: &TermCollectionSpec,
2569 term_indices: &[usize],
2570 options: &SpatialLengthScaleOptimizationOptions,
2571 ) -> Self {
2572 let mut out = Array1::<f64>::zeros(term_indices.len());
2573 for (slot, &term_idx) in term_indices.iter().enumerate() {
2574 if let Some(cc) = constant_curvature_term_spec(spec, term_idx) {
2580 out[slot] = cc.kappa;
2581 continue;
2582 }
2583 let length_scale = get_spatial_length_scale(spec, term_idx)
2584 .unwrap_or(options.min_length_scale)
2585 .clamp(options.min_length_scale, options.max_length_scale);
2586 out[slot] = -length_scale.ln();
2587 }
2588 Self {
2589 values: out,
2590 dims_per_term: vec![1; term_indices.len()],
2591 }
2592 }
2593
2594 pub fn from_length_scales_aniso(
2612 spec: &TermCollectionSpec,
2613 term_indices: &[usize],
2614 options: &SpatialLengthScaleOptimizationOptions,
2615 ) -> Self {
2616 let mut vals = Vec::new();
2617 let mut dims = Vec::new();
2618 for &term_idx in term_indices {
2619 if let Some(mj) = measure_jet_term_spec(spec, term_idx) {
2623 let seed = measure_jet_psi_seed(mj);
2624 dims.push(seed.len());
2625 vals.extend(seed);
2626 continue;
2627 }
2628 if let Some(cc) = constant_curvature_term_spec(spec, term_idx) {
2634 vals.push(cc.kappa);
2635 dims.push(1);
2636 continue;
2637 }
2638 let length_scale = get_spatial_length_scale(spec, term_idx)
2639 .unwrap_or(options.min_length_scale)
2640 .clamp(options.min_length_scale, options.max_length_scale);
2641 let psi_bar = -length_scale.ln(); if spatial_term_uses_per_axis_psi(spec, term_idx) {
2644 let d = get_spatial_feature_dim(spec, term_idx).unwrap_or(1);
2649 let eta_raw = get_spatial_aniso_log_scales(spec, term_idx)
2650 .expect("predicate guarantees aniso_log_scales is Some");
2651 let eta = center_aniso_log_scales(&eta_raw);
2652 for &eta_a in &eta {
2653 vals.push(psi_bar + eta_a);
2654 }
2655 dims.push(d);
2656 } else {
2657 vals.push(psi_bar);
2664 dims.push(1);
2665 }
2666 }
2667 Self {
2668 values: Array1::from_vec(vals),
2669 dims_per_term: dims,
2670 }
2671 }
2672
2673 pub fn lower_bounds_from_data(
2677 data: ArrayView2<'_, f64>,
2678 spec: &TermCollectionSpec,
2679 term_indices: &[usize],
2680 options: &SpatialLengthScaleOptimizationOptions,
2681 ) -> Self {
2682 let mut values = Array1::<f64>::zeros(term_indices.len());
2683 for (slot, &term_idx) in term_indices.iter().enumerate() {
2684 values[slot] = spatial_term_psi_bounds(data, spec, term_idx, options).0;
2685 }
2686 Self {
2687 values,
2688 dims_per_term: vec![1; term_indices.len()],
2689 }
2690 }
2691
2692 pub fn upper_bounds_from_data(
2694 data: ArrayView2<'_, f64>,
2695 spec: &TermCollectionSpec,
2696 term_indices: &[usize],
2697 options: &SpatialLengthScaleOptimizationOptions,
2698 ) -> Self {
2699 let mut values = Array1::<f64>::zeros(term_indices.len());
2700 for (slot, &term_idx) in term_indices.iter().enumerate() {
2701 values[slot] = spatial_term_psi_bounds(data, spec, term_idx, options).1;
2702 }
2703 Self {
2704 values,
2705 dims_per_term: vec![1; term_indices.len()],
2706 }
2707 }
2708
2709 pub fn lower_bounds_aniso_from_data(
2726 data: ArrayView2<'_, f64>,
2727 spec: &TermCollectionSpec,
2728 term_indices: &[usize],
2729 dims_per_term: &[usize],
2730 options: &SpatialLengthScaleOptimizationOptions,
2731 ) -> Self {
2732 Self::aniso_bounds_from_data(
2733 data,
2734 spec,
2735 term_indices,
2736 dims_per_term,
2737 options,
2738 AnisoBoundEnd::Lower,
2739 )
2740 }
2741
2742 pub fn upper_bounds_aniso_from_data(
2746 data: ArrayView2<'_, f64>,
2747 spec: &TermCollectionSpec,
2748 term_indices: &[usize],
2749 dims_per_term: &[usize],
2750 options: &SpatialLengthScaleOptimizationOptions,
2751 ) -> Self {
2752 Self::aniso_bounds_from_data(
2753 data,
2754 spec,
2755 term_indices,
2756 dims_per_term,
2757 options,
2758 AnisoBoundEnd::Upper,
2759 )
2760 }
2761
2762 fn aniso_bounds_from_data(
2768 data: ArrayView2<'_, f64>,
2769 spec: &TermCollectionSpec,
2770 term_indices: &[usize],
2771 dims_per_term: &[usize],
2772 options: &SpatialLengthScaleOptimizationOptions,
2773 end: AnisoBoundEnd,
2774 ) -> Self {
2775 assert_eq!(term_indices.len(), dims_per_term.len());
2776 let total: usize = dims_per_term.iter().sum();
2777 let mut values = Array1::<f64>::zeros(total);
2778 let mut cursor = 0;
2779 for (slot, &term_idx) in term_indices.iter().enumerate() {
2780 let d = dims_per_term[slot];
2781 if let Some(mj) = measure_jet_term_spec(spec, term_idx) {
2784 let bounds = measure_jet_psi_bound_values(mj, matches!(end, AnisoBoundEnd::Upper));
2785 for (offset, bound) in bounds.into_iter().enumerate() {
2786 if offset < d {
2787 values[cursor + offset] = bound;
2788 }
2789 }
2790 cursor += d;
2791 continue;
2792 }
2793 if constant_curvature_term_spec(spec, term_idx).is_some() {
2796 let (lo, hi) = constant_curvature_kappa_bounds(data, spec, term_idx);
2797 if d >= 1 {
2798 values[cursor] = match end {
2799 AnisoBoundEnd::Lower => lo,
2800 AnisoBoundEnd::Upper => hi,
2801 };
2802 }
2803 cursor += d;
2804 continue;
2805 }
2806 let psi_bound = {
2807 let (lo, hi) = spatial_term_psi_bounds(data, spec, term_idx, options);
2808 match end {
2809 AnisoBoundEnd::Lower => lo,
2810 AnisoBoundEnd::Upper => hi,
2811 }
2812 };
2813 let axis_offsets = if d <= 1 {
2814 vec![0.0; d]
2815 } else {
2816 get_spatial_aniso_log_scales(spec, term_idx)
2817 .filter(|eta| eta.len() == d)
2818 .map(|eta| center_aniso_log_scales(&eta))
2819 .unwrap_or_else(|| vec![0.0; d])
2820 };
2821 for offset in 0..d {
2822 values[cursor + offset] = psi_bound + axis_offsets[offset];
2823 }
2824 cursor += d;
2825 }
2826 Self {
2827 values,
2828 dims_per_term: dims_per_term.to_vec(),
2829 }
2830 }
2831
2832 pub fn reseed_from_data(
2841 mut self,
2842 data: ArrayView2<'_, f64>,
2843 spec: &TermCollectionSpec,
2844 term_indices: &[usize],
2845 options: &SpatialLengthScaleOptimizationOptions,
2846 ) -> Self {
2847 assert_eq!(term_indices.len(), self.dims_per_term.len());
2848 let mut cursor = 0;
2849 for (slot, &term_idx) in term_indices.iter().enumerate() {
2850 let d = self.dims_per_term[slot];
2851 if measure_jet_term_spec(spec, term_idx).is_some() {
2854 cursor += d;
2855 continue;
2856 }
2857 if constant_curvature_term_spec(spec, term_idx).is_some() {
2861 cursor += d;
2862 continue;
2863 }
2864 let Some(psi_bar_new) = spatial_term_psi_seed(data, spec, term_idx, options) else {
2865 cursor += d;
2866 continue;
2867 };
2868 if d == 0 {
2869 continue;
2870 }
2871 let current: Vec<f64> = self.values.slice(s![cursor..cursor + d]).to_vec();
2872 let psi_bar_old = current.iter().sum::<f64>() / d as f64;
2873 for (offset, &old_value) in current.iter().enumerate() {
2874 self.values[cursor + offset] = psi_bar_new + (old_value - psi_bar_old);
2875 }
2876 cursor += d;
2877 }
2878 self
2879 }
2880
2881 pub fn clamp_to_bounds(
2892 mut self,
2893 lower: &SpatialLogKappaCoords,
2894 upper: &SpatialLogKappaCoords,
2895 ) -> Self {
2896 assert_eq!(self.values.len(), lower.values.len());
2897 assert_eq!(self.values.len(), upper.values.len());
2898 let mut n_projected = 0usize;
2899 let mut worst_delta = 0.0_f64;
2900 for idx in 0..self.values.len() {
2901 let lo = lower.values[idx];
2902 let hi = upper.values[idx];
2903 if !(lo.is_finite() && hi.is_finite()) {
2904 continue;
2905 }
2906 let v = self.values[idx];
2907 if v < lo {
2908 worst_delta = worst_delta.max(lo - v);
2909 self.values[idx] = lo;
2910 n_projected += 1;
2911 } else if v > hi {
2912 worst_delta = worst_delta.max(v - hi);
2913 self.values[idx] = hi;
2914 n_projected += 1;
2915 }
2916 }
2917 if n_projected > 0 {
2918 log::info!(
2919 "[spatial-kappa] projected {n_projected}/{} ψ seed coords into data-derived bounds \
2920 (worst excess={worst_delta:.3} log units); user length_scale falls outside \
2921 [{KERNEL_RANGE_MIN_DIAMETER_FRACTION}/r_max, {KERNEL_RANGE_MAX_SPACING_MULTIPLE}/r_min] geometry window",
2922 self.values.len()
2923 );
2924 }
2925 self
2926 }
2927
2928 pub fn from_theta_tail_with_dims(
2930 theta: &Array1<f64>,
2931 start: usize,
2932 dims_per_term: Vec<usize>,
2933 ) -> Self {
2934 let total: usize = dims_per_term.iter().sum();
2935 Self {
2936 values: theta.slice(s![start..start + total]).to_owned(),
2937 dims_per_term,
2938 }
2939 }
2940
2941 pub fn len(&self) -> usize {
2943 self.values.len()
2944 }
2945
2946 pub fn dims_per_term(&self) -> &[usize] {
2948 &self.dims_per_term
2949 }
2950
2951 fn term_offset(&self, term_idx: usize) -> usize {
2953 self.dims_per_term[..term_idx].iter().sum()
2954 }
2955
2956 pub fn term_slice(&self, term_idx: usize) -> &[f64] {
2958 let offset = self.term_offset(term_idx);
2959 let d = self.dims_per_term[term_idx];
2960 &self.values.as_slice().unwrap()[offset..offset + d]
2961 }
2962
2963 pub fn as_array(&self) -> &Array1<f64> {
2964 &self.values
2965 }
2966
2967 pub fn set_scalar_slot(&mut self, slot: usize, value: f64) -> bool {
2973 if slot >= self.dims_per_term.len() || self.dims_per_term[slot] != 1 {
2974 return false;
2975 }
2976 let offset = self.term_offset(slot);
2977 self.values[offset] = value;
2978 true
2979 }
2980
2981 pub fn split_at(&self, mid: usize) -> (Self, Self) {
2984 let flat_mid: usize = self.dims_per_term[..mid].iter().sum();
2985 (
2986 Self {
2987 values: self.values.slice(s![0..flat_mid]).to_owned(),
2988 dims_per_term: self.dims_per_term[..mid].to_vec(),
2989 },
2990 Self {
2991 values: self.values.slice(s![flat_mid..]).to_owned(),
2992 dims_per_term: self.dims_per_term[mid..].to_vec(),
2993 },
2994 )
2995 }
2996
2997 pub fn apply_tospec(
3004 &self,
3005 spec: &TermCollectionSpec,
3006 term_indices: &[usize],
3007 ) -> Result<TermCollectionSpec, EstimationError> {
3008 if term_indices.len() != self.dims_per_term.len() {
3009 crate::bail_invalid_estim!(
3010 "SpatialLogKappaCoords::apply_tospec: term count mismatch: \
3011 term_indices={} dims_per_term={}",
3012 term_indices.len(),
3013 self.dims_per_term.len()
3014 );
3015 }
3016 let mut updated = spec.clone();
3017 for (slot, &term_idx) in term_indices.iter().enumerate() {
3018 let psi = self.term_slice(slot);
3019 let d = self.dims_per_term[slot];
3020 if measure_jet_term_spec(&updated, term_idx).is_some() {
3023 set_measure_jet_psi_dials(&mut updated, term_idx, psi)?;
3024 continue;
3025 }
3026 if constant_curvature_term_spec(&updated, term_idx).is_some() {
3030 set_constant_curvature_kappa(&mut updated, term_idx, psi)?;
3031 continue;
3032 }
3033 let (next_length_scale, next_aniso) = spatial_term_psi_to_length_scale_and_aniso(psi);
3034 if (d == 1 || next_length_scale.is_some())
3035 && let Some(length_scale) = next_length_scale
3036 {
3037 set_spatial_length_scale(&mut updated, term_idx, length_scale)?;
3038 }
3039 if let Some(eta) = next_aniso {
3040 set_spatial_aniso_log_scales(&mut updated, term_idx, eta)?;
3041 }
3042 }
3043 Ok(updated)
3044 }
3045}
3046
3047pub fn center_aniso_log_scales(eta: &[f64]) -> Vec<f64> {
3048 if eta.len() <= 1 {
3049 return eta.to_vec();
3050 }
3051 let mean = eta.iter().sum::<f64>() / eta.len() as f64;
3052 eta.iter()
3053 .map(|&v| {
3054 let centered = v - mean;
3055 if centered.abs() <= 1e-15 {
3056 0.0
3057 } else {
3058 centered
3059 }
3060 })
3061 .collect()
3062}
3063
3064pub fn spatial_term_uses_per_axis_psi(resolvedspec: &TermCollectionSpec, term_idx: usize) -> bool {
3067 if let Some(mj) = measure_jet_term_spec(resolvedspec, term_idx) {
3068 return measure_jet_enrolls_psi(mj);
3069 }
3070 let Some(d) = get_spatial_feature_dim(resolvedspec, term_idx) else {
3071 return false;
3072 };
3073 if d <= 1 {
3074 return false;
3075 }
3076 let Some(eta) = get_spatial_aniso_log_scales(resolvedspec, term_idx) else {
3077 return false;
3078 };
3079 if eta.len() != d {
3080 return false;
3081 }
3082 !matches!(
3083 resolvedspec.smooth_terms.get(term_idx).map(|term| &term.basis),
3084 Some(SmoothBasisSpec::Duchon { .. })
3085 )
3086}
3087
3088pub fn set_spatial_length_scale(
3089 spec: &mut TermCollectionSpec,
3090 term_idx: usize,
3091 length_scale: f64,
3092) -> Result<(), EstimationError> {
3093 let Some(term) = spec.smooth_terms.get_mut(term_idx) else {
3094 crate::bail_invalid_estim!("spatial length-scale term index {term_idx} out of range");
3095 };
3096 match &mut term.basis {
3097 SmoothBasisSpec::ThinPlate { spec, .. } => {
3098 spec.length_scale = length_scale;
3099 Ok(())
3100 }
3101 SmoothBasisSpec::Matern { spec, .. } => {
3102 spec.length_scale = length_scale;
3103 Ok(())
3104 }
3105 SmoothBasisSpec::Duchon { spec, .. } => {
3106 spec.length_scale = Some(length_scale);
3107 Ok(())
3108 }
3109 _ => Err(EstimationError::InvalidInput(format!(
3110 "term '{}' does not expose a spatial length scale",
3111 term.name
3112 ))),
3113 }
3114}
3115
3116pub fn get_spatial_length_scale(spec: &TermCollectionSpec, term_idx: usize) -> Option<f64> {
3117 spec.smooth_terms
3118 .get(term_idx)
3119 .and_then(|term| match &term.basis {
3120 SmoothBasisSpec::ThinPlate { spec, .. } => Some(spec.length_scale),
3121 SmoothBasisSpec::Matern { spec, .. } => Some(spec.length_scale),
3122 SmoothBasisSpec::Duchon { spec, .. } => spec.length_scale,
3123 _ => None,
3124 })
3125}
3126
3127pub fn spatial_term_supports_hyper_optimization(spec: &TermCollectionSpec, term_idx: usize) -> bool {
3128 if let Some(term) = spec.smooth_terms.get(term_idx)
3134 && let SmoothBasisSpec::ThinPlate { .. } = &term.basis
3135 {
3136 return false;
3137 }
3138
3139 if let Some(term) = spec.smooth_terms.get(term_idx)
3164 && let SmoothBasisSpec::Matern { .. } = &term.basis
3165 {
3166 return true;
3167 }
3168
3169 if let Some(mj) = measure_jet_term_spec(spec, term_idx) {
3172 return measure_jet_enrolls_psi(mj);
3173 }
3174
3175 if constant_curvature_term_spec(spec, term_idx).is_some() {
3182 return true;
3183 }
3184
3185 get_spatial_length_scale(spec, term_idx).is_some()
3186}
3187
3188pub fn measure_jet_term_spec(
3191 spec: &TermCollectionSpec,
3192 term_idx: usize,
3193) -> Option<&crate::basis::MeasureJetBasisSpec> {
3194 spec.smooth_terms
3195 .get(term_idx)
3196 .and_then(|term| match &term.basis {
3197 SmoothBasisSpec::MeasureJet { spec, .. } => Some(spec),
3198 _ => None,
3199 })
3200}
3201
3202pub fn measure_jet_enrolls_psi(mj: &crate::basis::MeasureJetBasisSpec) -> bool {
3209 measure_jet_learns_length_scale(mj)
3218 || (mj.tau0 > 0.0 && crate::basis::measure_jet_multiscale_mode(mj))
3219}
3220
3221pub fn measure_jet_learns_length_scale(mj: &crate::basis::MeasureJetBasisSpec) -> bool {
3224 mj.learn_length_scale
3225}
3226
3227pub fn freeze_measure_jet_length_scale_learning(spec: &mut TermCollectionSpec) -> usize {
3228 let mut frozen = 0;
3229 for term in spec.smooth_terms.iter_mut() {
3230 if let SmoothBasisSpec::MeasureJet { spec: mj, .. } = &mut term.basis
3231 && mj.learn_length_scale
3232 {
3233 mj.learn_length_scale = false;
3234 frozen += 1;
3235 }
3236 }
3237 frozen
3238}
3239
3240pub const MEASURE_JET_PSI_ALPHA_BOUNDS: (f64, f64) = (-1.0, 3.0);
3248
3249pub const MEASURE_JET_PSI_LN_TAU_BOUNDS: (f64, f64) = (-18.420680743952367, 4.605170185988092);
3250
3251pub const MEASURE_JET_PSI_LN_LENGTH_SCALE_BOUNDS: (f64, f64) = (-6.907755278982137, 4.605170185988092);
3257
3258pub fn measure_jet_penalty_psi_dim(mj: &crate::basis::MeasureJetBasisSpec) -> usize {
3266 if crate::basis::measure_jet_multiscale_mode(mj) {
3267 2
3268 } else {
3269 0
3270 }
3271}
3272
3273pub fn measure_jet_psi_dim(mj: &crate::basis::MeasureJetBasisSpec) -> usize {
3277 usize::from(measure_jet_learns_length_scale(mj)) + measure_jet_penalty_psi_dim(mj)
3278}
3279
3280pub fn measure_jet_psi_seed(mj: &crate::basis::MeasureJetBasisSpec) -> Vec<f64> {
3285 let mut seed = Vec::with_capacity(measure_jet_psi_dim(mj));
3286 if measure_jet_learns_length_scale(mj) {
3287 let ell = if mj.length_scale > 0.0 {
3291 mj.length_scale
3292 } else {
3293 1.0
3294 };
3295 seed.push(ell.ln());
3296 }
3297 if measure_jet_penalty_psi_dim(mj) > 0 {
3298 let ln_tau = mj.tau0.max(f64::MIN_POSITIVE).ln();
3300 seed.extend_from_slice(&[mj.alpha, ln_tau]);
3301 }
3302 seed
3303}
3304
3305pub fn measure_jet_psi_bound_values(mj: &crate::basis::MeasureJetBasisSpec, upper: bool) -> Vec<f64> {
3308 let pick = |b: (f64, f64)| if upper { b.1 } else { b.0 };
3309 let mut bounds = Vec::with_capacity(measure_jet_psi_dim(mj));
3310 if measure_jet_learns_length_scale(mj) {
3311 bounds.push(pick(MEASURE_JET_PSI_LN_LENGTH_SCALE_BOUNDS));
3312 }
3313 if measure_jet_penalty_psi_dim(mj) > 0 {
3314 bounds.push(pick(MEASURE_JET_PSI_ALPHA_BOUNDS));
3316 bounds.push(pick(MEASURE_JET_PSI_LN_TAU_BOUNDS));
3317 }
3318 bounds
3319}
3320
3321pub fn apply_measure_jet_psi(
3326 mj: &mut crate::basis::MeasureJetBasisSpec,
3327 psi: &[f64],
3328) -> Result<bool, EstimationError> {
3329 if psi.len() != measure_jet_psi_dim(mj) {
3330 crate::bail_invalid_estim!(
3331 "measure-jet ψ write-back dimension mismatch: got {} values for a {}-dial term",
3332 psi.len(),
3333 measure_jet_psi_dim(mj)
3334 );
3335 }
3336 let mut changed = false;
3337 let mut cursor = 0usize;
3341 if measure_jet_learns_length_scale(mj) {
3342 let next_ell = psi[cursor].exp();
3343 cursor += 1;
3344 if !(next_ell.is_finite() && next_ell > 0.0) {
3345 crate::bail_invalid_estim!(
3346 "measure-jet ψ write-back produced a non-finite/non-positive length_scale (ℓ={next_ell})"
3347 );
3348 }
3349 if next_ell != mj.length_scale {
3350 mj.length_scale = next_ell;
3351 changed = true;
3352 }
3353 }
3354 if measure_jet_penalty_psi_dim(mj) > 0 {
3355 let next_alpha = psi[cursor];
3358 let next_tau = psi[cursor + 1].exp();
3359 if !(next_alpha.is_finite() && next_tau.is_finite() && next_tau > 0.0) {
3360 crate::bail_invalid_estim!(
3361 "measure-jet ψ write-back produced non-finite dials (alpha={next_alpha}, tau={next_tau})"
3362 );
3363 }
3364 if next_alpha != mj.alpha {
3365 mj.alpha = next_alpha;
3366 changed = true;
3367 }
3368 if next_tau != mj.tau0 {
3369 mj.tau0 = next_tau;
3370 changed = true;
3371 }
3372 }
3373 Ok(changed)
3374}
3375
3376pub fn set_measure_jet_psi_dials(
3379 spec: &mut TermCollectionSpec,
3380 term_idx: usize,
3381 psi: &[f64],
3382) -> Result<bool, EstimationError> {
3383 let Some(term) = spec.smooth_terms.get_mut(term_idx) else {
3384 crate::bail_invalid_estim!("measure-jet ψ write-back: term index {term_idx} out of range");
3385 };
3386 set_single_term_measure_jet_psi_dials(term, psi)
3387}
3388
3389pub fn set_single_term_measure_jet_psi_dials(
3394 term: &mut SmoothTermSpec,
3395 psi: &[f64],
3396) -> Result<bool, EstimationError> {
3397 let SmoothBasisSpec::MeasureJet { spec: mj, .. } = &mut term.basis else {
3398 crate::bail_invalid_estim!("measure-jet ψ write-back targeted a non-measure-jet term");
3399 };
3400 apply_measure_jet_psi(mj, psi)
3401}
3402
3403pub fn constant_curvature_term_spec(
3406 spec: &TermCollectionSpec,
3407 term_idx: usize,
3408) -> Option<&crate::basis::ConstantCurvatureBasisSpec> {
3409 spec.smooth_terms
3410 .get(term_idx)
3411 .and_then(|term| match &term.basis {
3412 SmoothBasisSpec::ConstantCurvature { spec, .. } => Some(spec),
3413 _ => None,
3414 })
3415}
3416
3417pub const CONSTANT_CURVATURE_KAPPA_CHART_FRACTION: f64 = 0.5;
3425
3426pub const CONSTANT_CURVATURE_MIN_CHART_RADIUS2: f64 = 1e-8;
3430
3431pub fn constant_curvature_kappa_bounds(
3436 data: ArrayView2<'_, f64>,
3437 spec: &TermCollectionSpec,
3438 term_idx: usize,
3439) -> (f64, f64) {
3440 let feature_cols = match spec.smooth_terms.get(term_idx).map(|t| &t.basis) {
3441 Some(SmoothBasisSpec::ConstantCurvature { feature_cols, .. }) => feature_cols,
3442 _ => return (-1.0, 1.0),
3443 };
3444 let mut max_r2 = CONSTANT_CURVATURE_MIN_CHART_RADIUS2;
3445 for row in data.outer_iter() {
3446 let mut r2 = 0.0_f64;
3447 for &c in feature_cols.iter() {
3448 if let Some(&v) = row.get(c)
3449 && v.is_finite()
3450 {
3451 r2 += v * v;
3452 }
3453 }
3454 if r2 > max_r2 {
3455 max_r2 = r2;
3456 }
3457 }
3458 let half = CONSTANT_CURVATURE_KAPPA_CHART_FRACTION / max_r2;
3459 (-half, half)
3460}
3461
3462pub fn set_constant_curvature_kappa(
3466 spec: &mut TermCollectionSpec,
3467 term_idx: usize,
3468 psi: &[f64],
3469) -> Result<bool, EstimationError> {
3470 let Some(term) = spec.smooth_terms.get_mut(term_idx) else {
3471 crate::bail_invalid_estim!(
3472 "constant-curvature κ write-back: term index {term_idx} out of range"
3473 );
3474 };
3475 set_single_term_constant_curvature_kappa(term, psi)
3476}
3477
3478pub fn set_single_term_constant_curvature_kappa(
3483 term: &mut SmoothTermSpec,
3484 psi: &[f64],
3485) -> Result<bool, EstimationError> {
3486 if psi.len() != 1 {
3487 crate::bail_invalid_estim!(
3488 "constant-curvature κ write-back expects exactly one value, got {}",
3489 psi.len()
3490 );
3491 }
3492 let next_kappa = psi[0];
3493 if !next_kappa.is_finite() {
3494 crate::bail_invalid_estim!(
3495 "constant-curvature κ write-back produced a non-finite κ = {next_kappa}"
3496 );
3497 }
3498 let SmoothBasisSpec::ConstantCurvature { spec: cc, .. } = &mut term.basis else {
3499 crate::bail_invalid_estim!(
3500 "constant-curvature κ write-back targeted a non-constant-curvature term"
3501 );
3502 };
3503 if cc.kappa != next_kappa {
3504 cc.kappa = next_kappa;
3505 Ok(true)
3506 } else {
3507 Ok(false)
3508 }
3509}
3510
3511pub fn spatial_term_has_locked_kappa(spec: &TermCollectionSpec, term_idx: usize) -> bool {
3522 get_spatial_length_scale(spec, term_idx).is_some()
3523 && !spatial_term_uses_per_axis_psi(spec, term_idx)
3524}
3525
3526pub fn all_spatial_terms_kappa_fixed(spec: &TermCollectionSpec) -> bool {
3527 spec.smooth_terms.iter().enumerate().all(|(idx, _)| {
3528 !spatial_term_supports_hyper_optimization(spec, idx)
3529 || spatial_term_has_locked_kappa(spec, idx)
3530 })
3531}
3532
3533pub fn spatial_identifiability_policy(termspec: &SmoothTermSpec) -> Option<&SpatialIdentifiability> {
3534 match &termspec.basis {
3535 SmoothBasisSpec::ThinPlate { spec, .. } => Some(&spec.identifiability),
3536 SmoothBasisSpec::Duchon { spec, .. } => Some(&spec.identifiability),
3537 _ => None,
3538 }
3539}
3540
3541pub const NULLSPACE_WELLDET_DEGENERACY_RHO_SD: f64 = 15.0;
3545
3546pub fn is_nullspace_degeneracy_prior(prior: &gam_spec::RhoPrior) -> bool {
3549 matches!(
3550 prior,
3551 gam_spec::RhoPrior::Normal { mean, sd }
3552 if *mean == 0.0 && *sd == NULLSPACE_WELLDET_DEGENERACY_RHO_SD
3553 )
3554}
3555
3556pub const KERNEL_RANGE_MIN_DIAMETER_FRACTION: f64 = 2.0;
3568
3569pub const KERNEL_RANGE_MAX_SPACING_MULTIPLE: f64 = 1e2;
3574
3575
3576pub fn spatial_term_psi_bounds(
3585 data: ArrayView2<'_, f64>,
3586 spec: &TermCollectionSpec,
3587 term_idx: usize,
3588 options: &SpatialLengthScaleOptimizationOptions,
3589) -> (f64, f64) {
3590 let fallback = (
3591 -options.max_length_scale.ln(),
3592 -options.min_length_scale.ln(),
3593 );
3594 if constant_curvature_term_spec(spec, term_idx).is_some() {
3599 return constant_curvature_kappa_bounds(data, spec, term_idx);
3600 }
3601 let Some(term) = spec.smooth_terms.get(term_idx) else {
3602 return fallback;
3603 };
3604 let aniso = get_spatial_aniso_log_scales(spec, term_idx);
3617 let r_bounds = match spatial_term_center_strategy(term) {
3618 Some(CenterStrategy::UserProvided(centers)) if centers.nrows() >= 2 => {
3619 match aniso.as_deref() {
3620 Some(eta) if eta.len() == centers.ncols() => {
3621 let y = points_in_aniso_y_space(centers.view(), eta);
3622 pairwise_distance_bounds(y.view())
3623 }
3624 _ => pairwise_distance_bounds(centers.view()),
3625 }
3626 }
3627 _ => standardized_spatial_term_data(data, term)
3628 .ok()
3629 .and_then(|x| match aniso.as_deref() {
3630 Some(eta) if eta.len() == x.ncols() => {
3631 let y = points_in_aniso_y_space(x.view(), eta);
3632 pairwise_distance_bounds_sampled(y.view())
3633 }
3634 _ => pairwise_distance_bounds_sampled(x.view()),
3635 }),
3636 };
3637 let Some((r_min, r_max)) = r_bounds else {
3638 return fallback;
3639 };
3640 let psi_lo_data = (KERNEL_RANGE_MIN_DIAMETER_FRACTION / r_max).ln();
3646 let psi_hi_data = (KERNEL_RANGE_MAX_SPACING_MULTIPLE / r_min).ln();
3647 let psi_lo = psi_lo_data.max(fallback.0);
3657 let psi_hi = psi_hi_data.min(fallback.1);
3658 if psi_lo >= psi_hi {
3659 return fallback;
3662 }
3663 (psi_lo, psi_hi)
3664}
3665
3666pub fn spatial_term_psi_seed(
3670 data: ArrayView2<'_, f64>,
3671 spec: &TermCollectionSpec,
3672 term_idx: usize,
3673 options: &SpatialLengthScaleOptimizationOptions,
3674) -> Option<f64> {
3675 if get_spatial_length_scale(spec, term_idx).is_some() {
3676 return None; }
3678 let (psi_lo, psi_hi) = spatial_term_psi_bounds(data, spec, term_idx, options);
3679 Some(0.5 * (psi_lo + psi_hi))
3680}
3681
3682pub fn spatial_term_psi_to_length_scale_and_aniso(psi: &[f64]) -> (Option<f64>, Option<Vec<f64>>) {
3683 if psi.len() <= 1 {
3684 (Some((-psi.first().copied().unwrap_or(0.0)).exp()), None)
3685 } else {
3686 let psi_bar = psi.iter().sum::<f64>() / psi.len() as f64;
3687 (
3688 Some((-psi_bar).exp()),
3689 Some(psi.iter().map(|&value| value - psi_bar).collect()),
3690 )
3691 }
3692}
3693
3694pub fn get_spatial_aniso_log_scales(
3696 spec: &TermCollectionSpec,
3697 term_idx: usize,
3698) -> Option<Vec<f64>> {
3699 spec.smooth_terms
3700 .get(term_idx)
3701 .and_then(|term| match &term.basis {
3702 SmoothBasisSpec::Matern { spec, .. } => spec.aniso_log_scales.clone(),
3703 SmoothBasisSpec::Duchon { spec, .. } => spec.aniso_log_scales.clone(),
3704 _ => None,
3705 })
3706}
3707
3708pub fn response_aware_axis_contrasts(
3728 x: ndarray::ArrayView2<'_, f64>,
3729 y: ndarray::ArrayView1<'_, f64>,
3730) -> Option<Vec<f64>> {
3731 let n = x.nrows();
3732 let d = x.ncols();
3733 if d <= 1 || n < 4 || y.len() != n {
3734 return None;
3735 }
3736 if x.iter().any(|v| !v.is_finite()) || y.iter().any(|v| !v.is_finite()) {
3737 return None;
3738 }
3739 let mut scores = Vec::with_capacity(d);
3740 for a in 0..d {
3741 let mut order: Vec<usize> = (0..n).collect();
3742 let col = x.column(a);
3743 order.sort_by(|&i, &j| {
3744 col[i]
3745 .partial_cmp(&col[j])
3746 .unwrap_or(std::cmp::Ordering::Equal)
3747 });
3748 let mut tv = 0.0_f64;
3749 for w in order.windows(2) {
3750 let diff = y[w[1]] - y[w[0]];
3751 tv += diff * diff;
3752 }
3753 scores.push(-0.5 * (tv + 1e-12).ln());
3755 }
3756 if scores.iter().any(|v| !v.is_finite()) {
3757 return None;
3758 }
3759 let mean = scores.iter().sum::<f64>() / d as f64;
3760 let centered: Vec<f64> = scores.iter().map(|&s| s - mean).collect();
3761 if centered.iter().all(|&v| v.abs() < 1e-9) {
3764 return None;
3765 }
3766 Some(centered)
3767}
3768
3769pub fn apply_response_aware_anisotropy_seed(
3778 data: ArrayView2<'_, f64>,
3779 y: ndarray::ArrayView1<'_, f64>,
3780 spec: &mut TermCollectionSpec,
3781 spatial_terms: &[usize],
3782) {
3783 const MAX_NUDGE: f64 = std::f64::consts::LN_2;
3788 for &term_idx in spatial_terms {
3789 let Some(current_eta) = get_spatial_aniso_log_scales(spec, term_idx) else {
3790 continue;
3791 };
3792 let d = current_eta.len();
3793 if d <= 1 {
3794 continue;
3795 }
3796 let Some(term) = spec.smooth_terms.get(term_idx) else {
3797 continue;
3798 };
3799 let feature_cols = term.basis.structural_feature_cols();
3800 if feature_cols.len() != d {
3801 continue;
3802 }
3803 let Ok(x) = select_columns(data, &feature_cols) else {
3804 continue;
3805 };
3806 let Some(contrast) = response_aware_axis_contrasts(x.view(), y) else {
3807 continue;
3808 };
3809 let nudged: Vec<f64> = current_eta
3810 .iter()
3811 .zip(contrast.iter())
3812 .map(|(&eta_a, &c_a)| eta_a + c_a.clamp(-MAX_NUDGE, MAX_NUDGE))
3813 .collect();
3814 if let Err(err) = set_spatial_aniso_log_scales(spec, term_idx, nudged) {
3817 log::debug!(
3818 "[spatial-kappa] response-aware anisotropy seed skipped for term {term_idx}: {err}"
3819 );
3820 }
3821 }
3822}
3823
3824pub fn get_spatial_feature_dim(spec: &TermCollectionSpec, term_idx: usize) -> Option<usize> {
3826 spec.smooth_terms
3827 .get(term_idx)
3828 .and_then(|term| match &term.basis {
3829 SmoothBasisSpec::ThinPlate { feature_cols, .. } => Some(feature_cols.len()),
3830 SmoothBasisSpec::Matern { feature_cols, .. } => Some(feature_cols.len()),
3831 SmoothBasisSpec::Duchon { feature_cols, .. } => Some(feature_cols.len()),
3832 _ => None,
3833 })
3834}
3835
3836pub fn log_spatial_aniso_scales(spec: &TermCollectionSpec) {
3843 for (term_idx, term) in spec.smooth_terms.iter().enumerate() {
3844 let (aniso, length_scale) = match &term.basis {
3845 SmoothBasisSpec::Matern { spec, .. } => {
3846 (spec.aniso_log_scales.as_ref(), Some(spec.length_scale))
3847 }
3848 SmoothBasisSpec::Duchon { spec, .. } => {
3849 (spec.aniso_log_scales.as_ref(), spec.length_scale)
3850 }
3851 _ => (None, None),
3852 };
3853 let Some(eta) = aniso else { continue };
3854 if eta.is_empty() {
3855 continue;
3856 }
3857 let mut lines = match length_scale {
3858 Some(ls) => format!(
3859 "[spatial-kappa] term {} (\"{}\"): anisotropic length scales optimized (global length_scale={:.4})",
3860 term_idx, term.name, ls
3861 ),
3862 None => format!(
3863 "[spatial-kappa] term {} (\"{}\"): pure Duchon shape anisotropy optimized",
3864 term_idx, term.name
3865 ),
3866 };
3867 for (a, &eta_a) in eta.iter().enumerate() {
3868 if let Some(ls) = length_scale {
3869 let length_a = ls * (-eta_a).exp();
3870 let kappa_a = (1.0 / ls) * eta_a.exp();
3871 lines.push_str(&format!(
3872 "\n axis {}: eta={:+.4}, length={:.4}, kappa={:.4}",
3873 a, eta_a, length_a, kappa_a
3874 ));
3875 } else {
3876 lines.push_str(&format!("\n axis {}: eta={:+.4}", a, eta_a));
3877 }
3878 }
3879 log::info!("{}", lines);
3880 }
3881}
3882
3883pub fn set_spatial_aniso_log_scales(
3885 spec: &mut TermCollectionSpec,
3886 term_idx: usize,
3887 eta: Vec<f64>,
3888) -> Result<(), EstimationError> {
3889 let eta = center_aniso_log_scales(&eta);
3890 let Some(term) = spec.smooth_terms.get_mut(term_idx) else {
3891 crate::bail_invalid_estim!("spatial aniso_log_scales term index {term_idx} out of range");
3892 };
3893 match &mut term.basis {
3894 SmoothBasisSpec::Matern { spec, .. } => {
3895 spec.aniso_log_scales = Some(eta);
3896 Ok(())
3897 }
3898 SmoothBasisSpec::Duchon { spec, .. } => {
3899 spec.aniso_log_scales = Some(eta);
3900 Ok(())
3901 }
3902 _ => Err(EstimationError::InvalidInput(format!(
3903 "term '{}' does not support aniso_log_scales",
3904 term.name
3905 ))),
3906 }
3907}
3908
3909pub fn sync_aniso_contrasts_from_metadata(
3916 spec: &mut TermCollectionSpec,
3917 design: &SmoothDesign,
3918) {
3919 for (term_idx, term) in design.terms.iter().enumerate() {
3920 let meta_aniso = match &term.metadata {
3921 BasisMetadata::Matern {
3922 aniso_log_scales, ..
3923 } => aniso_log_scales.clone(),
3924 BasisMetadata::Duchon {
3925 aniso_log_scales, ..
3926 } => aniso_log_scales.clone(),
3927 _ => None,
3928 };
3929 if let Some(eta) = meta_aniso
3930 && eta.len() > 1
3931 {
3932 set_spatial_aniso_log_scales(spec, term_idx, eta).ok();
3933 }
3934 }
3935}
3936
3937#[derive(Debug, Clone)]
3938pub struct SpatialLengthScaleOptimizationOptions {
3939 pub enabled: bool,
3943 pub max_outer_iter: usize,
3945 pub rel_tol: f64,
3947 pub log_step: f64,
3949 pub min_length_scale: f64,
3951 pub max_length_scale: f64,
3953 pub pilot_subsample_threshold: usize,
3966}
3967
3968impl Default for SpatialLengthScaleOptimizationOptions {
3969 fn default() -> Self {
3970 Self {
3971 enabled: true,
3972 max_outer_iter: 80,
3973 rel_tol: 1e-4,
3974 log_step: std::f64::consts::LN_2,
3975 min_length_scale: 1e-3,
3976 max_length_scale: 1e3,
3977 pilot_subsample_threshold: 10_000,
3978 }
3979 }
3980}
3981
3982impl SpatialLengthScaleOptimizationOptions {
3983 pub fn validate(&self) -> Result<(), String> {
4001 if !self.min_length_scale.is_finite() || self.min_length_scale <= 0.0 {
4002 return Err(SmoothError::invalid_config(format!(
4003 "SpatialLengthScaleOptimizationOptions::min_length_scale must be > 0 and finite, got {}",
4004 self.min_length_scale
4005 ))
4006 .into());
4007 }
4008 if !self.max_length_scale.is_finite() || self.max_length_scale <= 0.0 {
4009 return Err(SmoothError::invalid_config(format!(
4010 "SpatialLengthScaleOptimizationOptions::max_length_scale must be > 0 and finite, got {}",
4011 self.max_length_scale
4012 ))
4013 .into());
4014 }
4015 if self.min_length_scale >= self.max_length_scale {
4016 return Err(SmoothError::invalid_config(format!(
4017 "SpatialLengthScaleOptimizationOptions requires min_length_scale < max_length_scale, got min={} max={}",
4018 self.min_length_scale, self.max_length_scale
4019 ))
4020 .into());
4021 }
4022 if !self.rel_tol.is_finite() || self.rel_tol <= 0.0 {
4023 return Err(SmoothError::invalid_config(format!(
4024 "SpatialLengthScaleOptimizationOptions::rel_tol must be > 0 and finite, got {}",
4025 self.rel_tol
4026 ))
4027 .into());
4028 }
4029 if !self.log_step.is_finite() || self.log_step <= 0.0 {
4030 return Err(SmoothError::invalid_config(format!(
4031 "SpatialLengthScaleOptimizationOptions::log_step must be > 0 and finite, got {}",
4032 self.log_step
4033 ))
4034 .into());
4035 }
4036 Ok(())
4037 }
4038}
4039
4040#[derive(Debug, Clone)]
4041pub struct RandomEffectBlock {
4042 pub name: String,
4043 pub group_ids: Vec<Option<usize>>,
4046 pub num_groups: usize,
4047 pub kept_levels: Vec<u64>,
4048}
4049
4050pub const BLOCK_SPARSE_ZERO_EPS: f64 = 1e-12;
4051
4052pub const BLOCK_SPARSE_MAX_DENSITY: f64 = 0.20;
4053
4054pub fn blocks_have_intrinsic_sparse_structure(blocks: &[DesignBlock]) -> bool {
4055 blocks
4056 .iter()
4057 .any(|block| matches!(block, DesignBlock::Sparse(_) | DesignBlock::RandomEffect(_)))
4058}
4059
4060pub fn sparse_compatible_block_nnz(block: &DesignBlock) -> Option<usize> {
4061 match block {
4062 DesignBlock::Intercept(n) => Some(*n),
4063 DesignBlock::RandomEffect(op) => {
4064 Some(op.group_ids.iter().filter(|gid| gid.is_some()).count())
4065 }
4066 DesignBlock::Sparse(sparse) => Some(sparse.val().len()),
4067 DesignBlock::Dense(dense) => dense.as_dense_ref().map(|matrix| {
4068 matrix
4069 .iter()
4070 .filter(|&&value| value.abs() > BLOCK_SPARSE_ZERO_EPS)
4071 .count()
4072 }),
4073 }
4074}
4075
4076pub fn try_build_sparse_design_from_blocks(
4077 blocks: &[DesignBlock],
4078) -> Result<Option<DesignMatrix>, BasisError> {
4079 if blocks.is_empty() {
4080 return Ok(None);
4081 }
4082 let nrows = blocks[0].nrows();
4083 let ncols: usize = blocks.iter().map(DesignBlock::ncols).sum();
4084 if nrows == 0 || ncols == 0 || ncols <= 32 {
4085 return Ok(None);
4086 }
4087
4088 let preserve_sparse_storage = blocks_have_intrinsic_sparse_structure(blocks);
4089 let sparse_nnz_limit = if preserve_sparse_storage {
4090 usize::MAX
4091 } else {
4092 let total_cells = nrows.saturating_mul(ncols);
4093 ((total_cells as f64) * BLOCK_SPARSE_MAX_DENSITY).floor() as usize
4094 };
4095 let mut nnz = 0usize;
4096 for block in blocks {
4097 let block_nnz = if let Some(block_nnz) = sparse_compatible_block_nnz(block) {
4098 block_nnz
4099 } else {
4100 return Ok(None);
4101 };
4102 nnz = nnz.saturating_add(block_nnz);
4103 if nnz > sparse_nnz_limit {
4104 return Ok(None);
4105 }
4106 }
4107
4108 let mut triplets = Vec::<Triplet<usize, usize, f64>>::with_capacity(nnz);
4109 let mut col_offset = 0usize;
4110 for block in blocks {
4111 match block {
4112 DesignBlock::Intercept(n) => {
4113 for row in 0..*n {
4114 triplets.push(Triplet::new(row, col_offset, 1.0));
4115 }
4116 }
4117 DesignBlock::RandomEffect(op) => {
4118 for (row, group_id) in op.group_ids.iter().enumerate() {
4119 if let Some(group) = group_id {
4120 triplets.push(Triplet::new(row, col_offset + group, 1.0));
4121 }
4122 }
4123 }
4124 DesignBlock::Sparse(sparse) => {
4125 let (symbolic, values) = sparse.parts();
4126 let col_ptr = symbolic.col_ptr();
4127 let row_idx = symbolic.row_idx();
4128 for col in 0..sparse.ncols() {
4129 for idx in col_ptr[col]..col_ptr[col + 1] {
4130 let value = values[idx];
4131 if value.abs() > BLOCK_SPARSE_ZERO_EPS {
4132 triplets.push(Triplet::new(row_idx[idx], col_offset + col, value));
4133 }
4134 }
4135 }
4136 }
4137 DesignBlock::Dense(dense) => {
4138 let matrix = dense.as_dense_ref().ok_or_else(|| {
4139 BasisError::InvalidInput(
4140 "sparse-compatible block assembly requires materialized dense blocks"
4141 .to_string(),
4142 )
4143 })?;
4144 for row in 0..matrix.nrows() {
4145 for col in 0..matrix.ncols() {
4146 let value = matrix[[row, col]];
4147 if value.abs() > BLOCK_SPARSE_ZERO_EPS {
4148 triplets.push(Triplet::new(row, col_offset + col, value));
4149 }
4150 }
4151 }
4152 }
4153 }
4154 col_offset += block.ncols();
4155 }
4156
4157 let sparse = SparseColMat::try_new_from_triplets(nrows, ncols, &triplets).map_err(|_| {
4158 BasisError::SparseCreation("failed to assemble sparse term-collection design".to_string())
4159 })?;
4160 Ok(Some(DesignMatrix::Sparse(
4161 gam_linalg::matrix::SparseDesignMatrix::new(sparse),
4162 )))
4163}
4164
4165pub fn assemble_term_collection_design_matrix(
4166 blocks: Vec<DesignBlock>,
4167) -> Result<DesignMatrix, BasisError> {
4168 if let Some(sparse) = try_build_sparse_design_from_blocks(&blocks)? {
4169 return Ok(sparse);
4170 }
4171 let block_op = BlockDesignOperator::new(blocks).map_err(|e| {
4172 BasisError::InvalidInput(format!("failed to build block design operator: {e}"))
4173 })?;
4174 Ok(DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(
4175 Arc::new(block_op),
4176 )))
4177}
4178
4179pub fn select_columns(data: ArrayView2<'_, f64>, cols: &[usize]) -> Result<Array2<f64>, BasisError> {
4180 let n = data.nrows();
4181 let p = data.ncols();
4182 for &c in cols {
4183 if c >= p {
4184 crate::bail_dim_basis!("feature column {c} is out of bounds for data with {p} columns");
4185 }
4186 }
4187 let mut out = Array2::<f64>::zeros((n, cols.len()));
4188 for (j, &c) in cols.iter().enumerate() {
4189 out.column_mut(j).assign(&data.column(c));
4190 }
4191 Ok(out)
4192}
4193
4194pub fn nonfinite_value_label(value: f64) -> &'static str {
4195 if value.is_nan() {
4196 "NaN"
4197 } else if value.is_sign_positive() {
4198 "+Inf"
4199 } else {
4200 "-Inf"
4201 }
4202}
4203
4204pub fn validate_term_feature_column_finite(
4205 data: ArrayView2<'_, f64>,
4206 term_kind: &str,
4207 term_name: &str,
4208 feature_col: usize,
4209) -> Result<(), BasisError> {
4210 let p = data.ncols();
4211 if feature_col >= p {
4212 crate::bail_dim_basis!(
4213 "{term_kind} term '{term_name}' feature column {feature_col} out of bounds for {p} columns"
4214 );
4215 }
4216 for (row, &value) in data.column(feature_col).iter().enumerate() {
4217 if !value.is_finite() {
4218 crate::bail_invalid_basis!(
4219 "{term_kind} term '{term_name}' feature column {feature_col} row {row} contains non-finite value {}",
4220 nonfinite_value_label(value)
4221 );
4222 }
4223 }
4224 Ok(())
4225}
4226
4227pub fn validate_smooth_terms_finite_inputs(
4228 data: ArrayView2<'_, f64>,
4229 terms: &[SmoothTermSpec],
4230) -> Result<(), BasisError> {
4231 for term in terms {
4232 for feature_col in smooth_term_feature_cols(term) {
4233 validate_term_feature_column_finite(data, "smooth", &term.name, feature_col)?;
4234 }
4235 }
4236 Ok(())
4237}
4238
4239pub fn validate_term_collection_finite_inputs(
4240 data: ArrayView2<'_, f64>,
4241 spec: &TermCollectionSpec,
4242) -> Result<(), BasisError> {
4243 for term in &spec.linear_terms {
4244 validate_term_feature_column_finite(data, "linear", &term.name, term.feature_col)?;
4245 }
4246 for term in &spec.random_effect_terms {
4247 validate_term_feature_column_finite(data, "random-effect", &term.name, term.feature_col)?;
4248 }
4249 validate_smooth_terms_finite_inputs(data, &spec.smooth_terms)
4250}
4251
4252#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord)]
4253pub struct JointSpatialCenterGroupKey {
4254 feature_cols: Vec<usize>,
4255 strategy_kind: CenterStrategyKind,
4256 strategy_aux: usize,
4257 requested_num_centers: usize,
4258 input_scale_bits: Option<Vec<u64>>,
4259}
4260
4261pub fn spatial_term_min_center_count(term: &SmoothTermSpec) -> usize {
4262 match &term.basis {
4263 SmoothBasisSpec::ThinPlate { feature_cols, .. } => feature_cols.len() + 1,
4264 SmoothBasisSpec::Duchon {
4265 feature_cols, spec, ..
4266 } => match spec.nullspace_order {
4267 crate::basis::DuchonNullspaceOrder::Zero => 1,
4268 crate::basis::DuchonNullspaceOrder::Linear => feature_cols.len() + 1,
4269 crate::basis::DuchonNullspaceOrder::Degree(degree) => {
4270 crate::basis::duchon_nullspace_dimension(feature_cols.len(), degree)
4271 }
4272 },
4273 SmoothBasisSpec::Matern { .. } => 1,
4274 _ => 1,
4275 }
4276}
4277
4278pub fn spatial_term_group_key(term: &SmoothTermSpec) -> Option<JointSpatialCenterGroupKey> {
4279 let (feature_cols, strategy, input_scales) = match &term.basis {
4280 SmoothBasisSpec::ThinPlate {
4281 feature_cols,
4282 spec,
4283 input_scales,
4284 } => (feature_cols, &spec.center_strategy, input_scales.as_ref()),
4285 SmoothBasisSpec::Matern {
4286 feature_cols,
4287 spec,
4288 input_scales,
4289 } => (feature_cols, &spec.center_strategy, input_scales.as_ref()),
4290 SmoothBasisSpec::Duchon {
4291 feature_cols,
4292 spec,
4293 input_scales,
4294 } => (feature_cols, &spec.center_strategy, input_scales.as_ref()),
4295 _ => return None,
4296 };
4297 let strategy_kind = center_strategy_kind(strategy);
4298 let strategy_aux = match strategy {
4299 CenterStrategy::Auto(inner) => match inner.as_ref() {
4300 CenterStrategy::KMeans { max_iter, .. } => *max_iter,
4301 CenterStrategy::UniformGrid { points_per_dim } => *points_per_dim,
4302 _ => 0,
4303 },
4304 CenterStrategy::KMeans { max_iter, .. } => *max_iter,
4305 CenterStrategy::UniformGrid { points_per_dim } => *points_per_dim,
4306 _ => 0,
4307 };
4308 Some(JointSpatialCenterGroupKey {
4309 feature_cols: feature_cols.clone(),
4310 strategy_kind,
4311 strategy_aux,
4312 requested_num_centers: center_strategy_num_centers(strategy)?,
4313 input_scale_bits: input_scales
4314 .map(|values| values.iter().map(|value| value.to_bits()).collect()),
4315 })
4316}
4317
4318pub fn spatial_term_center_strategy(term: &SmoothTermSpec) -> Option<&CenterStrategy> {
4319 match &term.basis {
4320 SmoothBasisSpec::ThinPlate { spec, .. } => Some(&spec.center_strategy),
4321 SmoothBasisSpec::Matern { spec, .. } => Some(&spec.center_strategy),
4322 SmoothBasisSpec::Duchon { spec, .. } => Some(&spec.center_strategy),
4323 _ => None,
4324 }
4325}
4326
4327pub fn set_spatial_term_centers(
4328 term: &mut SmoothTermSpec,
4329 centers: Array2<f64>,
4330) -> Result<(), BasisError> {
4331 match &mut term.basis {
4332 SmoothBasisSpec::ThinPlate { spec, .. } => {
4333 spec.center_strategy = CenterStrategy::UserProvided(centers);
4334 Ok(())
4335 }
4336 SmoothBasisSpec::Matern { spec, .. } => {
4337 spec.center_strategy = CenterStrategy::UserProvided(centers);
4338 Ok(())
4339 }
4340 SmoothBasisSpec::Duchon { spec, .. } => {
4341 spec.center_strategy = CenterStrategy::UserProvided(centers);
4342 Ok(())
4343 }
4344 _ => Err(BasisError::InvalidInput(format!(
4345 "term '{}' does not support spatial center planning",
4346 term.name
4347 ))),
4348 }
4349}
4350
4351pub fn standardized_spatial_term_data(
4352 data: ArrayView2<'_, f64>,
4353 term: &SmoothTermSpec,
4354) -> Result<Array2<f64>, BasisError> {
4355 let (feature_cols, input_scales) = match &term.basis {
4356 SmoothBasisSpec::ThinPlate {
4357 feature_cols,
4358 input_scales,
4359 ..
4360 }
4361 | SmoothBasisSpec::Matern {
4362 feature_cols,
4363 input_scales,
4364 ..
4365 }
4366 | SmoothBasisSpec::Duchon {
4367 feature_cols,
4368 input_scales,
4369 ..
4370 } => (feature_cols, input_scales.as_ref()),
4371 _ => {
4372 crate::bail_invalid_basis!("term '{}' is not a spatial smooth", term.name);
4373 }
4374 };
4375 let mut x = select_columns(data, feature_cols)?;
4376 if let Some(scales) = input_scales {
4377 apply_input_standardization(&mut x, scales);
4378 } else if let Some(scales) = compute_spatial_input_scales(x.view()) {
4379 apply_input_standardization(&mut x, &scales);
4380 }
4381 Ok(x)
4382}
4383
4384pub fn plan_joint_spatial_centers_for_term_blocks(
4385 data: ArrayView2<'_, f64>,
4386 term_blocks: &[Vec<SmoothTermSpec>],
4387) -> Result<Vec<Vec<SmoothTermSpec>>, BasisError> {
4388 let mut planned_blocks = term_blocks.to_vec();
4389 let n = data.nrows();
4390 let mut groups: BTreeMap<JointSpatialCenterGroupKey, Vec<(usize, usize)>> = BTreeMap::new();
4391
4392 for (block_idx, terms) in planned_blocks.iter().enumerate() {
4393 for (term_idx, term) in terms.iter().enumerate() {
4394 let Some(strategy) = spatial_term_center_strategy(term) else {
4395 continue;
4396 };
4397 if !center_strategy_is_auto(strategy) {
4398 continue;
4399 }
4400 let Some(group_key) = spatial_term_group_key(term) else {
4401 continue;
4402 };
4403 if !matches!(
4404 group_key.strategy_kind,
4405 CenterStrategyKind::EqualMass
4406 | CenterStrategyKind::EqualMassCovarRepresentative
4407 | CenterStrategyKind::FarthestPoint
4408 | CenterStrategyKind::KMeans
4409 ) {
4410 continue;
4411 }
4412 if center_strategy_num_centers(strategy).is_none() {
4413 continue;
4414 }
4415 groups
4416 .entry(group_key)
4417 .or_default()
4418 .push((block_idx, term_idx));
4419 }
4420 }
4421
4422 for (group_key, members) in groups {
4423 if members.len() < 2 {
4424 continue;
4425 }
4426 let min_required = members
4427 .iter()
4428 .map(|&(block_idx, term_idx)| {
4429 spatial_term_min_center_count(&planned_blocks[block_idx][term_idx])
4430 })
4431 .max()
4432 .unwrap_or(1);
4433 let joint_centers = group_key
4434 .requested_num_centers
4435 .max(min_required)
4436 .min(n.max(1));
4437 let (first_block_idx, first_term_idx) = members[0];
4438 let prototype = &planned_blocks[first_block_idx][first_term_idx];
4439 let standardized = standardized_spatial_term_data(data, prototype)?;
4440 let strategy = spatial_term_center_strategy(prototype).ok_or_else(|| {
4441 BasisError::InvalidInput(format!(
4442 "term '{}' lost its spatial center strategy during joint planning",
4443 prototype.name
4444 ))
4445 })?;
4446 let joint_strategy = center_strategy_with_num_centers(strategy, joint_centers)?;
4447 let shared_centers = select_centers_by_strategy(standardized.view(), &joint_strategy)?;
4448 log::info!(
4449 "sharing {} spatial centers across {} smooth terms over columns {:?} (requested {} centers)",
4450 shared_centers.nrows(),
4451 members.len(),
4452 group_key.feature_cols,
4453 group_key.requested_num_centers,
4454 );
4455 for (block_idx, term_idx) in members {
4456 set_spatial_term_centers(
4457 &mut planned_blocks[block_idx][term_idx],
4458 shared_centers.clone(),
4459 )?;
4460 }
4461 }
4462
4463 for block in planned_blocks.iter_mut() {
4470 for term in block.iter_mut() {
4471 auto_init_length_scale_in_place(data, term);
4472 }
4473 }
4474
4475 Ok(planned_blocks)
4476}
4477
4478const AUTO_LENGTH_SCALE_FLOOR: f64 = 1e-6;
4481
4482fn feature_columns_max_range(data: ArrayView2<'_, f64>, feature_cols: &[usize]) -> Option<f64> {
4485 let mut max_range = 0.0_f64;
4486 for &c in feature_cols {
4487 if c >= data.ncols() {
4488 continue;
4489 }
4490 let col = data.column(c);
4491 let mut lo = f64::INFINITY;
4492 let mut hi = f64::NEG_INFINITY;
4493 for &v in col.iter() {
4494 if v.is_finite() {
4495 if v < lo {
4496 lo = v;
4497 }
4498 if v > hi {
4499 hi = v;
4500 }
4501 }
4502 }
4503 if hi > lo {
4504 let r = hi - lo;
4505 if r > max_range {
4506 max_range = r;
4507 }
4508 }
4509 }
4510 if max_range.is_finite() && max_range > 0.0 {
4511 Some(max_range)
4512 } else {
4513 None
4514 }
4515}
4516
4517pub fn auto_initial_length_scale(data: ArrayView2<'_, f64>, feature_cols: &[usize]) -> f64 {
4524 let n = data.nrows();
4525 if n == 0 || feature_cols.is_empty() {
4526 return 1.0;
4527 }
4528 let Some(max_range) = feature_columns_max_range(data, feature_cols) else {
4529 return 1.0;
4530 };
4531 let init = max_range / (n as f64).sqrt();
4532 init.max(AUTO_LENGTH_SCALE_FLOOR).min(max_range)
4533}
4534
4535pub fn auto_initial_length_scale_for_centers(
4558 data: ArrayView2<'_, f64>,
4559 feature_cols: &[usize],
4560 num_centers: usize,
4561) -> f64 {
4562 let n = data.nrows();
4563 if n == 0 || feature_cols.is_empty() {
4564 return 1.0;
4565 }
4566 let Some(max_range) = feature_columns_max_range(data, feature_cols) else {
4567 return 1.0;
4568 };
4569 let resolution_points = n.max(num_centers).max(1) as f64;
4575 let spacing = max_range / resolution_points.sqrt();
4576 spacing.max(AUTO_LENGTH_SCALE_FLOOR).min(max_range)
4577}
4578
4579pub fn auto_initial_length_scale_for_low_rank_centers(
4589 data: ArrayView2<'_, f64>,
4590 feature_cols: &[usize],
4591 num_centers: usize,
4592) -> f64 {
4593 if data.nrows() == 0 || feature_cols.is_empty() {
4594 return 1.0;
4595 }
4596 let Some(max_range) = feature_columns_max_range(data, feature_cols) else {
4597 return 1.0;
4598 };
4599 let resolution_points = num_centers.max(1) as f64;
4600 let spacing = max_range / resolution_points.sqrt();
4601 spacing.max(AUTO_LENGTH_SCALE_FLOOR).min(max_range)
4602}
4603
4604fn center_strategy_requested_count(strategy: &CenterStrategy) -> Option<usize> {
4607 match strategy {
4608 CenterStrategy::Auto(inner) => center_strategy_requested_count(inner),
4609 CenterStrategy::UserProvided(centers) => Some(centers.nrows()),
4610 CenterStrategy::EqualMass { num_centers }
4611 | CenterStrategy::EqualMassCovarRepresentative { num_centers }
4612 | CenterStrategy::FarthestPoint { num_centers }
4613 | CenterStrategy::KMeans { num_centers, .. } => Some(*num_centers),
4614 CenterStrategy::UniformGrid { .. } => None,
4615 }
4616}
4617
4618pub fn auto_init_length_scale_in_place(data: ArrayView2<'_, f64>, term: &mut SmoothTermSpec) {
4622 auto_init_length_scale_in_basis(data, &mut term.basis);
4623}
4624
4625pub fn auto_init_length_scale_in_basis(data: ArrayView2<'_, f64>, basis: &mut SmoothBasisSpec) {
4638 match basis {
4639 SmoothBasisSpec::Matern {
4640 feature_cols, spec, ..
4641 } => {
4642 if spec.length_scale == 0.0 {
4643 spec.length_scale = match center_strategy_requested_count(&spec.center_strategy) {
4652 Some(k) => auto_initial_length_scale_for_centers(data, feature_cols, k),
4653 None => auto_initial_length_scale(data, feature_cols),
4654 };
4655 }
4656 }
4657 SmoothBasisSpec::ThinPlate {
4658 feature_cols, spec, ..
4659 } => {
4660 if spec.length_scale == 0.0 {
4661 spec.length_scale = match center_strategy_requested_count(&spec.center_strategy) {
4662 Some(k) => {
4663 auto_initial_length_scale_for_low_rank_centers(data, feature_cols, k)
4664 }
4665 None => auto_initial_length_scale(data, feature_cols),
4666 };
4667 }
4668 }
4669 SmoothBasisSpec::ByVariable { inner, .. }
4670 | SmoothBasisSpec::FactorSumToZero { inner, .. } => {
4671 auto_init_length_scale_in_basis(data, inner);
4672 }
4673 SmoothBasisSpec::BySmooth { smooth, .. } => {
4674 auto_init_length_scale_in_basis(data, smooth);
4675 }
4676 _ => {}
4677 }
4678}
4679
4680impl LinearFitConditioning {
4681 pub fn from_columns(design: &TermCollectionDesign, selected_cols: &[usize]) -> Self {
4682 const SCALE_EPS: f64 = 1e-12;
4683 let n = design.design.nrows();
4684 let p = design.design.ncols();
4685 let mut columns = Vec::with_capacity(selected_cols.len());
4686 if n == 0 || selected_cols.is_empty() {
4687 return Self {
4688 intercept_idx: design.intercept_range.start,
4689 columns,
4690 };
4691 }
4692 let chunk_rows = gam_linalg::utils::row_chunk_for_byte_budget(n, p);
4693 let mut sums = vec![0.0_f64; selected_cols.len()];
4699 for start in (0..n).step_by(chunk_rows) {
4700 let end = (start + chunk_rows).min(n);
4701 let chunk = design
4702 .design
4703 .try_row_chunk(start..end)
4704 .expect("LinearFitConditioning::from_columns row chunk failed");
4705 for (k, &col_idx) in selected_cols.iter().enumerate() {
4706 let column = chunk.column(col_idx);
4707 for &v in column.iter() {
4708 sums[k] += v;
4709 }
4710 }
4711 }
4712 let inv_n = 1.0_f64 / n as f64;
4713 let means: Vec<f64> = sums.iter().map(|&s| s * inv_n).collect();
4714 let mut sq_devs = vec![0.0_f64; selected_cols.len()];
4715 for start in (0..n).step_by(chunk_rows) {
4716 let end = (start + chunk_rows).min(n);
4717 let chunk = design
4718 .design
4719 .try_row_chunk(start..end)
4720 .expect("LinearFitConditioning::from_columns row chunk failed");
4721 for (k, &col_idx) in selected_cols.iter().enumerate() {
4722 let mean_k = means[k];
4723 let column = chunk.column(col_idx);
4724 for &v in column.iter() {
4725 let d = v - mean_k;
4726 sq_devs[k] += d * d;
4727 }
4728 }
4729 }
4730 for (k, &col_idx) in selected_cols.iter().enumerate() {
4731 let mean = means[k];
4732 let var = sq_devs[k] * inv_n;
4733 let (mean, scale) = if var.is_finite() && var > SCALE_EPS * SCALE_EPS {
4734 (mean, var.sqrt())
4735 } else {
4736 (0.0, 1.0)
4739 };
4740 columns.push(LinearColumnConditioning {
4741 col_idx,
4742 mean,
4743 scale,
4744 });
4745 }
4746 Self {
4747 intercept_idx: design.intercept_range.start,
4748 columns,
4749 }
4750 }
4751
4752 pub fn apply_to_design(&self, design: &Array2<f64>) -> Array2<f64> {
4753 let mut out = design.clone();
4754 for col in &self.columns {
4755 {
4756 let mut dst = out.column_mut(col.col_idx);
4757 dst -= col.mean;
4758 }
4759 if col.scale != 1.0 {
4760 out.column_mut(col.col_idx).mapv_inplace(|v| v / col.scale);
4761 }
4762 }
4763 out
4764 }
4765
4766 fn transform_matrix_columnswith_a(&self, mat: &Array2<f64>) -> Array2<f64> {
4767 let mut out = mat.clone();
4768 let intercept = self.intercept_idx;
4769 for col in &self.columns {
4770 let intercept_col = out.column(intercept).to_owned();
4771 let mut target = out.column_mut(col.col_idx);
4772 target -= &(intercept_col * col.mean);
4773 if col.scale != 1.0 {
4774 target.mapv_inplace(|v| v / col.scale);
4775 }
4776 }
4777 out
4778 }
4779
4780 fn transform_matrixrowswith_a_transpose(&self, mat: &Array2<f64>) -> Array2<f64> {
4781 let mut out = mat.clone();
4782 let intercept = self.intercept_idx;
4783 for col in &self.columns {
4784 let interceptrow = out.row(intercept).to_owned();
4785 let mut target = out.row_mut(col.col_idx);
4786 target -= &(interceptrow * col.mean);
4787 if col.scale != 1.0 {
4788 target.mapv_inplace(|v| v / col.scale);
4789 }
4790 }
4791 out
4792 }
4793
4794 fn left_multiply_by_m_inv_transpose(&self, mat_internal: &Array2<f64>) -> Array2<f64> {
4799 let mut out = mat_internal.clone();
4800 let intercept = self.intercept_idx;
4801 let interceptrow_snapshot = mat_internal.row(intercept).to_owned();
4802 for col in &self.columns {
4803 if col.scale != 1.0 {
4804 out.row_mut(col.col_idx).mapv_inplace(|v| v * col.scale);
4805 }
4806 if col.mean != 0.0 {
4807 let mut target = out.row_mut(col.col_idx);
4808 target += &(&interceptrow_snapshot * col.mean);
4809 }
4810 }
4811 out
4812 }
4813
4814 fn right_multiply_by_m_inv(&self, mat_internal: &Array2<f64>) -> Array2<f64> {
4817 let mut out = mat_internal.clone();
4818 let intercept = self.intercept_idx;
4819 let intercept_col_snapshot = mat_internal.column(intercept).to_owned();
4820 for col in &self.columns {
4821 if col.scale != 1.0 {
4822 out.column_mut(col.col_idx).mapv_inplace(|v| v * col.scale);
4823 }
4824 if col.mean != 0.0 {
4825 let mut target = out.column_mut(col.col_idx);
4826 target += &(&intercept_col_snapshot * col.mean);
4827 }
4828 }
4829 out
4830 }
4831
4832 pub fn transform_blockwise_penalties_to_internal(
4839 &self,
4840 penalties: &[BlockwisePenalty],
4841 p: usize,
4842 ) -> Vec<crate::penalty_spec::PenaltySpec> {
4843 let conditioning_cols: std::collections::HashSet<usize> =
4844 self.columns.iter().map(|c| c.col_idx).collect();
4845 penalties
4846 .iter()
4847 .map(|bp| {
4848 let overlaps =
4849 (bp.col_range.start..bp.col_range.end).any(|j| conditioning_cols.contains(&j));
4850 if overlaps {
4851 let global = bp.to_global(p);
4854 let right = self.transform_matrix_columnswith_a(&global);
4855 let transformed = self.transform_matrixrowswith_a_transpose(&right);
4856 crate::penalty_spec::PenaltySpec::Dense(transformed)
4857 } else {
4858 crate::penalty_spec::PenaltySpec::from_blockwise(bp.clone())
4861 }
4862 })
4863 .collect()
4864 }
4865
4866 pub fn backtransform_beta(&self, beta_internal: &Array1<f64>) -> Array1<f64> {
4867 let mut beta = beta_internal.clone();
4868 let intercept = self.intercept_idx;
4869 for col in &self.columns {
4870 beta[intercept] -= beta_internal[col.col_idx] * col.mean / col.scale;
4871 beta[col.col_idx] = beta_internal[col.col_idx] / col.scale;
4872 }
4873 beta
4874 }
4875
4876 pub fn transform_penalized_hessian_to_original(&self, h_internal: &Array2<f64>) -> Array2<f64> {
4879 let right = self.right_multiply_by_m_inv(h_internal);
4880 self.left_multiply_by_m_inv_transpose(&right)
4881 }
4882
4883 pub fn internal_bounds_for(&self, col_idx: usize, min: f64, max: f64) -> (f64, f64) {
4884 if let Some(col) = self.columns.iter().find(|c| c.col_idx == col_idx) {
4885 (min * col.scale, max * col.scale)
4886 } else {
4887 (min, max)
4888 }
4889 }
4890}
4891
4892pub fn freeze_raw_spatial_metadata(metadata: BasisMetadata, raw_cols: usize) -> BasisMetadata {
4893 match metadata {
4894 BasisMetadata::ThinPlate {
4895 centers,
4896 length_scale,
4897 periodic,
4898 identifiability_transform: None,
4899 input_scales,
4900 radial_reparam,
4901 } => BasisMetadata::ThinPlate {
4902 centers,
4903 length_scale,
4904 periodic,
4905 identifiability_transform: Some(Array2::eye(raw_cols)),
4906 input_scales,
4907 radial_reparam,
4908 },
4909 BasisMetadata::Duchon {
4910 centers,
4911 length_scale,
4912 periodic,
4913 power,
4914 nullspace_order,
4915 identifiability_transform: None,
4916 input_scales,
4917 aniso_log_scales,
4918 operator_collocation_points,
4919 radial_reparam,
4920 } => BasisMetadata::Duchon {
4921 centers,
4922 length_scale,
4923 periodic,
4924 power,
4925 nullspace_order,
4926 identifiability_transform: Some(Array2::eye(raw_cols)),
4927 input_scales,
4928 aniso_log_scales,
4929 operator_collocation_points,
4930 radial_reparam,
4931 },
4932 other => other,
4933 }
4934}
4935
4936pub fn matern_operator_penalty_triplet_from_metadata(
4937 metadata: &BasisMetadata,
4938) -> Result<(Vec<Array2<f64>>, Vec<usize>, Vec<PenaltyInfo>), BasisError> {
4939 let BasisMetadata::Matern {
4940 centers,
4941 length_scale,
4942 periodic,
4943 nu,
4944 include_intercept,
4945 identifiability_transform,
4946 aniso_log_scales,
4947 input_scales,
4948 ..
4949 } = metadata
4950 else {
4951 crate::bail_invalid_basis!("Matérn operator penalties require Matérn metadata");
4952 };
4953 let penalty_length_scale = match input_scales.as_deref() {
4965 Some(scales) => compensate_length_scale_for_standardization(*length_scale, scales),
4966 None => *length_scale,
4967 };
4968 matern_operator_penalty_triplet_at_length_scale(
4969 centers.view(),
4970 periodic.as_deref(),
4971 identifiability_transform.as_ref(),
4972 *nu,
4973 *include_intercept,
4974 aniso_log_scales.as_deref(),
4975 penalty_length_scale,
4976 )
4977}
4978
4979pub fn matern_operator_penalty_triplet_at_length_scale(
4997 centers: ArrayView2<'_, f64>,
4998 periodic: Option<&[Option<f64>]>,
4999 identifiability_transform: Option<&Array2<f64>>,
5000 nu: crate::basis::MaternNu,
5001 include_intercept: bool,
5002 aniso_log_scales: Option<&[f64]>,
5003 effective_length_scale: f64,
5004) -> Result<(Vec<Array2<f64>>, Vec<usize>, Vec<PenaltyInfo>), BasisError> {
5005 let penalty_centers = crate::basis::expand_periodic_centers(¢ers.to_owned(), periodic)?;
5006 let ops = build_matern_collocation_operator_matrices(
5007 penalty_centers.view(),
5008 None,
5009 effective_length_scale,
5010 nu,
5011 include_intercept,
5012 identifiability_transform.map(|z| z.view()),
5013 aniso_log_scales,
5014 )?;
5015 const ORDER_EPS: f64 = 1e-9;
5023 let d = penalty_centers.ncols();
5024 let m = nu.half_integer_value() + 0.5 * d as f64;
5025 let mut candidates = Vec::with_capacity(3);
5026 for (raw, source, min_order) in [
5027 (ops.d0.t().dot(&ops.d0), PenaltySource::OperatorMass, 0.0),
5028 (ops.d1.t().dot(&ops.d1), PenaltySource::OperatorTension, 1.0),
5029 (
5030 ops.d2.t().dot(&ops.d2),
5031 PenaltySource::OperatorStiffness,
5032 2.0,
5033 ),
5034 ] {
5035 if min_order > 0.0 && m <= min_order + ORDER_EPS {
5036 continue;
5037 }
5038 let sym = (&raw + &raw.t()) * 0.5;
5039 let (matrix, normalization_scale) = normalize_penalty_in_constrained_space(&sym);
5040 candidates.push(PenaltyCandidate {
5041 matrix,
5042 nullspace_dim_hint: 0,
5043 source,
5044 normalization_scale,
5045 kronecker_factors: None,
5046 op: None,
5047 });
5048 }
5049 filter_active_penalty_candidates(candidates)
5050}
5051
5052pub fn normalize_penalty_in_constrained_space(matrix: &Array2<f64>) -> (Array2<f64>, f64) {
5053 let matrix = (matrix + &matrix.t().to_owned()) * 0.5;
5058 let matrix = crate::basis::project_penalty_to_psd_cone(&matrix);
5060 let c = matrix.iter().map(|v| v * v).sum::<f64>().sqrt();
5061 if c.is_finite() && c > 0.0 {
5062 (matrix.mapv(|v| v / c), c)
5063 } else {
5064 (matrix, 1.0)
5065 }
5066}
5067
5068pub fn tensor_product_design_from_sparse_marginals(
5069 marginal_sparse: &[&SparseColMat<usize, f64>],
5070) -> Result<SparseColMat<usize, f64>, BasisError> {
5071 if marginal_sparse.is_empty() {
5072 crate::bail_invalid_basis!("TensorBSpline requires at least one marginal basis");
5073 }
5074 let n = marginal_sparse[0].nrows();
5075 for (i, m) in marginal_sparse.iter().enumerate().skip(1) {
5076 if m.nrows() != n {
5077 crate::bail_dim_basis!(
5078 "tensor sparse marginal row mismatch at dim {i}: expected {n}, got {}",
5079 m.nrows()
5080 );
5081 }
5082 }
5083 let dims: Vec<usize> = marginal_sparse.iter().map(|m| m.ncols()).collect();
5084 let total_cols = dims.iter().try_fold(1usize, |acc, &q| {
5085 acc.checked_mul(q)
5086 .ok_or_else(|| BasisError::DimensionMismatch("tensor basis too large".to_string()))
5087 })?;
5088 let mut strides = vec![1usize; dims.len()];
5089 for d in (0..dims.len().saturating_sub(1)).rev() {
5090 strides[d] = strides[d + 1]
5091 .checked_mul(dims[d + 1])
5092 .ok_or_else(|| BasisError::DimensionMismatch("tensor basis too large".to_string()))?;
5093 }
5094
5095 use faer::sparse::SparseRowMat;
5096 let csrs: Vec<SparseRowMat<usize, f64>> = marginal_sparse
5097 .iter()
5098 .enumerate()
5099 .map(|(d, m)| {
5100 m.as_ref().to_row_major().map_err(|e| {
5101 BasisError::SparseCreation(format!(
5102 "tensor sparse marginal {d} CSR conversion failed: {e:?}"
5103 ))
5104 })
5105 })
5106 .collect::<Result<Vec<_>, _>>()?;
5107 let row_ptrs: Vec<&[usize]> = csrs.iter().map(|c| c.symbolic().row_ptr()).collect();
5108 let col_idxs: Vec<&[usize]> = csrs.iter().map(|c| c.symbolic().col_idx()).collect();
5109 let vals: Vec<&[f64]> = csrs.iter().map(|c| c.val()).collect();
5110
5111 use rayon::prelude::*;
5112 const CHUNK: usize = 1024;
5113 let num_chunks = n.div_ceil(CHUNK);
5114 let per_chunk: Vec<Vec<Triplet<usize, usize, f64>>> = (0..num_chunks)
5115 .into_par_iter()
5116 .map(|chunk_idx| {
5117 let row_start = chunk_idx * CHUNK;
5118 let row_end = (row_start + CHUNK).min(n);
5119 let mut chunk_triplets = Vec::<Triplet<usize, usize, f64>>::new();
5120 let mut cur_cols = Vec::<usize>::with_capacity(64);
5121 let mut cur_vals = Vec::<f64>::with_capacity(64);
5122 let mut next_cols = Vec::<usize>::with_capacity(64);
5123 let mut next_vals = Vec::<f64>::with_capacity(64);
5124 for i in row_start..row_end {
5125 cur_cols.clear();
5126 cur_vals.clear();
5127 cur_cols.push(0);
5128 cur_vals.push(1.0);
5129 let mut row_is_zero = false;
5130 for d in 0..dims.len() {
5131 let row_start_d = row_ptrs[d][i];
5132 let row_end_d = row_ptrs[d][i + 1];
5133 if row_start_d == row_end_d {
5134 row_is_zero = true;
5135 break;
5136 }
5137 let stride = strides[d];
5138 next_cols.clear();
5139 next_vals.clear();
5140 next_cols.reserve(cur_cols.len() * (row_end_d - row_start_d));
5141 next_vals.reserve(cur_vals.len() * (row_end_d - row_start_d));
5142 for (&prev_col, &prev_val) in cur_cols.iter().zip(cur_vals.iter()) {
5143 for ptr in row_start_d..row_end_d {
5144 let cj = col_idxs[d][ptr];
5145 let vj = vals[d][ptr];
5146 next_cols.push(prev_col + cj * stride);
5147 next_vals.push(prev_val * vj);
5148 }
5149 }
5150 std::mem::swap(&mut cur_cols, &mut next_cols);
5151 std::mem::swap(&mut cur_vals, &mut next_vals);
5152 }
5153 if row_is_zero {
5154 continue;
5155 }
5156 for (&col, &val) in cur_cols.iter().zip(cur_vals.iter()) {
5157 chunk_triplets.push(Triplet::new(i, col, val));
5158 }
5159 }
5160 chunk_triplets
5161 })
5162 .collect();
5163 let total_nnz: usize = per_chunk.iter().map(Vec::len).sum();
5164 let mut triplets = Vec::<Triplet<usize, usize, f64>>::with_capacity(total_nnz);
5165 for chunk in per_chunk {
5166 triplets.extend(chunk);
5167 }
5168 SparseColMat::try_new_from_triplets(n, total_cols, &triplets).map_err(|e| {
5169 BasisError::SparseCreation(format!(
5170 "failed to assemble sparse tensor product design: {e:?}"
5171 ))
5172 })
5173}
5174
5175pub fn dense_local_margin_to_sparse(
5176 dense: &Array2<f64>,
5177) -> Result<SparseColMat<usize, f64>, BasisError> {
5178 let expected_row_nnz = dense.ncols().min(4);
5179 let mut triplets =
5180 Vec::<Triplet<usize, usize, f64>>::with_capacity(dense.nrows() * expected_row_nnz);
5181 for ((row, col), &value) in dense.indexed_iter() {
5182 if value != 0.0 {
5183 triplets.push(Triplet::new(row, col, value));
5184 }
5185 }
5186 SparseColMat::try_new_from_triplets(dense.nrows(), dense.ncols(), &triplets).map_err(|e| {
5187 BasisError::SparseCreation(format!(
5188 "failed to convert tensor marginal design to sparse form: {e:?}"
5189 ))
5190 })
5191}
5192
5193pub struct TensorMarginRangeNullProjectors {
5194 range: Array2<f64>,
5195 null: Array2<f64>,
5196}
5197
5198pub fn projector_from_columns(columns: &Array2<f64>, indices: &[usize]) -> Array2<f64> {
5199 if indices.is_empty() {
5200 return Array2::<f64>::zeros((columns.nrows(), columns.nrows()));
5201 }
5202 let basis = columns.select(Axis(1), indices);
5203 basis.dot(&basis.t())
5204}
5205
5206pub fn tensor_margin_range_null_projectors(
5207 normalized_marginal_penalties: &[(Array2<f64>, f64)],
5208) -> Result<Vec<TensorMarginRangeNullProjectors>, BasisError> {
5209 normalized_marginal_penalties
5210 .iter()
5211 .enumerate()
5212 .map(|(dim, (penalty, _))| {
5213 let analysis = crate::basis::analyze_penalty_block(penalty)?;
5214 if analysis.rank == 0 {
5215 crate::bail_invalid_basis!(
5216 "t2 separable tensor penalty margin {dim} has rank-zero penalty; \
5217 cannot split penalized and null subspaces"
5218 );
5219 }
5220 let mut range_idx = Vec::<usize>::new();
5221 let mut null_idx = Vec::<usize>::new();
5222 for (idx, &ev) in analysis.eigenvalues.iter().enumerate() {
5223 if ev > analysis.tol {
5224 range_idx.push(idx);
5225 } else {
5226 null_idx.push(idx);
5227 }
5228 }
5229 Ok(TensorMarginRangeNullProjectors {
5230 range: projector_from_columns(&analysis.eigenvectors, &range_idx),
5231 null: projector_from_columns(&analysis.eigenvectors, &null_idx),
5232 })
5233 })
5234 .collect()
5235}
5236
5237pub fn build_tensor_bspline_basis(
5238 data: ArrayView2<'_, f64>,
5239 feature_cols: &[usize],
5240 spec: &TensorBSplineSpec,
5241) -> Result<BasisBuildResult, BasisError> {
5242 if feature_cols.is_empty() {
5243 crate::bail_invalid_basis!("TensorBSpline requires at least one feature column");
5244 }
5245 if feature_cols.len() != spec.marginalspecs.len() {
5246 crate::bail_dim_basis!(
5247 "TensorBSpline feature/spec mismatch: feature_cols={}, marginalspecs={}",
5248 feature_cols.len(),
5249 spec.marginalspecs.len()
5250 );
5251 }
5252 if !spec.periods.is_empty() && spec.periods.len() != feature_cols.len() {
5253 crate::bail_dim_basis!(
5254 "TensorBSpline periods length {} does not match feature count {}",
5255 spec.periods.len(),
5256 feature_cols.len()
5257 );
5258 }
5259 let p = data.ncols();
5260 for &c in feature_cols {
5261 if c >= p {
5262 crate::bail_dim_basis!(
5263 "tensor feature column {c} is out of bounds for data with {p} columns"
5264 );
5265 }
5266 }
5267
5268 let mut marginal_knots = Vec::<Array1<f64>>::with_capacity(feature_cols.len());
5269 let mut marginal_is_cr_flags = Vec::<bool>::with_capacity(feature_cols.len());
5272 let mut marginal_degrees = Vec::<usize>::with_capacity(feature_cols.len());
5273 let mut marginalnum_basis = Vec::<usize>::with_capacity(feature_cols.len());
5274 let mut marginal_penalties = Vec::<Array2<f64>>::with_capacity(feature_cols.len());
5275 let mut marginal_designs = Vec::<Array2<f64>>::with_capacity(feature_cols.len());
5276 let mut marginal_effective_periods = Vec::<Option<f64>>::with_capacity(feature_cols.len());
5284 let mut marginal_sparse =
5292 Vec::<Option<SparseColMat<usize, f64>>>::with_capacity(feature_cols.len());
5293
5294 for (dim, (&col, marginalspec)) in feature_cols
5297 .iter()
5298 .zip(spec.marginalspecs.iter())
5299 .enumerate()
5300 {
5301 let mut marginal_unconstrained = marginalspec.clone();
5306 marginal_unconstrained.identifiability = BSplineIdentifiability::None;
5307 let built = build_bspline_basis_1d(data.column(col), &marginal_unconstrained)?;
5308 let (knots, marginal_is_cr) = match built.metadata {
5313 BasisMetadata::BSpline1D { knots, .. } => (knots, false),
5314 BasisMetadata::CubicRegression1D { knots, .. } => (knots, true),
5315 _ => {
5316 crate::bail_invalid_basis!(
5317 "internal TensorBSpline error at dim {dim}: expected BSpline1D or CubicRegression1D metadata"
5318 );
5319 }
5320 };
5321 let metadata_knots = match marginalspec.knotspec {
5322 BSplineKnotSpec::PeriodicUniform {
5323 data_range,
5324 num_basis,
5325 } => Array1::linspace(data_range.0, data_range.1, num_basis),
5326 _ => knots,
5327 };
5328 marginal_knots.push(metadata_knots);
5329 marginal_is_cr_flags.push(marginal_is_cr);
5330 marginal_degrees.push(marginalspec.degree);
5331 marginalnum_basis.push(built.design.ncols());
5332 let dense_marginal = built.design.to_dense();
5337 let sparse_view: Option<SparseColMat<usize, f64>> = match built.design.as_sparse() {
5338 Some(sd) => {
5339 let inner: &SparseColMat<usize, f64> = sd;
5340 Some(inner.clone())
5341 }
5342 None => match marginalspec.knotspec {
5343 BSplineKnotSpec::PeriodicUniform { .. } => {
5344 Some(dense_local_margin_to_sparse(&dense_marginal)?)
5345 }
5346 _ => None,
5347 },
5348 };
5349 marginal_sparse.push(sparse_view);
5350 marginal_designs.push(dense_marginal);
5351 marginal_penalties.push(
5352 built
5353 .penalties
5354 .first()
5355 .ok_or_else(|| {
5356 BasisError::InvalidInput(format!(
5357 "internal TensorBSpline error at dim {dim}: missing marginal penalty"
5358 ))
5359 })?
5360 .clone(),
5361 );
5362 built.nullspace_dims.first().ok_or_else(|| {
5363 BasisError::InvalidInput(format!(
5364 "internal TensorBSpline error at dim {dim}: missing marginal nullspace dim"
5365 ))
5366 })?;
5367 let implied_period = match marginalspec.knotspec {
5375 BSplineKnotSpec::PeriodicUniform { data_range, .. } => {
5376 Some(data_range.1 - data_range.0)
5377 }
5378 _ => spec.periods.get(dim).and_then(|p| *p),
5379 };
5380 marginal_effective_periods.push(implied_period);
5381 }
5382
5383 let total_cols: usize = marginalnum_basis.iter().product();
5384 let mut dense_design = (!matches!(spec.identifiability, TensorBSplineIdentifiability::None))
5385 .then(|| tensor_product_design_from_marginals(&marginal_designs))
5386 .transpose()?;
5387 let mut candidates = Vec::<PenaltyCandidate>::with_capacity(
5388 match spec.penalty_decomposition {
5389 TensorBSplinePenaltyDecomposition::MarginalKroneckerSum => marginal_penalties.len(),
5390 TensorBSplinePenaltyDecomposition::Separable => marginal_penalties.len() * 2,
5391 } + if spec.double_penalty { 1 } else { 0 },
5392 );
5393
5394 let normalized_marginal_penalties: Vec<(Array2<f64>, f64)> = marginal_penalties
5402 .iter()
5403 .map(normalize_penalty_in_constrained_space)
5404 .collect();
5405 let mut kronecker_marginal_penalties =
5406 Vec::<Array2<f64>>::with_capacity(normalized_marginal_penalties.len());
5407
5408 match spec.penalty_decomposition {
5409 TensorBSplinePenaltyDecomposition::MarginalKroneckerSum => {
5410 let mut marginal_kron_sum = Array2::<f64>::zeros((total_cols, total_cols));
5416
5417 for dim in 0..normalized_marginal_penalties.len() {
5418 let mut s_dim = Array2::<f64>::eye(1);
5419 let mut factors = Vec::<Array2<f64>>::with_capacity(marginalnum_basis.len());
5420 for (j, &qj) in marginalnum_basis.iter().enumerate() {
5421 let factor = if j == dim {
5422 normalized_marginal_penalties[j].0.clone()
5423 } else {
5424 Array2::<f64>::eye(qj)
5425 };
5426 factors.push(factor.clone());
5427 s_dim = kronecker_product(&s_dim, &factor);
5428 }
5429 if dim == kronecker_marginal_penalties.len() {
5430 kronecker_marginal_penalties.push(normalized_marginal_penalties[dim].0.clone());
5431 }
5432 marginal_kron_sum += &s_dim;
5433
5434 candidates.push(PenaltyCandidate {
5435 matrix: s_dim,
5436 nullspace_dim_hint: 0,
5437 source: PenaltySource::TensorMarginal { dim },
5438 normalization_scale: normalized_marginal_penalties[dim].1,
5439 kronecker_factors: Some(factors),
5440 op: None,
5441 });
5442 }
5443
5444 if spec.double_penalty
5445 && let Some(shrink) =
5446 crate::basis::build_nullspace_shrinkage_penalty(&marginal_kron_sum)?
5447 {
5448 let (matrix, normalization_scale) =
5449 normalize_penalty_in_constrained_space(&shrink.sym_penalty);
5450 candidates.push(PenaltyCandidate {
5451 matrix,
5452 nullspace_dim_hint: 0,
5453 source: PenaltySource::TensorGlobalRidge,
5454 normalization_scale,
5455 kronecker_factors: None,
5456 op: None,
5457 });
5458 }
5459 }
5460 TensorBSplinePenaltyDecomposition::Separable => {
5461 let projectors = tensor_margin_range_null_projectors(&normalized_marginal_penalties)?;
5462 let n_masks = 1usize.checked_shl(projectors.len() as u32).ok_or_else(|| {
5463 BasisError::InvalidInput(format!(
5464 "t2 separable tensor penalty supports at most {} margins, got {}",
5465 usize::BITS - 1,
5466 projectors.len()
5467 ))
5468 })?;
5469 for mask in 1..n_masks {
5470 let mut matrix = Array2::<f64>::eye(1);
5471 let mut factors = Vec::<Array2<f64>>::with_capacity(projectors.len());
5472 let mut penalized_margins = Vec::<usize>::new();
5473 for (dim, projector) in projectors.iter().enumerate() {
5474 let use_range = ((mask >> dim) & 1) == 1;
5475 let factor = if use_range {
5476 penalized_margins.push(dim);
5477 projector.range.clone()
5478 } else {
5479 projector.null.clone()
5480 };
5481 matrix = kronecker_product(&matrix, &factor);
5482 factors.push(factor);
5483 }
5484 let (matrix, normalization_scale) = normalize_penalty_in_constrained_space(&matrix);
5485 candidates.push(PenaltyCandidate {
5486 matrix,
5487 nullspace_dim_hint: 0,
5488 source: PenaltySource::TensorSeparable { penalized_margins },
5489 normalization_scale,
5490 kronecker_factors: Some(factors),
5491 op: None,
5492 });
5493 }
5494
5495 if spec.double_penalty {
5496 let mut matrix = Array2::<f64>::eye(1);
5497 let mut factors = Vec::<Array2<f64>>::with_capacity(projectors.len());
5498 for projector in &projectors {
5499 matrix = kronecker_product(&matrix, &projector.null);
5500 factors.push(projector.null.clone());
5501 }
5502 let (matrix, normalization_scale) = normalize_penalty_in_constrained_space(&matrix);
5503 candidates.push(PenaltyCandidate {
5504 matrix,
5505 nullspace_dim_hint: 0,
5506 source: PenaltySource::TensorGlobalRidge,
5507 normalization_scale,
5508 kronecker_factors: Some(factors),
5509 op: None,
5510 });
5511 }
5512 }
5513 }
5514
5515 let z_opt = match &spec.identifiability {
5516 TensorBSplineIdentifiability::None => None,
5517 TensorBSplineIdentifiability::SumToZero => {
5518 if total_cols < 2 {
5519 crate::bail_invalid_basis!(
5520 "TensorBSpline requires at least 2 basis coefficients to enforce sum-to-zero identifiability"
5521 );
5522 }
5523 let dense_design_ref = dense_design.as_ref().ok_or_else(|| {
5524 BasisError::InvalidInput(
5525 "tensor sum-to-zero identifiability requires a realized basis".to_string(),
5526 )
5527 })?;
5528 let (_, z) = apply_sum_to_zero_constraint(dense_design_ref.view(), None)?;
5529 let gauge = gam_problem::Gauge::sum_to_zero(z);
5530 Some(gauge.block_transform(0))
5531 }
5532 TensorBSplineIdentifiability::MarginalSumToZero => {
5533 if marginal_designs.len() < 2 {
5544 crate::bail_invalid_basis!(
5545 "tensor interaction (ti) identifiability requires at least 2 margins"
5546 );
5547 }
5548 let mut z = Array2::<f64>::eye(1);
5549 for (dim, marginal) in marginal_designs.iter().enumerate() {
5550 if marginal.ncols() < 2 {
5551 crate::bail_invalid_basis!(
5552 "tensor interaction (ti) margin {dim} has fewer than 2 basis functions; \
5553 cannot remove its marginal main effect"
5554 );
5555 }
5556 let (_, z_dim) = apply_sum_to_zero_constraint(marginal.view(), None)?;
5557 let gauge_dim = gam_problem::Gauge::sum_to_zero(z_dim);
5558 let z_dim = gauge_dim.block_transform(0);
5559 z = kronecker_product(&z, &z_dim);
5560 }
5561 Some(z)
5562 }
5563 TensorBSplineIdentifiability::FrozenTransform { transform } => {
5564 if transform.nrows() != total_cols {
5565 crate::bail_dim_basis!(
5566 "frozen tensor identifiability transform mismatch: design has {} columns but transform has {} rows",
5567 total_cols,
5568 transform.nrows()
5569 );
5570 }
5571 Some(transform.clone())
5572 }
5573 };
5574
5575 if let Some(z) = z_opt.as_ref() {
5576 let gauge = gam_problem::Gauge::from_block_transforms(&[z.clone()]);
5577 let dense = dense_design.as_mut().ok_or_else(|| {
5578 BasisError::InvalidInput(
5579 "tensor identifiability transform requires a realized basis".to_string(),
5580 )
5581 })?;
5582 let restricted_design = gauge.restrict_design(dense);
5583 *dense = restricted_design;
5584 candidates = candidates
5585 .into_iter()
5586 .map(|candidate| -> Result<PenaltyCandidate, BasisError> {
5587 let matrix = gauge.restrict_penalty(&candidate.matrix);
5588 let (matrix, c_new) = normalize_penalty_in_constrained_space(&matrix);
5596 Ok(PenaltyCandidate {
5597 nullspace_dim_hint: candidate.nullspace_dim_hint,
5598 matrix,
5599 source: candidate.source,
5600 normalization_scale: candidate.normalization_scale * c_new,
5601 kronecker_factors: None,
5607 op: candidate.op.clone(),
5608 })
5609 })
5610 .collect::<Result<Vec<_>, _>>()?;
5611 }
5612
5613 let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
5614 filter_active_penalty_candidates_with_ops(candidates)?;
5615 let identifiability_is_none =
5616 matches!(spec.identifiability, TensorBSplineIdentifiability::None);
5617 let all_marginals_sparse = marginal_sparse.iter().all(Option::is_some);
5625 let design = if let Some(dense_design) = dense_design {
5626 DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(dense_design))
5627 } else if identifiability_is_none && all_marginals_sparse {
5628 let sparse_marginals: Vec<&SparseColMat<usize, f64>> = marginal_sparse
5634 .iter()
5635 .map(|m| m.as_ref().expect("all_marginals_sparse just verified"))
5636 .collect();
5637 let sparse_design = tensor_product_design_from_sparse_marginals(&sparse_marginals)?;
5638 DesignMatrix::Sparse(gam_linalg::matrix::SparseDesignMatrix::new(sparse_design))
5639 } else {
5640 let marginals: Vec<Arc<Array2<f64>>> = marginal_designs
5641 .iter()
5642 .map(|m| Arc::new(m.clone()))
5643 .collect();
5644 let op = TensorProductDesignOperator::new(marginals).map_err(|e| {
5645 BasisError::InvalidInput(format!("TensorProductDesignOperator build failed: {e}"))
5646 })?;
5647 DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(Arc::new(op)))
5648 };
5649
5650 Ok(BasisBuildResult {
5651 design,
5652 penalties,
5653 nullspace_dims,
5654 penaltyinfo,
5655 ops,
5656 null_eigenvectors,
5657 joint_null_rotation: None,
5658 metadata: BasisMetadata::TensorBSpline {
5659 feature_cols: feature_cols.to_vec(),
5660 knots: marginal_knots,
5661 degrees: marginal_degrees,
5662 periods: marginal_effective_periods,
5669 is_cr: marginal_is_cr_flags,
5670 identifiability_transform: z_opt,
5671 },
5672 kronecker_factored: if matches!(spec.identifiability, TensorBSplineIdentifiability::None)
5673 && matches!(
5674 spec.penalty_decomposition,
5675 TensorBSplinePenaltyDecomposition::MarginalKroneckerSum
5676 ) {
5677 Some(KroneckerFactoredBasis::new(
5678 marginal_designs,
5679 kronecker_marginal_penalties,
5680 marginalnum_basis.clone(),
5681 spec.double_penalty,
5682 ))
5683 } else {
5684 None
5685 },
5686 })
5687}
5688
5689pub fn tensor_product_design_from_marginals(
5690 marginal_designs: &[Array2<f64>],
5691) -> Result<Array2<f64>, BasisError> {
5692 if marginal_designs.is_empty() {
5693 crate::bail_invalid_basis!("TensorBSpline requires at least one marginal basis");
5694 }
5695 let n = marginal_designs[0].nrows();
5696 for (i, b) in marginal_designs.iter().enumerate().skip(1) {
5697 if b.nrows() != n {
5698 crate::bail_dim_basis!(
5699 "tensor marginal row mismatch at dim {i}: expected {n}, got {}",
5700 b.nrows()
5701 );
5702 }
5703 }
5704 let total_cols = marginal_designs.iter().try_fold(1usize, |acc, b| {
5705 acc.checked_mul(b.ncols())
5706 .ok_or_else(|| BasisError::DimensionMismatch("tensor basis too large".to_string()))
5707 })?;
5708 use ndarray::parallel::prelude::*;
5714 use rayon::iter::{IntoParallelIterator, ParallelIterator};
5715 let mut design = Array2::<f64>::zeros((n, total_cols));
5716 design
5717 .axis_chunks_iter_mut(ndarray::Axis(0), 1024)
5718 .into_par_iter()
5719 .enumerate()
5720 .for_each(|(chunk_idx, mut block)| {
5721 let row_offset = chunk_idx * 1024;
5722 let mut cur = Vec::<f64>::with_capacity(total_cols);
5724 let mut next = Vec::<f64>::with_capacity(total_cols);
5725 for (local_i, mut out_row) in block.outer_iter_mut().enumerate() {
5726 let i = row_offset + local_i;
5727 cur.clear();
5728 cur.push(1.0);
5729 for b in marginal_designs {
5730 let q = b.ncols();
5731 next.clear();
5732 next.resize(cur.len() * q, 0.0);
5733 let b_row = b.row(i);
5737 let b_slice = b_row
5738 .as_slice()
5739 .expect("Array2 row from outer_iter is contiguous");
5740 for (a_idx, &aval) in cur.iter().enumerate() {
5741 let off = a_idx * q;
5742 let dst = &mut next[off..off + q];
5743 for col in 0..q {
5744 dst[col] = aval * b_slice[col];
5745 }
5746 }
5747 std::mem::swap(&mut cur, &mut next);
5748 }
5749 let out_slice = out_row
5754 .as_slice_mut()
5755 .expect("design row is contiguous in C-major Array2");
5756 out_slice.copy_from_slice(&cur);
5757 }
5758 });
5759 Ok(design)
5760}
5761
5762pub fn build_random_effect_block(
5763 data: ArrayView2<'_, f64>,
5764 spec: &RandomEffectTermSpec,
5765) -> Result<RandomEffectBlock, BasisError> {
5766 let n = data.nrows();
5767 let p = data.ncols();
5768 if spec.feature_col >= p {
5769 crate::bail_dim_basis!(
5770 "random-effect term '{}' feature column {} out of bounds for {} columns",
5771 spec.name,
5772 spec.feature_col,
5773 p
5774 );
5775 }
5776
5777 let col = data.column(spec.feature_col);
5778 if col.iter().any(|v| !v.is_finite()) {
5779 crate::bail_invalid_basis!(
5780 "random-effect term '{}' contains non-finite group values",
5781 spec.name
5782 );
5783 }
5784
5785 let kept_levels: Vec<u64> = if let Some(levels) = spec.frozen_levels.as_ref() {
5786 if levels.is_empty() {
5787 crate::bail_invalid_basis!(
5788 "random-effect term '{}' has empty frozen_levels",
5789 spec.name
5790 );
5791 }
5792 levels.clone()
5793 } else {
5794 let mut seen = BTreeSet::<u64>::new();
5795 let mut levels = Vec::<u64>::new();
5796 for &v in col {
5797 let bits = v.to_bits();
5798 if seen.insert(bits) {
5799 levels.push(bits);
5800 }
5801 }
5802 if levels.is_empty() {
5803 crate::bail_invalid_basis!("random-effect term '{}' has no observed levels", spec.name);
5804 }
5805 let start_idx = if spec.drop_first_level && levels.len() > 1 {
5806 1usize
5807 } else {
5808 0usize
5809 };
5810 levels[start_idx..].to_vec()
5811 };
5812
5813 if kept_levels.is_empty() {
5814 crate::bail_invalid_basis!(
5815 "random-effect term '{}' drops all levels; keep at least one level",
5816 spec.name
5817 );
5818 }
5819
5820 let q = kept_levels.len();
5821 let mut level_to_col = BTreeMap::<u64, usize>::new();
5822 for (idx, &bits) in kept_levels.iter().enumerate() {
5823 if level_to_col.insert(bits, idx).is_some() {
5824 crate::bail_invalid_basis!(
5825 "random-effect term '{}' has duplicate frozen level bits {bits}",
5826 spec.name
5827 );
5828 }
5829 }
5830 let mut group_ids = Vec::with_capacity(n);
5831 for &v in col {
5832 let bits = v.to_bits();
5833 group_ids.push(level_to_col.get(&bits).copied());
5834 }
5835
5836 Ok(RandomEffectBlock {
5837 name: spec.name.clone(),
5838 group_ids,
5839 num_groups: q,
5840 kept_levels,
5841 })
5842}
5843
5844impl SmoothDesign {
5845 pub fn map_term_coefficients(
5848 unconstrained: &Array1<f64>,
5849 shape: ShapeConstraint,
5850 ) -> Result<Array1<f64>, BasisError> {
5851 if unconstrained.is_empty() {
5852 crate::bail_invalid_basis!("unconstrained coefficient vector cannot be empty");
5853 }
5854 let mapped = match shape {
5855 ShapeConstraint::None => unconstrained.clone(),
5856 ShapeConstraint::MonotoneIncreasing => cumulative_exp(unconstrained, 1.0),
5857 ShapeConstraint::MonotoneDecreasing => cumulative_exp(unconstrained, -1.0),
5858 ShapeConstraint::Convex => second_cumulative_exp(unconstrained, 1.0),
5859 ShapeConstraint::Concave => second_cumulative_exp(unconstrained, -1.0),
5860 };
5861 Ok(mapped)
5862 }
5863}
5864
5865pub struct LocalSmoothTermBuild {
5866 pub dim: usize,
5867 pub design: DesignMatrix,
5868 pub penalties: Vec<Array2<f64>>,
5869 pub ops: Vec<Option<std::sync::Arc<dyn crate::analytic_penalties::PenaltyOp>>>,
5870 pub nullspaces: Vec<usize>,
5871 pub null_eigenvectors: Vec<Option<Array2<f64>>>,
5879 pub joint_null_rotation: Option<crate::basis::JointNullRotation>,
5886 pub penaltyinfo: Vec<PenaltyInfo>,
5887 pub pre_dropped_penaltyinfo: Vec<PenaltyInfo>,
5888 pub metadata: BasisMetadata,
5889 pub linear_constraints: Option<LinearInequalityConstraints>,
5890 pub box_reparam: bool,
5891 pub kronecker_factored: Option<KroneckerFactoredBasis>,
5892}
5893
5894#[derive(Clone)]
5895pub struct PcaScoresMemmapDesignOperator {
5896 mmap: Arc<memmap2::Mmap>,
5897 data_offset: usize,
5898 nrows: usize,
5899 ncols: usize,
5900 chunk_size: usize,
5901}
5902
5903impl PcaScoresMemmapDesignOperator {
5904 fn open(path: PathBuf, chunk_size: usize) -> Result<Self, BasisError> {
5905 let file = File::open(&path).map_err(|err| {
5906 BasisError::InvalidInput(format!(
5907 "failed to open lazy Pca .npy scores '{}': {err}",
5908 path.display()
5909 ))
5910 })?;
5911 let mmap = unsafe {
5917 memmap2::Mmap::map(&file).map_err(|err| {
5918 BasisError::InvalidInput(format!(
5919 "failed to memmap lazy Pca .npy scores '{}': {err}",
5920 path.display()
5921 ))
5922 })?
5923 };
5924 let (data_offset, nrows, ncols) = parse_f64_2d_npy_header(&mmap, &path)?;
5925 let expected = data_offset
5926 .checked_add(nrows.saturating_mul(ncols).saturating_mul(8))
5927 .ok_or_else(|| {
5928 BasisError::InvalidInput(format!(
5929 "lazy Pca .npy scores '{}' shape is too large",
5930 path.display()
5931 ))
5932 })?;
5933 if mmap.len() < expected {
5934 crate::bail_invalid_basis!(
5935 "lazy Pca .npy scores '{}' is truncated: header expects {} bytes, file has {}",
5936 path.display(),
5937 expected,
5938 mmap.len()
5939 );
5940 }
5941 Ok(Self {
5942 mmap: Arc::new(mmap),
5943 data_offset,
5944 nrows,
5945 ncols,
5946 chunk_size: chunk_size.max(1),
5947 })
5948 }
5949
5950 fn value(&self, row: usize, col: usize) -> f64 {
5951 let offset = self.data_offset + (row * self.ncols + col) * 8;
5952 let mut bytes = [0_u8; 8];
5953 bytes.copy_from_slice(&self.mmap[offset..offset + 8]);
5954 f64::from_le_bytes(bytes)
5955 }
5956
5957 fn chunk_rows(&self) -> usize {
5958 self.chunk_size.min(self.nrows.max(1))
5959 }
5960}
5961
5962impl LinearOperator for PcaScoresMemmapDesignOperator {
5963 fn nrows(&self) -> usize {
5964 self.nrows
5965 }
5966
5967 fn ncols(&self) -> usize {
5968 self.ncols
5969 }
5970
5971 fn apply(&self, vector: &Array1<f64>) -> Array1<f64> {
5972 assert_eq!(
5973 vector.len(),
5974 self.ncols,
5975 "lazy Pca apply vector length mismatch"
5976 );
5977 let mut out = Array1::<f64>::zeros(self.nrows);
5978 for start in (0..self.nrows).step_by(self.chunk_rows()) {
5979 let end = (start + self.chunk_rows()).min(self.nrows);
5980 for row in start..end {
5981 let mut acc = 0.0;
5982 for col in 0..self.ncols {
5983 acc += self.value(row, col) * vector[col];
5984 }
5985 out[row] = acc;
5986 }
5987 }
5988 out
5989 }
5990
5991 fn apply_transpose(&self, vector: &Array1<f64>) -> Array1<f64> {
5992 assert_eq!(
5993 vector.len(),
5994 self.nrows,
5995 "lazy Pca apply_transpose vector length mismatch"
5996 );
5997 let mut out = Array1::<f64>::zeros(self.ncols);
5998 for start in (0..self.nrows).step_by(self.chunk_rows()) {
5999 let end = (start + self.chunk_rows()).min(self.nrows);
6000 for row in start..end {
6001 let scale = vector[row];
6002 if scale == 0.0 {
6003 continue;
6004 }
6005 for col in 0..self.ncols {
6006 out[col] += scale * self.value(row, col);
6007 }
6008 }
6009 }
6010 out
6011 }
6012
6013 fn diag_xtw_x(&self, weights: &Array1<f64>) -> Result<Array2<f64>, String> {
6014 if weights.len() != self.nrows {
6015 return Err(format!(
6016 "lazy Pca diag_xtw_x weight length mismatch: weights={}, nrows={}",
6017 weights.len(),
6018 self.nrows
6019 ));
6020 }
6021 let mut gram = Array2::<f64>::zeros((self.ncols, self.ncols));
6022 for start in (0..self.nrows).step_by(self.chunk_rows()) {
6023 let end = (start + self.chunk_rows()).min(self.nrows);
6024 for row in start..end {
6025 let w = weights[row];
6026 if w == 0.0 {
6027 continue;
6028 }
6029 for a in 0..self.ncols {
6030 let xa = self.value(row, a);
6031 if xa == 0.0 {
6032 continue;
6033 }
6034 for b in a..self.ncols {
6035 gram[[a, b]] += w * xa * self.value(row, b);
6036 }
6037 }
6038 }
6039 }
6040 for a in 0..self.ncols {
6041 for b in 0..a {
6042 gram[[a, b]] = gram[[b, a]];
6043 }
6044 }
6045 Ok(gram)
6046 }
6047
6048 fn apply_weighted_normal(
6049 &self,
6050 weights: &Array1<f64>,
6051 vector: &Array1<f64>,
6052 penalty: Option<&Array2<f64>>,
6053 ridge: f64,
6054 ) -> Array1<f64> {
6055 assert_eq!(
6056 weights.len(),
6057 self.nrows,
6058 "lazy Pca weighted-normal weight mismatch"
6059 );
6060 assert_eq!(
6061 vector.len(),
6062 self.ncols,
6063 "lazy Pca weighted-normal vector mismatch"
6064 );
6065 let mut out = Array1::<f64>::zeros(self.ncols);
6066 for start in (0..self.nrows).step_by(self.chunk_rows()) {
6067 let end = (start + self.chunk_rows()).min(self.nrows);
6068 for row in start..end {
6069 let w = weights[row].max(0.0);
6070 if w == 0.0 {
6071 continue;
6072 }
6073 let mut row_dot = 0.0;
6074 for col in 0..self.ncols {
6075 row_dot += self.value(row, col) * vector[col];
6076 }
6077 if row_dot == 0.0 {
6078 continue;
6079 }
6080 let scaled = w * row_dot;
6081 for col in 0..self.ncols {
6082 out[col] += scaled * self.value(row, col);
6083 }
6084 }
6085 }
6086 if let Some(pen) = penalty {
6087 out += &pen.dot(vector);
6088 }
6089 if ridge > 0.0 {
6090 out += &vector.mapv(|x| ridge * x);
6091 }
6092 out
6093 }
6094}
6095
6096impl DenseDesignOperator for PcaScoresMemmapDesignOperator {
6097 fn compute_xtwy(&self, weights: &Array1<f64>, y: &Array1<f64>) -> Result<Array1<f64>, String> {
6098 if weights.len() != self.nrows || y.len() != self.nrows {
6099 return Err(format!(
6100 "lazy Pca compute_xtwy dimension mismatch: weights={}, y={}, nrows={}",
6101 weights.len(),
6102 y.len(),
6103 self.nrows
6104 ));
6105 }
6106 let mut out = Array1::<f64>::zeros(self.ncols);
6107 for start in (0..self.nrows).step_by(self.chunk_rows()) {
6108 let end = (start + self.chunk_rows()).min(self.nrows);
6109 for row in start..end {
6110 let scale = weights[row] * y[row];
6111 if scale == 0.0 {
6112 continue;
6113 }
6114 for col in 0..self.ncols {
6115 out[col] += scale * self.value(row, col);
6116 }
6117 }
6118 }
6119 Ok(out)
6120 }
6121
6122 fn row_chunk_into(
6123 &self,
6124 rows: Range<usize>,
6125 mut out: ArrayViewMut2<'_, f64>,
6126 ) -> Result<(), MatrixMaterializationError> {
6127 if rows.end > self.nrows || rows.start > rows.end {
6128 return Err(MatrixMaterializationError::MissingRowChunk {
6129 context: "lazy Pca row range out of bounds",
6130 });
6131 }
6132 if out.nrows() != rows.end - rows.start || out.ncols() != self.ncols {
6133 return Err(MatrixMaterializationError::MissingRowChunk {
6134 context: "lazy Pca row_chunk_into shape mismatch",
6135 });
6136 }
6137 for (local, row) in (rows.start..rows.end).enumerate() {
6138 for col in 0..self.ncols {
6139 out[[local, col]] = self.value(row, col);
6140 }
6141 }
6142 Ok(())
6143 }
6144
6145 fn to_dense(&self) -> Array2<f64> {
6146 let mut out = Array2::<f64>::zeros((self.nrows, self.ncols));
6147 self.row_chunk_into(0..self.nrows, out.view_mut())
6148 .expect("lazy Pca full materialization failed");
6149 out
6150 }
6151}
6152
6153pub fn parse_f64_2d_npy_header(
6154 bytes: &[u8],
6155 path: &PathBuf,
6156) -> Result<(usize, usize, usize), BasisError> {
6157 if bytes.len() < 10 || &bytes[0..6] != b"\x93NUMPY" {
6158 crate::bail_invalid_basis!("lazy Pca scores '{}' is not a .npy file", path.display());
6159 }
6160 let major = bytes[6];
6161 let header_len = match major {
6162 1 => u16::from_le_bytes([bytes[8], bytes[9]]) as usize,
6163 2 | 3 => {
6164 if bytes.len() < 12 {
6165 crate::bail_invalid_basis!(
6166 "lazy Pca scores '{}' has a truncated .npy header",
6167 path.display()
6168 );
6169 }
6170 u32::from_le_bytes([bytes[8], bytes[9], bytes[10], bytes[11]]) as usize
6171 }
6172 other => {
6173 crate::bail_invalid_basis!(
6174 "lazy Pca scores '{}' uses unsupported .npy version {}",
6175 path.display(),
6176 other
6177 );
6178 }
6179 };
6180 let header_start = if major == 1 { 10 } else { 12 };
6181 let data_offset = header_start + header_len;
6182 if bytes.len() < data_offset {
6183 crate::bail_invalid_basis!(
6184 "lazy Pca scores '{}' has a truncated .npy header",
6185 path.display()
6186 );
6187 }
6188 let header = std::str::from_utf8(&bytes[header_start..data_offset]).map_err(|err| {
6189 BasisError::InvalidInput(format!(
6190 "lazy Pca scores '{}' has a non-UTF8 .npy header: {err}",
6191 path.display()
6192 ))
6193 })?;
6194 if !(header.contains("'descr': '<f8'")
6195 || header.contains("\"descr\": \"<f8\"")
6196 || header.contains("'descr': '|f8'")
6197 || header.contains("\"descr\": \"|f8\""))
6198 {
6199 crate::bail_invalid_basis!(
6200 "lazy Pca scores '{}' must be float64 little-endian .npy",
6201 path.display()
6202 );
6203 }
6204 if header.contains("True") {
6205 crate::bail_invalid_basis!(
6206 "lazy Pca scores '{}' must be C-contiguous, not Fortran-ordered",
6207 path.display()
6208 );
6209 }
6210 let shape_pos = header.find("shape").ok_or_else(|| {
6211 BasisError::InvalidInput(format!(
6212 "lazy Pca scores '{}' .npy header is missing shape",
6213 path.display()
6214 ))
6215 })?;
6216 let open = header[shape_pos..].find('(').ok_or_else(|| {
6217 BasisError::InvalidInput(format!(
6218 "lazy Pca scores '{}' .npy header has malformed shape",
6219 path.display()
6220 ))
6221 })? + shape_pos;
6222 let close = header[open..].find(')').ok_or_else(|| {
6223 BasisError::InvalidInput(format!(
6224 "lazy Pca scores '{}' .npy header has malformed shape",
6225 path.display()
6226 ))
6227 })? + open;
6228 let dims = header[open + 1..close]
6229 .split(',')
6230 .map(str::trim)
6231 .filter(|part| !part.is_empty())
6232 .map(|part| part.parse::<usize>())
6233 .collect::<Result<Vec<_>, _>>()
6234 .map_err(|err| {
6235 BasisError::InvalidInput(format!(
6236 "lazy Pca scores '{}' .npy shape is not integral: {err}",
6237 path.display()
6238 ))
6239 })?;
6240 if dims.len() != 2 {
6241 crate::bail_invalid_basis!(
6242 "lazy Pca scores '{}' must have shape (N, K), got {:?}",
6243 path.display(),
6244 dims
6245 );
6246 }
6247 Ok((data_offset, dims[0], dims[1]))
6248}
6249
6250pub fn pca_center_mean(x: ArrayView2<'_, f64>) -> Result<Array1<f64>, BasisError> {
6251 if x.nrows() == 0 {
6252 crate::bail_invalid_basis!("Pca basis requires at least one row to compute center mean");
6253 }
6254 let mut mean = Array1::<f64>::zeros(x.ncols());
6255 for row in x.rows() {
6256 mean += &row;
6257 }
6258 mean.mapv_inplace(|v| v / x.nrows() as f64);
6259 Ok(mean)
6260}
6261
6262pub fn build_pca_smooth_basis(
6263 data: ArrayView2<'_, f64>,
6264 feature_cols: &[usize],
6265 basis_matrix: &Array2<f64>,
6266 centered: bool,
6267 smooth_penalty: f64,
6268 center_mean: Option<&Array1<f64>>,
6269 pca_basis_path: Option<&PathBuf>,
6270 chunk_size: usize,
6271) -> Result<BasisBuildResult, BasisError> {
6272 if let Some(path) = pca_basis_path {
6273 let op = PcaScoresMemmapDesignOperator::open(path.clone(), chunk_size)?;
6274 if op.nrows != data.nrows() {
6275 crate::bail_dim_basis!(
6276 "lazy Pca scores row mismatch: .npy has {}, data has {}",
6277 op.nrows,
6278 data.nrows()
6279 );
6280 }
6281 let k = op.ncols;
6282 let mut penalty = Array2::<f64>::eye(k);
6283 penalty.mapv_inplace(|v| v * smooth_penalty);
6284 let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
6285 filter_active_penalty_candidates_with_ops(vec![PenaltyCandidate {
6286 matrix: penalty,
6287 nullspace_dim_hint: 0,
6288 source: PenaltySource::Other("PcaRidge".to_string()),
6289 normalization_scale: 1.0,
6290 kronecker_factors: None,
6291 op: None,
6292 }])?;
6293 return Ok(BasisBuildResult {
6294 design: DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(Arc::new(op))),
6295 penalties,
6296 nullspace_dims,
6297 penaltyinfo,
6298 ops,
6299 null_eigenvectors,
6300 joint_null_rotation: None,
6301 metadata: BasisMetadata::Pca {
6302 feature_cols: feature_cols.to_vec(),
6303 basis_matrix: basis_matrix.clone(),
6304 centered,
6305 smooth_penalty,
6306 center_mean: center_mean.cloned(),
6307 pca_basis_path: Some(path.clone()),
6308 chunk_size: chunk_size.max(1),
6309 },
6310 kronecker_factored: None,
6311 });
6312 }
6313 if basis_matrix.nrows() != feature_cols.len() {
6314 crate::bail_dim_basis!(
6315 "Pca basis row mismatch: basis rows={}, feature columns={}",
6316 basis_matrix.nrows(),
6317 feature_cols.len()
6318 );
6319 }
6320 let mut x = select_columns(data, feature_cols)?;
6321 let mean = if centered {
6322 match center_mean {
6323 Some(mean) => mean.clone(),
6324 None => pca_center_mean(x.view())?,
6325 }
6326 } else {
6327 Array1::<f64>::zeros(feature_cols.len())
6328 };
6329 if centered {
6330 for mut row in x.rows_mut() {
6331 row -= &mean;
6332 }
6333 }
6334 let design = fast_ab(&x, basis_matrix);
6335 let k = basis_matrix.ncols();
6336 let mut penalty = Array2::<f64>::eye(k);
6337 penalty.mapv_inplace(|v| v * smooth_penalty);
6338 let (penalties, nullspace_dims, penaltyinfo, null_eigenvectors, ops) =
6339 filter_active_penalty_candidates_with_ops(vec![PenaltyCandidate {
6340 matrix: penalty,
6341 nullspace_dim_hint: 0,
6342 source: PenaltySource::Other("PcaRidge".to_string()),
6343 normalization_scale: 1.0,
6344 kronecker_factors: None,
6345 op: None,
6346 }])?;
6347 Ok(BasisBuildResult {
6348 design: DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(design)),
6349 penalties,
6350 nullspace_dims,
6351 penaltyinfo,
6352 ops,
6353 null_eigenvectors,
6354 joint_null_rotation: None,
6355 metadata: BasisMetadata::Pca {
6356 feature_cols: feature_cols.to_vec(),
6357 basis_matrix: basis_matrix.clone(),
6358 centered,
6359 smooth_penalty,
6360 center_mean: centered.then_some(mean),
6361 pca_basis_path: None,
6362 chunk_size: chunk_size.max(1),
6363 },
6364 kronecker_factored: None,
6365 })
6366}
6367
6368pub fn defer_inner_model_centering_to_factor_level_wrapper(basis: &mut SmoothBasisSpec) {
6384 if let SmoothBasisSpec::BSpline1D { spec, .. } = basis
6385 && matches!(
6386 spec.identifiability,
6387 BSplineIdentifiability::WeightedSumToZero { .. }
6388 )
6389 {
6390 spec.identifiability = BSplineIdentifiability::None;
6391 }
6392}
6393
6394pub fn apply_by_variable_to_local_build(
6395 mut built: LocalSmoothTermBuild,
6396 data: ArrayView2<'_, f64>,
6397 by_col: usize,
6398 by: &ByVariableSpec,
6399 term_name: &str,
6400) -> Result<LocalSmoothTermBuild, BasisError> {
6401 if by_col >= data.ncols() {
6402 crate::bail_dim_basis!(
6403 "by-variable smooth term '{term_name}' references column {by_col}, but data has {} columns",
6404 data.ncols()
6405 );
6406 }
6407 let weights = match by {
6408 ByVariableSpec::Numeric => data.column(by_col).to_owned(),
6409 ByVariableSpec::Level { value_bits, .. } => data.column(by_col).mapv(|value| {
6410 if value.to_bits() == *value_bits {
6411 1.0
6412 } else {
6413 0.0
6414 }
6415 }),
6416 };
6417 if weights.iter().any(|value| !value.is_finite()) {
6418 crate::bail_invalid_basis!(
6419 "by-variable smooth term '{term_name}' has non-finite by-column values"
6420 );
6421 }
6422
6423 let mut dense = built
6424 .design
6425 .try_to_dense_by_chunks("by-variable smooth row gating")
6426 .map_err(BasisError::InvalidInput)?;
6427 for (mut row, &weight) in dense.rows_mut().into_iter().zip(weights.iter()) {
6428 row.mapv_inplace(|value| value * weight);
6429 }
6430 built.design = DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(dense));
6431 built.kronecker_factored = None;
6432 Ok(built)
6433}
6434
6435pub fn build_by_smooth_local(
6446 data: ArrayView2<'_, f64>,
6447 term: &SmoothTermSpec,
6448 smooth: &SmoothBasisSpec,
6449 by_kind: &ByVarKind,
6450 workspace: &mut crate::basis::BasisWorkspace,
6451) -> Result<LocalSmoothTermBuild, BasisError> {
6452 let inner_term = SmoothTermSpec {
6453 name: term.name.clone(),
6454 basis: (*smooth).clone(),
6455 shape: term.shape,
6456 joint_null_rotation: None,
6457 };
6458 let inner = build_single_local_smooth_term(data, &inner_term, workspace)?;
6459
6460 match by_kind {
6461 ByVarKind::Numeric { feature_col } => {
6462 let inner_meta = inner.metadata.clone();
6463 let mut built = apply_by_variable_to_local_build(
6464 inner,
6465 data,
6466 *feature_col,
6467 &ByVariableSpec::Numeric,
6468 &term.name,
6469 )?;
6470 built.metadata = BasisMetadata::BySmooth {
6471 inner: Box::new(inner_meta),
6472 by_col: *feature_col,
6473 levels: None,
6474 ordered: false,
6475 };
6476 Ok(built)
6477 }
6478 ByVarKind::Factor {
6479 feature_col,
6480 frozen_levels,
6481 ordered,
6482 } => {
6483 let level_bits: Vec<u64> = if let Some(fl) = frozen_levels {
6486 fl.clone()
6487 } else {
6488 let col = data.column(*feature_col);
6489 let mut seen = BTreeSet::<u64>::new();
6490 for &v in col.iter() {
6491 if v.is_finite() {
6492 seen.insert(v.to_bits());
6493 }
6494 }
6495 seen.into_iter().collect()
6496 };
6497 let n_levels = level_bits.len();
6498 if n_levels == 0 {
6499 crate::bail_invalid_basis!(
6500 "by-factor smooth term '{}': factor column {} has no observed levels",
6501 term.name,
6502 feature_col
6503 );
6504 }
6505 let p = inner.dim;
6506 let q = n_levels * p;
6507 let n = data.nrows();
6508
6509 let inner_dense = inner
6510 .design
6511 .try_to_dense_by_chunks("by-factor smooth design gating")
6512 .map_err(BasisError::InvalidInput)?;
6513
6514 let mut combined = Array2::<f64>::zeros((n, q));
6516 for (lvl_idx, &bits) in level_bits.iter().enumerate() {
6517 let col_start = lvl_idx * p;
6518 for row in 0..n {
6519 if data[[row, *feature_col]].to_bits() == bits {
6520 combined
6521 .slice_mut(s![row, col_start..col_start + p])
6522 .assign(&inner_dense.row(row));
6523 }
6524 }
6525 }
6526
6527 let inner_meta = inner.metadata.clone();
6539 let n_penalties = inner.penalties.len();
6540 let n_blocks = n_penalties.saturating_mul(n_levels);
6541 let mut penalties = Vec::<Array2<f64>>::with_capacity(n_blocks);
6542 let mut penaltyinfo = Vec::<PenaltyInfo>::with_capacity(n_blocks);
6543 let mut nullspaces = Vec::<usize>::with_capacity(n_blocks);
6544 for (pen_pos, s_inner) in inner.penalties.iter().enumerate() {
6545 for lvl in 0..n_levels {
6546 let off = lvl * p;
6547 let mut s_big = Array2::<f64>::zeros((q, q));
6548 s_big
6549 .slice_mut(s![off..off + p, off..off + p])
6550 .assign(s_inner);
6551 let (s_big, scale) = normalize_penalty_in_constrained_space(&s_big);
6552 let mut info = inner.penaltyinfo[pen_pos].clone();
6553 info.original_index = pen_pos * n_levels + lvl;
6556 info.normalization_scale *= scale;
6557 info.kronecker_factors = None;
6560 penalties.push(s_big);
6561 penaltyinfo.push(info);
6562 nullspaces.push(inner.nullspaces[pen_pos]);
6563 }
6564 }
6565
6566 let null_eigenvectors = vec![None; penalties.len()];
6567 let ops = vec![None; penalties.len()];
6568
6569 Ok(LocalSmoothTermBuild {
6570 dim: q,
6571 design: DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(combined)),
6572 penalties,
6573 ops,
6574 nullspaces,
6575 null_eigenvectors,
6576 joint_null_rotation: None,
6577 penaltyinfo,
6578 pre_dropped_penaltyinfo: inner.pre_dropped_penaltyinfo,
6579 metadata: BasisMetadata::BySmooth {
6580 inner: Box::new(inner_meta),
6581 by_col: *feature_col,
6582 levels: Some(level_bits),
6583 ordered: *ordered,
6584 },
6585 linear_constraints: None,
6586 box_reparam: false,
6587 kronecker_factored: None,
6588 })
6589 }
6590 }
6591}
6592
6593pub fn ensure_by_variable_specs_match(
6594 kind: &BySmoothKind,
6595 by: &ByVariableSpec,
6596 term_name: &str,
6597) -> Result<(), BasisError> {
6598 match (kind, by) {
6599 (BySmoothKind::Numeric, ByVariableSpec::Numeric) => Ok(()),
6600 (BySmoothKind::Level { level_bits }, ByVariableSpec::Level { value_bits, .. })
6601 if level_bits == value_bits =>
6602 {
6603 Ok(())
6604 }
6605 _ => Err(BasisError::InvalidInput(format!(
6606 "by-variable smooth term '{term_name}' has inconsistent by-variable specifications"
6607 ))),
6608 }
6609}
6610
6611pub fn build_factor_smooth(
6639 data: ArrayView2<'_, f64>,
6640 spec: &FactorSmoothSpec,
6641 term_name: &str,
6642 workspace: &mut crate::basis::BasisWorkspace,
6643) -> Result<LocalSmoothTermBuild, BasisError> {
6644 if spec.continuous_cols.len() != 1 {
6645 crate::bail_invalid_basis!(
6646 "factor smooth term '{}' currently supports exactly one continuous covariate; found {}",
6647 term_name,
6648 spec.continuous_cols.len()
6649 );
6650 }
6651 let feature_col = spec.continuous_cols[0];
6652 let group_col = spec.group_col;
6653 if feature_col >= data.ncols() || group_col >= data.ncols() {
6654 crate::bail_dim_basis!(
6655 "factor smooth term '{}' references columns ({}, {}) out of bounds for {} columns",
6656 term_name,
6657 feature_col,
6658 group_col,
6659 data.ncols()
6660 );
6661 }
6662
6663 if matches!(spec.flavour, FactorSmoothFlavour::Sz) {
6666 let levels = resolve_factor_smooth_levels(data, group_col, spec, term_name)?;
6667 let inner = SmoothBasisSpec::BSpline1D {
6668 feature_col,
6669 spec: factor_smooth_marginal_for_replay(&spec.marginal),
6670 };
6671 let sz_term = SmoothTermSpec {
6672 name: term_name.to_string(),
6673 basis: SmoothBasisSpec::FactorSumToZero {
6674 inner: Box::new(inner),
6675 by_col: group_col,
6676 levels: levels.clone(),
6677 frozen_global_orthogonality: None,
6678 },
6679 shape: ShapeConstraint::None,
6680 joint_null_rotation: None,
6681 };
6682 let mut built = build_single_local_smooth_term(data, &sz_term, workspace)?;
6683 let (knots, degree, periodic, marginal_is_cr) = match &built.metadata {
6704 BasisMetadata::BSpline1D {
6705 knots,
6706 periodic,
6707 degree,
6708 ..
6709 } => (
6710 knots.clone(),
6711 degree.unwrap_or(spec.marginal.degree),
6712 *periodic,
6713 false,
6714 ),
6715 BasisMetadata::CubicRegression1D { knots, .. } => {
6716 (knots.clone(), spec.marginal.degree, None, true)
6717 }
6718 other => {
6719 crate::bail_invalid_basis!(
6720 "sz factor smooth term '{}' produced an unexpected marginal metadata variant {:?}",
6721 term_name,
6722 other
6723 );
6724 }
6725 };
6726 built.metadata = BasisMetadata::FactorSmooth {
6727 continuous_cols: spec.continuous_cols.clone(),
6728 group_col,
6729 knots,
6730 degree,
6731 periodic,
6732 group_levels: levels,
6733 flavour: "sz".to_string(),
6734 marginal_is_cr,
6735 };
6736 return Ok(built);
6737 }
6738
6739 let levels = resolve_factor_smooth_levels(data, group_col, spec, term_name)?;
6740 let n_levels = levels.len();
6741 if n_levels < 2 {
6742 crate::bail_invalid_basis!(
6743 "factor smooth term '{}' requires at least two grouping levels; found {}",
6744 term_name,
6745 n_levels
6746 );
6747 }
6748
6749 let use_per_dim_null = matches!(
6757 &spec.flavour,
6758 FactorSmoothFlavour::Fs { m_null_penalty_orders }
6759 if m_null_penalty_orders.iter().copied().max().unwrap_or(0) >= 1
6760 );
6761
6762 let mut marginal_spec = factor_smooth_marginal_for_replay(&spec.marginal);
6768 if use_per_dim_null {
6769 marginal_spec.double_penalty = false;
6770 }
6771 let inner_term = SmoothTermSpec {
6772 name: format!("{term_name}::marginal"),
6773 basis: SmoothBasisSpec::BSpline1D {
6774 feature_col,
6775 spec: marginal_spec,
6776 },
6777 shape: ShapeConstraint::None,
6778 joint_null_rotation: None,
6779 };
6780 let inner = build_single_local_smooth_term(data, &inner_term, workspace)?;
6781 let mut base = inner
6782 .design
6783 .try_to_dense_by_chunks("factor smooth marginal")
6784 .map_err(BasisError::InvalidInput)?;
6785 if matches!(spec.flavour, FactorSmoothFlavour::Re) {
6786 let center = match &inner.metadata {
6796 BasisMetadata::BSpline1D { knots, .. } if !knots.is_empty() => {
6797 0.5 * (knots[0] + knots[knots.len() - 1])
6798 }
6799 _ => 0.0,
6800 };
6801 let mut linear = Array2::<f64>::ones((data.nrows(), 2));
6802 linear
6803 .column_mut(1)
6804 .assign(&data.column(feature_col).mapv(|x| x - center));
6805 base = linear;
6806 }
6807 let n = base.nrows();
6808 let p = base.ncols();
6809 let q = p * n_levels;
6810
6811 let mut dense = Array2::<f64>::zeros((n, q));
6814 for i in 0..n {
6815 let bits = data[[i, group_col]].to_bits();
6816 let level_idx = levels.iter().position(|b| *b == bits).ok_or_else(|| {
6817 BasisError::InvalidInput(format!(
6818 "factor smooth term '{term_name}' saw an unseen grouping level at row {}",
6819 i + 1
6820 ))
6821 })?;
6822 let start = level_idx * p;
6823 dense
6824 .slice_mut(s![i, start..start + p])
6825 .assign(&base.row(i));
6826 }
6827
6828 let marginal_penalties: Vec<Array2<f64>> = if matches!(spec.flavour, FactorSmoothFlavour::Re) {
6834 (0..p)
6835 .map(|j| {
6836 let mut s = Array2::<f64>::zeros((p, p));
6837 s[[j, j]] = 1.0;
6838 s
6839 })
6840 .collect()
6841 } else {
6842 inner.penalties.clone()
6843 };
6844 let marginal_penaltyinfo: Vec<PenaltyInfo> = if matches!(spec.flavour, FactorSmoothFlavour::Re)
6845 {
6846 (0..p)
6847 .map(|j| PenaltyInfo {
6848 source: PenaltySource::Primary,
6849 original_index: j,
6850 active: true,
6851 effective_rank: 1,
6852 dropped_reason: None,
6853 nullspace_dim_hint: p.saturating_sub(1),
6854 normalization_scale: 1.0,
6855 kronecker_factors: None,
6856 })
6857 .collect()
6858 } else {
6859 inner.penaltyinfo.clone()
6860 };
6861 if marginal_penalties.len() != marginal_penaltyinfo.len() {
6862 crate::bail_invalid_basis!(
6863 "internal factor-smooth penalty metadata mismatch for term '{}': penalties={}, infos={}",
6864 term_name,
6865 marginal_penalties.len(),
6866 marginal_penaltyinfo.len()
6867 );
6868 }
6869
6870 let mut penalties = Vec::<Array2<f64>>::with_capacity(marginal_penalties.len());
6871 let mut penaltyinfo = Vec::<PenaltyInfo>::with_capacity(marginal_penalties.len());
6872 for (penalty_pos, s_inner) in marginal_penalties.iter().enumerate() {
6873 let mut s_big = Array2::<f64>::zeros((q, q));
6874 for level in 0..n_levels {
6875 let start = level * p;
6876 s_big
6877 .slice_mut(s![start..start + p, start..start + p])
6878 .assign(s_inner);
6879 }
6880 let (s_big, factor_smooth_scale) = normalize_penalty_in_constrained_space(&s_big);
6881 let mut info = marginal_penaltyinfo[penalty_pos].clone();
6882 info.original_index = penalty_pos;
6883 info.normalization_scale *= factor_smooth_scale;
6884 info.nullspace_dim_hint = info.nullspace_dim_hint.saturating_mul(n_levels);
6885 info.kronecker_factors = None;
6886 penalties.push(s_big);
6887 penaltyinfo.push(info);
6888 }
6889
6890 let mut nullspaces: Vec<usize> = if matches!(spec.flavour, FactorSmoothFlavour::Re) {
6891 vec![q.saturating_sub(n_levels); p]
6892 } else {
6893 inner
6894 .nullspaces
6895 .iter()
6896 .map(|ns| ns.saturating_mul(n_levels))
6897 .collect()
6898 };
6899
6900 if use_per_dim_null
6930 && let Some(Some(z)) = inner.null_eigenvectors.first()
6931 && z.nrows() == p
6932 {
6933 for k in 0..z.ncols() {
6934 let zk = z.column(k);
6939 let mut p_k = Array2::<f64>::zeros((p, p));
6940 for a in 0..p {
6941 for b in 0..p {
6942 p_k[[a, b]] = zk[a] * zk[b];
6943 }
6944 }
6945 let mut s_null = Array2::<f64>::zeros((q, q));
6946 for level in 0..n_levels {
6947 let start = level * p;
6948 s_null
6949 .slice_mut(s![start..start + p, start..start + p])
6950 .assign(&p_k);
6951 }
6952 let (s_null, null_scale) = normalize_penalty_in_constrained_space(&s_null);
6953 let null_block = crate::basis::analyze_penalty_block_with_op(&s_null, None)?;
6954 if null_block.rank > 0 {
6955 let original_index = penalties.len();
6956 penalties.push(null_block.sym_penalty);
6957 nullspaces.push(null_block.nullity);
6958 penaltyinfo.push(PenaltyInfo {
6959 source: PenaltySource::Primary,
6960 original_index,
6961 active: true,
6962 effective_rank: null_block.rank,
6963 dropped_reason: None,
6964 nullspace_dim_hint: null_block.nullity,
6965 normalization_scale: null_scale,
6966 kronecker_factors: None,
6967 });
6968 }
6969 }
6970 }
6971 let null_eigenvectors = crate::basis::recompute_null_eigenvectors(&penalties)?;
6972 let joint_null_rotation = crate::basis::compute_joint_null_rotation(&penalties)?;
6973
6974 let (knots, degree, periodic) = match &inner.metadata {
6977 BasisMetadata::BSpline1D {
6978 knots,
6979 periodic,
6980 degree,
6981 ..
6982 } => (
6983 knots.clone(),
6984 degree.unwrap_or(spec.marginal.degree),
6985 *periodic,
6986 ),
6987 other => {
6988 crate::bail_invalid_basis!(
6989 "factor smooth term '{}' produced an unexpected marginal metadata variant {:?}",
6990 term_name,
6991 other
6992 );
6993 }
6994 };
6995 let flavour_tag = match &spec.flavour {
6996 FactorSmoothFlavour::Fs { .. } => "fs",
6997 FactorSmoothFlavour::Sz => "sz",
6998 FactorSmoothFlavour::Re => "re",
6999 }
7000 .to_string();
7001 let metadata = BasisMetadata::FactorSmooth {
7002 continuous_cols: spec.continuous_cols.clone(),
7003 group_col,
7004 knots,
7005 degree,
7006 periodic,
7007 group_levels: levels,
7008 flavour: flavour_tag,
7009 marginal_is_cr: false,
7012 };
7013
7014 let ops = vec![None; penalties.len()];
7015 Ok(LocalSmoothTermBuild {
7016 dim: q,
7017 design: DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(dense)),
7018 penalties,
7019 ops,
7020 nullspaces,
7021 null_eigenvectors,
7022 joint_null_rotation,
7023 penaltyinfo,
7024 pre_dropped_penaltyinfo: Vec::new(),
7025 metadata,
7026 linear_constraints: None,
7027 box_reparam: false,
7028 kronecker_factored: None,
7029 })
7030}
7031
7032pub fn resolve_factor_smooth_levels(
7036 data: ArrayView2<'_, f64>,
7037 group_col: usize,
7038 spec: &FactorSmoothSpec,
7039 term_name: &str,
7040) -> Result<Vec<u64>, BasisError> {
7041 if let Some(frozen) = &spec.group_frozen_levels {
7042 if frozen.is_empty() {
7043 crate::bail_invalid_basis!(
7044 "factor smooth term '{}' has an empty frozen level list",
7045 term_name
7046 );
7047 }
7048 return Ok(frozen.clone());
7049 }
7050 let mut bits: Vec<u64> = data.column(group_col).iter().map(|v| v.to_bits()).collect();
7051 bits.sort_by(|a, b| {
7052 f64::from_bits(*a)
7053 .partial_cmp(&f64::from_bits(*b))
7054 .unwrap_or(std::cmp::Ordering::Equal)
7055 });
7056 bits.dedup();
7057 Ok(bits)
7058}
7059
7060pub fn factor_smooth_marginal_for_replay(marginal: &BSplineBasisSpec) -> BSplineBasisSpec {
7067 let mut m = marginal.clone();
7068 m.identifiability = BSplineIdentifiability::None;
7069 m
7070}
7071
7072pub fn build_single_local_smooth_term(
7073 data: ArrayView2<'_, f64>,
7074 term: &SmoothTermSpec,
7075 workspace: &mut crate::basis::BasisWorkspace,
7076) -> Result<LocalSmoothTermBuild, BasisError> {
7077 if term.shape != ShapeConstraint::None && !shape_supports_basis(term) {
7078 crate::bail_invalid_basis!(
7079 "ShapeConstraint::{:?} is unsupported for term '{}'",
7080 term.shape,
7081 term.name
7082 );
7083 }
7084 if let SmoothBasisSpec::ByVariable {
7085 inner,
7086 by_col,
7087 kind,
7088 by,
7089 } = &term.basis
7090 {
7091 ensure_by_variable_specs_match(kind, by, &term.name)?;
7092 let mut inner_basis = (**inner).clone();
7093 if matches!(by, ByVariableSpec::Level { .. }) {
7100 defer_inner_model_centering_to_factor_level_wrapper(&mut inner_basis);
7101 }
7102 let inner_term = SmoothTermSpec {
7103 name: term.name.clone(),
7104 basis: inner_basis,
7105 shape: term.shape,
7106 joint_null_rotation: None,
7107 };
7108 let built = build_single_local_smooth_term(data, &inner_term, workspace)?;
7109 return apply_by_variable_to_local_build(built, data, *by_col, by, &term.name);
7110 }
7111
7112 if let SmoothBasisSpec::BySmooth { smooth, by_kind } = &term.basis {
7115 return build_by_smooth_local(data, term, smooth, by_kind, workspace);
7116 }
7117
7118 let mut shape_axis_col: Option<usize> = None;
7119 let mut built: BasisBuildResult = match &term.basis {
7120 SmoothBasisSpec::FactorSumToZero {
7121 inner,
7122 by_col,
7123 levels,
7124 ..
7125 } => {
7126 if *by_col >= data.ncols() {
7127 crate::bail_dim_basis!(
7128 "term '{}' by column {} out of bounds for {} columns",
7129 term.name,
7130 by_col,
7131 data.ncols()
7132 );
7133 }
7134 if levels.len() < 2 {
7135 crate::bail_invalid_basis!(
7136 "sum-to-zero factor smooth term '{}' requires at least two levels",
7137 term.name
7138 );
7139 }
7140 if term.shape != ShapeConstraint::None {
7141 crate::bail_invalid_basis!(
7142 "ShapeConstraint::{:?} is unsupported for sum-to-zero factor smooth term '{}'",
7143 term.shape,
7144 term.name
7145 );
7146 }
7147 let inner_term = SmoothTermSpec {
7148 name: format!("{}::inner", term.name),
7149 basis: (**inner).clone(),
7150 shape: ShapeConstraint::None,
7151 joint_null_rotation: None,
7152 };
7153 let mut inner_built = build_single_local_smooth_term(data, &inner_term, workspace)?;
7154 let inner_null_eigenvectors = inner_built.null_eigenvectors.clone();
7158 let base = inner_built
7159 .design
7160 .try_to_dense_by_chunks("sum-to-zero factor smooth")
7161 .map_err(BasisError::InvalidInput)?;
7162 let n = base.nrows();
7163 let p = base.ncols();
7164 let l_minus_one = levels.len() - 1;
7165 let mut dense = Array2::<f64>::zeros((n, p * l_minus_one));
7166 for i in 0..n {
7167 let bits = data[[i, *by_col]].to_bits();
7168 let level_idx = levels.iter().position(|b| *b == bits).ok_or_else(|| {
7169 BasisError::InvalidInput(format!(
7170 "sum-to-zero factor smooth term '{}' saw an unseen level at row {}",
7171 term.name,
7172 i + 1
7173 ))
7174 })?;
7175 if level_idx < l_minus_one {
7176 let start = level_idx * p;
7177 dense
7178 .slice_mut(s![i, start..start + p])
7179 .assign(&base.row(i));
7180 } else {
7181 for level in 0..l_minus_one {
7182 let start = level * p;
7183 dense
7184 .slice_mut(s![i, start..start + p])
7185 .assign(&base.row(i).mapv(|v| -v));
7186 }
7187 }
7188 }
7189 let mut penalties = Vec::<Array2<f64>>::with_capacity(inner_built.penalties.len());
7190 let active_penalty_indices = inner_built
7191 .penaltyinfo
7192 .iter()
7193 .enumerate()
7194 .filter_map(|(idx, info)| info.active.then_some(idx))
7195 .collect::<Vec<_>>();
7196 if active_penalty_indices.len() != inner_built.penalties.len() {
7197 crate::bail_invalid_basis!(
7198 "internal sz penalty metadata mismatch: activeinfos={}, penalties={}",
7199 active_penalty_indices.len(),
7200 inner_built.penalties.len()
7201 );
7202 }
7203 let stz_per_group_penalty = |s_inner: &Array2<f64>, which_level: usize| -> Array2<f64> {
7238 let mut s_big = Array2::<f64>::zeros((p * l_minus_one, p * l_minus_one));
7239 if which_level < l_minus_one {
7240 let k = which_level;
7242 let mut block = s_big.slice_mut(s![k * p..(k + 1) * p, k * p..(k + 1) * p]);
7243 block.assign(s_inner);
7244 } else {
7245 for a in 0..l_minus_one {
7247 for b in 0..l_minus_one {
7248 let mut block =
7249 s_big.slice_mut(s![a * p..(a + 1) * p, b * p..(b + 1) * p]);
7250 block.assign(s_inner);
7251 }
7252 }
7253 }
7254 s_big
7255 };
7256 let mut nullspaces = Vec::<usize>::with_capacity(penalties.capacity());
7262 for (penalty_pos, s_inner) in inner_built.penalties.iter().enumerate() {
7263 let info_idx = active_penalty_indices[penalty_pos];
7264 let base_info = inner_built.penaltyinfo[info_idx].clone();
7265 let marginal_nullity = inner_built.nullspaces.get(penalty_pos).copied().unwrap_or(0);
7266 for which_level in 0..=l_minus_one {
7268 let raw = stz_per_group_penalty(s_inner, which_level);
7269 let (s_big, group_scale) = normalize_penalty_in_constrained_space(&raw);
7270 let block = crate::basis::analyze_penalty_block_with_op(&s_big, None)?;
7271 if block.rank == 0 {
7272 continue;
7273 }
7274 if which_level == 0 {
7275 inner_built.penaltyinfo[info_idx].normalization_scale *= group_scale;
7278 inner_built.penaltyinfo[info_idx].original_index = penalties.len();
7279 inner_built.penaltyinfo[info_idx].effective_rank = block.rank;
7280 inner_built.penaltyinfo[info_idx].nullspace_dim_hint = block.nullity;
7281 } else {
7282 let mut info = base_info.clone();
7283 info.original_index = penalties.len();
7284 info.normalization_scale = base_info.normalization_scale * group_scale;
7285 info.effective_rank = block.rank;
7286 info.nullspace_dim_hint = block.nullity;
7287 info.kronecker_factors = None;
7288 inner_built.penaltyinfo.push(info);
7289 }
7290 penalties.push(block.sym_penalty);
7291 nullspaces.push(marginal_nullity);
7297 }
7298 }
7299
7300 if let Some(Some(z)) = inner_null_eigenvectors.first()
7318 && z.nrows() == p
7319 {
7320 for k in 0..z.ncols() {
7321 let zk = z.column(k);
7322 let mut p_k = Array2::<f64>::zeros((p, p));
7323 for a in 0..p {
7324 for b in 0..p {
7325 p_k[[a, b]] = zk[a] * zk[b];
7326 }
7327 }
7328 let stz_pooled_null = {
7333 let mut s_big = Array2::<f64>::zeros((p * l_minus_one, p * l_minus_one));
7334 for a in 0..l_minus_one {
7335 for b in 0..l_minus_one {
7336 let factor = if a == b { 2.0 } else { 1.0 };
7337 let mut block =
7338 s_big.slice_mut(s![a * p..(a + 1) * p, b * p..(b + 1) * p]);
7339 block.assign(&p_k.mapv(|v| v * factor));
7340 }
7341 }
7342 s_big
7343 };
7344 let (s_null, null_scale) =
7345 normalize_penalty_in_constrained_space(&stz_pooled_null);
7346 let null_block = crate::basis::analyze_penalty_block_with_op(&s_null, None)?;
7347 if null_block.rank > 0 {
7348 let original_index = penalties.len();
7349 penalties.push(null_block.sym_penalty);
7350 nullspaces.push(null_block.nullity);
7351 inner_built.penaltyinfo.push(PenaltyInfo {
7352 source: PenaltySource::Primary,
7353 original_index,
7354 active: true,
7355 effective_rank: null_block.rank,
7356 dropped_reason: None,
7357 nullspace_dim_hint: null_block.nullity,
7358 normalization_scale: null_scale,
7359 kronecker_factors: None,
7360 });
7361 }
7362 }
7363 }
7364 inner_built.dim = p * l_minus_one;
7365 inner_built.design = DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(dense));
7366 inner_built.penalties = penalties;
7367 inner_built.ops = vec![None; inner_built.penalties.len()];
7368 inner_built.nullspaces = nullspaces;
7369 inner_built.null_eigenvectors =
7376 crate::basis::recompute_null_eigenvectors(&inner_built.penalties)?;
7377 inner_built.joint_null_rotation =
7378 crate::basis::compute_joint_null_rotation(&inner_built.penalties)?;
7379 inner_built.kronecker_factored = None;
7380 return Ok(inner_built);
7381 }
7382 SmoothBasisSpec::BSpline1D { feature_col, spec } => {
7383 if *feature_col >= data.ncols() {
7384 crate::bail_dim_basis!(
7385 "term '{}' feature column {} out of bounds for {} columns",
7386 term.name,
7387 feature_col,
7388 data.ncols()
7389 );
7390 }
7391 let mut spec_local = spec.clone();
7392 if term.shape != ShapeConstraint::None {
7393 spec_local.identifiability = BSplineIdentifiability::None;
7396 }
7397 build_bspline_basis_1d(data.column(*feature_col), &spec_local)?
7401 }
7402 SmoothBasisSpec::ThinPlate {
7403 feature_cols,
7404 spec,
7405 input_scales,
7406 } => {
7407 if term.shape != ShapeConstraint::None {
7408 if feature_cols.len() != 1 {
7409 crate::bail_invalid_basis!(
7410 "ShapeConstraint::{:?} for term '{}' on ThinPlate basis requires exactly 1 feature axis; found {}",
7411 term.shape,
7412 term.name,
7413 feature_cols.len()
7414 );
7415 }
7416 shape_axis_col = Some(feature_cols[0]);
7417 }
7418 let mut x = select_columns(data, feature_cols)?;
7419 let (scales, length_scale_eff) = if let Some(s) = input_scales {
7425 apply_input_standardization(&mut x, s);
7426 (
7427 Some(s.clone()),
7428 compensate_length_scale_for_standardization(spec.length_scale, s),
7429 )
7430 } else if let Some(s) = compute_spatial_input_scales(x.view()) {
7431 apply_input_standardization(&mut x, &s);
7432 let l_eff = compensate_length_scale_for_standardization(spec.length_scale, &s);
7433 (Some(s), l_eff)
7434 } else {
7435 (None, spec.length_scale)
7436 };
7437 let mut spec_local = spec.clone();
7438 spec_local.length_scale = length_scale_eff;
7439 if matches!(
7440 spec_local.identifiability,
7441 SpatialIdentifiability::OrthogonalToParametric
7442 ) {
7443 spec_local.identifiability = SpatialIdentifiability::None;
7444 }
7445 let mut result = build_thin_plate_basis(x.view(), &spec_local).map_err(|err| {
7446 rewrite_thin_plate_knots_error(err, &term.name, feature_cols.len(), spec)
7447 })?;
7448 match &mut result.metadata {
7456 BasisMetadata::ThinPlate {
7457 input_scales: ms,
7458 length_scale,
7459 ..
7460 } => {
7461 *ms = scales;
7462 *length_scale = spec.length_scale;
7463 }
7464 BasisMetadata::Duchon {
7465 input_scales: ms,
7466 length_scale,
7467 ..
7468 } => {
7469 if let (Some(s), Some(realized)) = (scales.as_ref(), *length_scale) {
7494 let inv_sigma_geom =
7495 compensate_length_scale_for_standardization(1.0, s);
7496 if inv_sigma_geom.is_finite() && inv_sigma_geom > 0.0 {
7497 *length_scale = Some(realized / inv_sigma_geom);
7498 }
7499 }
7500 *ms = scales;
7501 }
7502 _ => {}
7503 }
7504 result
7505 }
7506 SmoothBasisSpec::Sphere { feature_cols, spec } => {
7507 if term.shape != ShapeConstraint::None {
7508 crate::bail_invalid_basis!(
7509 "ShapeConstraint::{:?} for term '{}' is not supported on spherical splines",
7510 term.shape,
7511 term.name
7512 );
7513 }
7514 let x = select_columns(data, feature_cols)?;
7515 build_spherical_spline_basis(x.view(), spec)?
7516 }
7517 SmoothBasisSpec::ConstantCurvature { feature_cols, spec } => {
7518 if term.shape != ShapeConstraint::None {
7519 crate::bail_invalid_basis!(
7520 "ShapeConstraint::{:?} for term '{}' is not supported on constant-curvature smooths",
7521 term.shape,
7522 term.name
7523 );
7524 }
7525 let x = select_columns(data, feature_cols)?;
7532 build_constant_curvature_basis(x.view(), spec)?
7533 }
7534 SmoothBasisSpec::MeasureJet {
7535 feature_cols,
7536 spec,
7537 input_scales,
7538 } => {
7539 if term.shape != ShapeConstraint::None {
7540 crate::bail_invalid_basis!(
7541 "ShapeConstraint::{:?} for term '{}' is not supported on measure-jet smooths",
7542 term.shape,
7543 term.name
7544 );
7545 }
7546 let mut x = select_columns(data, feature_cols)?;
7547 let (scales, length_scale_eff) = if let Some(s) = input_scales {
7559 apply_input_standardization(&mut x, s);
7560 (Some(s.clone()), spec.length_scale)
7561 } else if let Some(s) = compute_spatial_input_scales(x.view()) {
7562 apply_input_standardization(&mut x, &s);
7563 let l_eff = if spec.length_scale > 0.0 {
7564 compensate_length_scale_for_standardization(spec.length_scale, &s)
7565 } else {
7566 spec.length_scale
7567 };
7568 (Some(s), l_eff)
7569 } else {
7570 (None, spec.length_scale)
7571 };
7572 let mut spec_local = spec.clone();
7573 spec_local.length_scale = length_scale_eff;
7574 let mut result = build_measure_jet_basis(x.view(), &spec_local)?;
7575 if let BasisMetadata::MeasureJet {
7576 input_scales: ms, ..
7577 } = &mut result.metadata
7578 {
7579 *ms = scales;
7580 }
7581 result
7582 }
7583 SmoothBasisSpec::Matern {
7584 feature_cols,
7585 spec,
7586 input_scales,
7587 } => {
7588 if term.shape != ShapeConstraint::None {
7589 if feature_cols.len() != 1 {
7590 crate::bail_invalid_basis!(
7591 "ShapeConstraint::{:?} for term '{}' on Matern basis requires exactly 1 feature axis; found {}",
7592 term.shape,
7593 term.name,
7594 feature_cols.len()
7595 );
7596 }
7597 shape_axis_col = Some(feature_cols[0]);
7598 }
7599 let mut x = select_columns(data, feature_cols)?;
7600 let (scales, length_scale_eff) = if let Some(s) = input_scales {
7615 apply_input_standardization(&mut x, s);
7616 (
7617 Some(s.clone()),
7618 compensate_length_scale_for_standardization(spec.length_scale, s),
7619 )
7620 } else if let Some(s) = compute_spatial_input_scales(x.view()) {
7621 apply_input_standardization(&mut x, &s);
7622 let l_eff = compensate_length_scale_for_standardization(spec.length_scale, &s);
7623 (Some(s), l_eff)
7624 } else {
7625 (None, spec.length_scale)
7626 };
7627 let mut spec_local = spec.clone();
7628 spec_local.length_scale = length_scale_eff;
7629 let mut result = build_matern_basiswithworkspace(x.view(), &spec_local, workspace)?;
7630 if let BasisMetadata::Matern {
7631 input_scales,
7632 length_scale,
7633 ..
7634 } = &mut result.metadata
7635 {
7636 *input_scales = scales;
7637 *length_scale = spec.length_scale;
7638 }
7639 result
7640 }
7641 SmoothBasisSpec::Duchon {
7642 feature_cols,
7643 spec,
7644 input_scales,
7645 } => {
7646 if term.shape != ShapeConstraint::None {
7647 if feature_cols.len() != 1 {
7648 crate::bail_invalid_basis!(
7649 "ShapeConstraint::{:?} for term '{}' on Duchon basis requires exactly 1 feature axis; found {}",
7650 term.shape,
7651 term.name,
7652 feature_cols.len()
7653 );
7654 }
7655 shape_axis_col = Some(feature_cols[0]);
7656 }
7657 let mut x = select_columns(data, feature_cols)?;
7658 let (scales, length_scale_eff) = if let Some(s) = input_scales {
7669 apply_input_standardization(&mut x, s);
7670 (
7671 Some(s.clone()),
7672 compensate_optional_length_scale_for_standardization(spec.length_scale, s),
7673 )
7674 } else if let Some(s) = compute_spatial_input_scales(x.view()) {
7675 apply_input_standardization(&mut x, &s);
7676 let l_eff =
7677 compensate_optional_length_scale_for_standardization(spec.length_scale, &s);
7678 (Some(s), l_eff)
7679 } else {
7680 (None, spec.length_scale)
7681 };
7682 let mut spec_local = spec.clone();
7683 spec_local.length_scale = length_scale_eff;
7684 if let (Some(s), crate::basis::OneDimensionalBoundary::Cyclic { start, end }) =
7695 (scales.as_ref(), spec_local.boundary.clone())
7696 && s.len() == 1
7697 && s[0] > 0.0
7698 {
7699 spec_local.boundary = crate::basis::OneDimensionalBoundary::Cyclic {
7700 start: start / s[0],
7701 end: end / s[0],
7702 };
7703 }
7704 if matches!(
7705 spec_local.identifiability,
7706 SpatialIdentifiability::OrthogonalToParametric
7707 ) {
7708 spec_local.identifiability = SpatialIdentifiability::None;
7709 }
7710 let mut result = build_duchon_basiswithworkspace(x.view(), &spec_local, workspace)?;
7711 if let BasisMetadata::Duchon {
7712 input_scales,
7713 length_scale,
7714 ..
7715 } = &mut result.metadata
7716 {
7717 *input_scales = scales;
7718 *length_scale = spec.length_scale;
7719 }
7720 result
7721 }
7722 SmoothBasisSpec::Pca {
7723 feature_cols,
7724 basis_matrix,
7725 centered,
7726 smooth_penalty,
7727 center_mean,
7728 pca_basis_path,
7729 chunk_size,
7730 } => {
7731 if term.shape != ShapeConstraint::None {
7732 crate::bail_invalid_basis!(
7733 "ShapeConstraint::{:?} for term '{}' is not supported on Pca basis",
7734 term.shape,
7735 term.name
7736 );
7737 }
7738 build_pca_smooth_basis(
7739 data,
7740 feature_cols,
7741 basis_matrix,
7742 *centered,
7743 *smooth_penalty,
7744 center_mean.as_ref(),
7745 pca_basis_path.as_ref(),
7746 *chunk_size,
7747 )?
7748 }
7749 SmoothBasisSpec::TensorBSpline { feature_cols, spec } => {
7750 build_tensor_bspline_basis(data, feature_cols, spec)?
7751 }
7752 SmoothBasisSpec::ByVariable { .. } => {
7753 crate::bail_invalid_basis!(
7754 "internal: ByVariable smooths must return before inner basis dispatch"
7755 );
7756 }
7757 SmoothBasisSpec::BySmooth { .. } => {
7758 crate::bail_invalid_basis!("internal: BySmooth smooths must be lowered to ByVariable before inner basis dispatch"
7759 .to_string(),);
7760 }
7761 SmoothBasisSpec::FactorSmooth { spec } => {
7762 if term.shape != ShapeConstraint::None {
7763 crate::bail_invalid_basis!(
7764 "ShapeConstraint::{:?} is unsupported for factor smooth term '{}'",
7765 term.shape,
7766 term.name
7767 );
7768 }
7769 return build_factor_smooth(data, spec, &term.name, workspace);
7770 }
7771 };
7772
7773 if let SmoothBasisSpec::Matern { .. } = &term.basis {
7789 let (penalties, nullspace_dims, penaltyinfo) =
7790 matern_operator_penalty_triplet_from_metadata(&built.metadata)?;
7791 built.penalties = penalties;
7792 built.nullspace_dims = nullspace_dims;
7793 built.penaltyinfo = penaltyinfo;
7794 }
7795
7796 let p_local = built.design.ncols();
7797 let mut metadata = built.metadata.clone();
7798 let kron_factored = if term.shape == ShapeConstraint::None {
7801 built.kronecker_factored
7802 } else {
7803 None
7804 };
7805 let mut design_t = built.design;
7806 let mut penalties_t: Vec<Array2<f64>> = built.penalties;
7807 let mut ops_t: Vec<Option<std::sync::Arc<dyn crate::analytic_penalties::PenaltyOp>>> =
7812 built.ops;
7813 if matches!(
7814 spatial_identifiability_policy(term),
7815 Some(SpatialIdentifiability::OrthogonalToParametric)
7816 ) {
7817 metadata = freeze_raw_spatial_metadata(metadata, design_t.ncols());
7818 }
7819
7820 let active_penaltyinfo_t = built
7821 .penaltyinfo
7822 .iter()
7823 .filter(|info| info.active)
7824 .cloned()
7825 .collect::<Vec<_>>();
7826 let pre_dropped_penaltyinfo_t = built
7827 .penaltyinfo
7828 .iter()
7829 .filter(|info| !info.active)
7830 .cloned()
7831 .collect::<Vec<_>>();
7832 let use_box_reparam =
7833 term.shape != ShapeConstraint::None && shape_uses_box_reparameterization(&term.basis);
7834 if let Some((order, sign)) = shape_order_and_sign(term.shape)
7835 && use_box_reparam
7836 {
7837 let t = if order == 2 {
7852 let bspline_meta = match &metadata {
7853 BasisMetadata::BSpline1D {
7854 knots,
7855 degree,
7856 periodic,
7857 ..
7858 } if periodic.is_none() => Some((knots.clone(), degree.unwrap_or(0))),
7859 _ => None,
7860 };
7861 match bspline_meta {
7862 Some((knots, degree)) if degree >= 1 => {
7863 let greville = crate::basis::compute_greville_abscissae(&knots, degree)?;
7864 if greville.len() != p_local {
7865 crate::bail_invalid_basis!(
7866 "shape-constraint Greville abscissae count {} does not match basis dim {} for term '{}'",
7867 greville.len(),
7868 p_local,
7869 term.name
7870 );
7871 }
7872 convex_divided_difference_transform_matrix(&greville, sign)?
7873 }
7874 _ => cumulative_sum_transform_matrix(p_local, order, sign),
7875 }
7876 } else {
7877 cumulative_sum_transform_matrix(p_local, order, sign)
7878 };
7879 let inner_dense = match design_t {
7883 DesignMatrix::Dense(d) => d,
7884 DesignMatrix::Sparse(sp) => gam_linalg::matrix::DenseDesignMatrix::from(
7885 sp.try_to_dense_arc("shape-constrained coefficient transform")
7886 .map_err(BasisError::InvalidInput)?,
7887 ),
7888 };
7889 let coeff_op = gam_linalg::matrix::CoefficientTransformOperator::new(inner_dense, t.clone())
7890 .map_err(|e| BasisError::InvalidInput(format!("CoefficientTransformOperator: {e}")))?;
7891 design_t = DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(Arc::new(coeff_op)));
7892 if penalties_t.len() != active_penaltyinfo_t.len() {
7893 crate::bail_invalid_basis!(
7894 "internal box-reparam penalty/info mismatch for term '{}': penalties={}, infos={}",
7895 term.name,
7896 penalties_t.len(),
7897 active_penaltyinfo_t.len()
7898 );
7899 }
7900 let transformed_wiggliness = penalties_t
7916 .iter()
7917 .zip(active_penaltyinfo_t.iter())
7918 .find(|(_, info)| !matches!(info.source, PenaltySource::DoublePenaltyNullspace))
7919 .map(|(s_local, _)| {
7920 let tt_s = fast_atb(&t, s_local);
7921 fast_ab(&tt_s, &t)
7922 });
7923 let mut rebuilt = Vec::with_capacity(penalties_t.len());
7924 for (s_local, info) in penalties_t.iter().zip(active_penaltyinfo_t.iter()) {
7925 if matches!(info.source, PenaltySource::DoublePenaltyNullspace) {
7926 if order == 2 {
7961 let tt_s = fast_atb(&t, s_local);
7962 rebuilt.push(fast_ab(&tt_s, &t));
7963 } else {
7964 let s_wiggle_t = transformed_wiggliness.as_ref().ok_or_else(|| {
7965 BasisError::InvalidInput(format!(
7966 "box-reparam term '{}' has a double-penalty ridge but no primary wiggliness penalty to derive its nullspace from",
7967 term.name
7968 ))
7969 })?;
7970 let ridge = crate::basis::build_nullspace_shrinkage_penalty(s_wiggle_t)?
7971 .map(|shrink| shrink.sym_penalty)
7972 .unwrap_or_else(|| Array2::<f64>::zeros((p_local, p_local)));
7973 rebuilt.push(ridge);
7974 }
7975 } else {
7976 let tt_s = fast_atb(&t, s_local);
7977 rebuilt.push(fast_ab(&tt_s, &t));
7978 }
7979 }
7980 penalties_t = rebuilt;
7981 ops_t = vec![None; penalties_t.len()];
7984 }
7985 if penalties_t.len() != active_penaltyinfo_t.len() {
7986 crate::bail_invalid_basis!(
7987 "internal penalty metadata mismatch for term '{}': active penalties={}, active infos={}",
7988 term.name,
7989 penalties_t.len(),
7990 active_penaltyinfo_t.len()
7991 );
7992 }
7993 if ops_t.len() != penalties_t.len() {
7994 ops_t = vec![None; penalties_t.len()];
7995 }
7996 let penalty_candidates = penalties_t
7997 .into_iter()
7998 .zip(active_penaltyinfo_t.into_iter())
7999 .zip(ops_t.into_iter())
8000 .map(
8001 |((matrix, info), op_in)| -> Result<PenaltyCandidate, BasisError> {
8002 let (matrix, c_new) = normalize_penalty_in_constrained_space(&matrix);
8003 let normalization_scale = info.normalization_scale * c_new;
8004 let op_scale = 1.0 / c_new;
8005 let kronecker_scale = 1.0 / c_new;
8006 let scaled_op = if op_scale > 0.0 && op_scale.is_finite() {
8009 op_in.map(|op| {
8010 std::sync::Arc::new(crate::analytic_penalties::ScaledPenaltyOp::new(
8011 op, op_scale,
8012 ))
8013 as std::sync::Arc<dyn crate::analytic_penalties::PenaltyOp>
8014 })
8015 } else {
8016 None
8017 };
8018 let kronecker_factors = info.kronecker_factors.map(|mut factors| {
8019 if let Some(first) = factors.first_mut() {
8020 first.mapv_inplace(|v| v * kronecker_scale);
8021 }
8022 factors
8023 });
8024 Ok(PenaltyCandidate {
8025 nullspace_dim_hint: info.nullspace_dim_hint,
8026 matrix,
8027 source: info.source,
8028 normalization_scale,
8029 kronecker_factors,
8030 op: scaled_op,
8031 })
8032 },
8033 )
8034 .collect::<Result<Vec<_>, _>>()?;
8035 let (penalties_t, nullspaces_t, penaltyinfo_t, null_eigenvectors_t, ops_t) =
8036 crate::basis::filter_active_penalty_candidates_with_ops(penalty_candidates)?;
8037 let shape_linear_constraints = if term.shape != ShapeConstraint::None && !use_box_reparam {
8038 let axis = shape_axis_col.ok_or_else(|| {
8039 BasisError::InvalidInput(format!(
8040 "internal shape-constraint axis missing for term '{}'",
8041 term.name
8042 ))
8043 })?;
8044 let (x_shape_eval, design_shape_eval) =
8045 build_shape_constraint_design_1d(data, term, &metadata, axis)?;
8046 build_shape_linear_constraints_1d(
8047 x_shape_eval.view(),
8048 design_shape_eval.view(),
8049 term.shape,
8050 )?
8051 } else {
8052 None
8053 };
8054 let linear_constraints_local = merge_linear_constraints_global(shape_linear_constraints, None);
8055
8056 let joint_null_rotation = match term.joint_null_rotation.clone() {
8075 Some(persisted) => Some(persisted),
8076 None if smooth_has_frozen_identifiability(term) => None,
8077 None if kron_factored.is_some() => None,
8078 None => crate::basis::compute_joint_null_rotation(&penalties_t)?,
8079 };
8080
8081 Ok(LocalSmoothTermBuild {
8082 dim: p_local,
8083 design: design_t,
8084 penalties: penalties_t,
8085 ops: ops_t,
8086 nullspaces: nullspaces_t,
8087 null_eigenvectors: null_eigenvectors_t,
8088 joint_null_rotation,
8089 penaltyinfo: penaltyinfo_t,
8090 pre_dropped_penaltyinfo: pre_dropped_penaltyinfo_t,
8091 metadata,
8092 linear_constraints: linear_constraints_local,
8093 box_reparam: use_box_reparam,
8094 kronecker_factored: kron_factored,
8095 })
8096}
8097
8098pub fn build_smooth_design(
8099 data: ArrayView2<'_, f64>,
8100 terms: &[SmoothTermSpec],
8101) -> Result<RawSmoothDesign, BasisError> {
8102 let mut ws = crate::basis::BasisWorkspace::new();
8103 build_smooth_design_withworkspace(data, terms, &mut ws)
8104}
8105
8106pub fn build_smooth_design_withworkspace(
8113 data: ArrayView2<'_, f64>,
8114 terms: &[SmoothTermSpec],
8115 workspace: &mut crate::basis::BasisWorkspace,
8116) -> Result<RawSmoothDesign, BasisError> {
8117 validate_smooth_terms_finite_inputs(data, terms)?;
8118 build_smooth_design_withworkspace_unvalidated(data, terms, workspace)
8119}
8120
8121pub fn build_smooth_design_withworkspace_unvalidated(
8122 data: ArrayView2<'_, f64>,
8123 terms: &[SmoothTermSpec],
8124 workspace: &mut crate::basis::BasisWorkspace,
8125) -> Result<RawSmoothDesign, BasisError> {
8126 let mut planned_blocks = plan_joint_spatial_centers_for_term_blocks(data, &[terms.to_vec()])?;
8127 let planned_terms = planned_blocks.pop().ok_or_else(|| {
8128 BasisError::InvalidInput(
8129 "joint spatial center planner returned no smooth blocks".to_string(),
8130 )
8131 })?;
8132 let policy = workspace.policy().clone();
8133 let local_builds: Vec<LocalSmoothTermBuild> = {
8134 use rayon::iter::{IntoParallelIterator, ParallelIterator};
8135 planned_terms
8136 .into_par_iter()
8137 .map(|term| {
8138 let mut term_workspace = crate::basis::BasisWorkspace::with_policy(policy.clone());
8139 build_single_local_smooth_term(data, &term, &mut term_workspace)
8140 })
8141 .collect::<Result<Vec<_>, _>>()?
8142 };
8143
8144 let total_p: usize = local_builds.iter().map(|built| built.dim).sum();
8145
8146 let mut local_designs: Vec<DesignMatrix> = Vec::with_capacity(local_builds.len());
8147 let mut terms_out = Vec::<SmoothTerm>::with_capacity(terms.len());
8148 let mut penalties_global = Vec::<BlockwisePenalty>::new();
8149 let mut nullspace_dims_global = Vec::<usize>::new();
8150 let mut penaltyinfo_global = Vec::<PenaltyBlockInfo>::new();
8151 let mut dropped_penaltyinfo_global = Vec::<DroppedPenaltyBlockInfo>::new();
8152 let mut coefficient_lower_bounds = Array1::<f64>::from_elem(total_p, f64::NEG_INFINITY);
8153 let mut any_bounds = false;
8154 let mut linear_constraintsrows: Vec<(usize, usize, Array1<f64>)> = Vec::new();
8159 let mut linear_constraints_b: Vec<f64> = Vec::new();
8160
8161 let mut col_start = 0usize;
8162 for (term, mut built) in terms.iter().zip(local_builds.into_iter()) {
8163 let p_local = built.dim;
8164 let col_end = col_start + p_local;
8165 let lb_local = if built.box_reparam {
8166 shape_lower_bounds_local(term.shape, p_local)
8167 } else {
8168 None
8169 };
8170
8171 let applied_rotation: Option<crate::basis::JointNullRotation> = match (
8203 built.joint_null_rotation.take(),
8204 lb_local.is_some(),
8205 built.linear_constraints.is_some(),
8206 ) {
8207 (Some(rot), false, false) => {
8208 let q = &rot.rotation;
8209 let dense = built
8210 .design
8211 .try_to_dense_by_chunks("joint-null absorption rotation")
8212 .map_err(BasisError::InvalidInput)?;
8213 let rotated = gam_linalg::faer_ndarray::fast_ab(&dense, q);
8214 built.design = DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(rotated));
8215 built.penalties = built
8216 .penalties
8217 .into_iter()
8218 .map(|s_local| {
8219 let qt_s = gam_linalg::faer_ndarray::fast_atb(q, &s_local);
8220 gam_linalg::faer_ndarray::fast_ab(&qt_s, q)
8221 })
8222 .collect();
8223 built.ops = vec![None; built.penalties.len()];
8224 built.kronecker_factored = None;
8225 Some(rot)
8226 }
8227 (Some(_), _, _) => None,
8228 (None, _, _) => None,
8229 };
8230
8231 let activeinfos = built
8232 .penaltyinfo
8233 .iter()
8234 .filter(|info| info.active)
8235 .collect::<Vec<_>>();
8236 if activeinfos.len() != built.penalties.len() {
8237 crate::bail_invalid_basis!(
8238 "internal penalty info mismatch for term '{}': activeinfos={}, penalties={}",
8239 term.name,
8240 activeinfos.len(),
8241 built.penalties.len()
8242 );
8243 }
8244 for (((s_local, &ns), info), op_local) in built
8245 .penalties
8246 .iter()
8247 .zip(built.nullspaces.iter())
8248 .zip(activeinfos.into_iter())
8249 .zip(built.ops.iter())
8250 {
8251 let global_index = penalties_global.len();
8252 penalties_global.push(
8253 BlockwisePenalty::new(col_start..col_end, s_local.clone())
8254 .with_op(op_local.clone()),
8255 );
8256 nullspace_dims_global.push(ns);
8257 let mut penalty = info.clone();
8258 penalty.nullspace_dim_hint = ns;
8259 penaltyinfo_global.push(PenaltyBlockInfo {
8260 global_index,
8261 termname: Some(term.name.clone()),
8262 penalty,
8263 });
8264 }
8265 for info in built.penaltyinfo.iter().filter(|info| !info.active) {
8266 dropped_penaltyinfo_global.push(DroppedPenaltyBlockInfo {
8267 termname: Some(term.name.clone()),
8268 penalty: info.clone(),
8269 });
8270 }
8271 for info in &built.pre_dropped_penaltyinfo {
8272 dropped_penaltyinfo_global.push(DroppedPenaltyBlockInfo {
8273 termname: Some(term.name.clone()),
8274 penalty: info.clone(),
8275 });
8276 }
8277
8278 if let Some(lin_local) = &built.linear_constraints {
8279 for r in 0..lin_local.a.nrows() {
8280 linear_constraintsrows.push((col_start, col_end, lin_local.a.row(r).to_owned()));
8281 linear_constraints_b.push(lin_local.b[r]);
8282 }
8283 }
8284 if let Some(lb_local) = &lb_local {
8285 coefficient_lower_bounds
8286 .slice_mut(s![col_start..col_end])
8287 .assign(lb_local);
8288 any_bounds = true;
8289 }
8290
8291 local_designs.push(built.design);
8293
8294 terms_out.push(SmoothTerm {
8295 name: term.name.clone(),
8296 coeff_range: col_start..col_end,
8297 shape: term.shape,
8298 penalties_local: built.penalties,
8299 nullspace_dims: built.nullspaces,
8300 penaltyinfo_local: built.penaltyinfo,
8301 metadata: built.metadata,
8302 lower_bounds_local: lb_local,
8303 linear_constraints_local: built.linear_constraints,
8304 kronecker_factored: built.kronecker_factored.take(),
8305 joint_null_rotation: applied_rotation,
8306 unabsorbed_global_orthogonality: None,
8307 });
8308
8309 col_start = col_end;
8310 }
8311
8312 assert_eq!(
8313 penalties_global.len(),
8314 nullspace_dims_global.len(),
8315 "global smooth penalty/nullspace bookkeeping diverged"
8316 );
8317 assert_eq!(
8318 penalties_global.len(),
8319 penaltyinfo_global.len(),
8320 "global smooth penalty metadata bookkeeping diverged"
8321 );
8322
8323 Ok(RawSmoothDesign {
8324 term_designs: local_designs,
8325 penalties: penalties_global,
8326 nullspace_dims: nullspace_dims_global,
8327 penaltyinfo: penaltyinfo_global,
8328 dropped_penaltyinfo: dropped_penaltyinfo_global,
8329 terms: terms_out,
8330 coefficient_lower_bounds: if any_bounds {
8331 Some(coefficient_lower_bounds)
8332 } else {
8333 None
8334 },
8335 linear_constraints: if linear_constraintsrows.is_empty() {
8336 None
8337 } else {
8338 let mut a = Array2::<f64>::zeros((linear_constraintsrows.len(), total_p));
8339 for (i, (cs, ce, values)) in linear_constraintsrows.iter().enumerate() {
8340 a.row_mut(i).slice_mut(s![*cs..*ce]).assign(values);
8341 }
8342 Some(LinearInequalityConstraints {
8343 a,
8344 b: Array1::from_vec(linear_constraints_b),
8345 })
8346 },
8347 })
8348}