1pub fn build_term_collection_designs_joint(
9 data: ArrayView2<'_, f64>,
10 specs: &[TermCollectionSpec],
11) -> Result<Vec<TermCollectionDesign>, BasisError> {
12 for spec in specs {
13 validate_term_collection_finite_inputs(data, spec)?;
14 }
15 let smooth_blocks = specs
16 .iter()
17 .map(|spec| spec.smooth_terms.clone())
18 .collect::<Vec<_>>();
19 let planned_blocks = plan_joint_spatial_centers_for_term_blocks(data, &smooth_blocks)?;
20 let mut out = Vec::with_capacity(specs.len());
21 for (spec, planned_terms) in specs.iter().zip(planned_blocks.into_iter()) {
22 let mut planned_spec = spec.clone();
23 planned_spec.smooth_terms = planned_terms;
24 out.push(build_term_collection_design_inner(data, &planned_spec)?);
25 }
26 Ok(out)
27}
28
29pub fn build_term_collection_designs_and_freeze_joint(
30 data: ArrayView2<'_, f64>,
31 specs: &[TermCollectionSpec],
32) -> Result<(Vec<TermCollectionDesign>, Vec<TermCollectionSpec>), EstimationError> {
33 let designs = build_term_collection_designs_joint(data, specs)?;
34 let mut resolved_specs = Vec::with_capacity(specs.len());
35 for (spec, design) in specs.iter().zip(designs.iter()) {
36 resolved_specs.push(freeze_term_collection_from_design(spec, design)?);
37 }
38 Ok((designs, resolved_specs))
39}
40
41pub fn fit_term_collection_forspec(
42 data: ArrayView2<'_, f64>,
43 y: ArrayView1<'_, f64>,
44 weights: ArrayView1<'_, f64>,
45 offset: ArrayView1<'_, f64>,
46 spec: &TermCollectionSpec,
47 family: LikelihoodSpec,
48 options: &FitOptions,
49) -> Result<FittedTermCollection, EstimationError> {
50 fit_term_collection_forspecwith_heuristic_lambdas(
51 data, y, weights, offset, spec, None, family, options,
52 )
53}
54
55pub fn fit_term_collection_with_coefficient_groups(
56 data: ArrayView2<'_, f64>,
57 y: ArrayView1<'_, f64>,
58 weights: ArrayView1<'_, f64>,
59 offset: ArrayView1<'_, f64>,
60 spec: &TermCollectionSpec,
61 groups: &[CoefficientGroupSpec],
62 family: LikelihoodSpec,
63 options: &FitOptions,
64) -> Result<FittedTermCollection, EstimationError> {
65 if groups.is_empty() {
66 return fit_term_collection_forspec(data, y, weights, offset, spec, family, options);
67 }
68 let design = build_term_collection_design_with_policy(data, spec, &options.resource_policy)?;
69 let base_fit_opts = adaptive_fit_options_base(options, &design);
70 let realized = design
71 .realize_coefficient_groups(groups, &base_fit_opts.rho_prior)
72 .map_err(EstimationError::BasisError)?;
73 let effective_offset = design
74 .compose_offset(offset, "coefficient-group fit")
75 .map_err(EstimationError::BasisError)?;
76 let mut grouped_options = base_fit_opts.clone();
77 grouped_options.rho_prior = realized.rho_prior;
78 let fitted = FittedTermCollection {
79 fit: gam_solve::estimate::fit_gam_with_penalty_specs(
80 design.design.clone(),
81 y,
82 weights,
83 effective_offset.view(),
84 realized.penalty_specs,
85 realized.nullspace_dims,
86 family.clone(),
87 &grouped_options,
88 )?,
89 design,
90 adaptive_diagnostics: None,
91 };
92 enforce_term_constraint_feasibility(&fitted.design, &fitted.fit)?;
93 Ok(fitted)
94}
95
96pub fn fit_term_collection_with_penalty_block_gamma_prior_callback<F>(
97 data: ArrayView2<'_, f64>,
98 y: ArrayView1<'_, f64>,
99 weights: ArrayView1<'_, f64>,
100 offset: ArrayView1<'_, f64>,
101 spec: &TermCollectionSpec,
102 callback: F,
103 family: LikelihoodSpec,
104 options: &FitOptions,
105) -> Result<FittedTermCollection, EstimationError>
106where
107 F: FnMut(&PenaltyBlockGammaPriorMetadata<'_>) -> Option<(f64, f64)>,
108{
109 let design = build_term_collection_design_with_policy(data, spec, &options.resource_policy)?;
110 let effective_offset = design
111 .compose_offset(offset, "penalty-prior callback fit")
112 .map_err(EstimationError::BasisError)?;
113 let mut fit_opts = adaptive_fit_options_base(options, &design);
114 fit_opts.rho_prior = realize_penalty_block_gamma_priors(&design, callback)
115 .map_err(EstimationError::BasisError)?;
116 let fitted = FittedTermCollection {
117 fit: fit_gamwith_heuristic_lambdas(
118 design.design.clone(),
119 y,
120 weights,
121 effective_offset.view(),
122 &design.penalties,
123 None,
124 family.clone(),
125 &fit_opts,
126 )?,
127 design,
128 adaptive_diagnostics: None,
129 };
130 enforce_term_constraint_feasibility(&fitted.design, &fitted.fit)?;
131 Ok(fitted)
132}
133
134pub fn fit_term_collection_with_penalty_block_gamma_priors(
135 data: ArrayView2<'_, f64>,
136 y: ArrayView1<'_, f64>,
137 weights: ArrayView1<'_, f64>,
138 offset: ArrayView1<'_, f64>,
139 spec: &TermCollectionSpec,
140 priors: &[(String, f64, f64)],
141 family: LikelihoodSpec,
142 options: &FitOptions,
143) -> Result<FittedTermCollection, EstimationError> {
144 let design = build_term_collection_design_with_policy(data, spec, &options.resource_policy)?;
145 let effective_offset = design
146 .compose_offset(offset, "penalty-prior fit")
147 .map_err(EstimationError::BasisError)?;
148 let mut fit_opts = adaptive_fit_options_base(options, &design);
149 fit_opts.rho_prior = realize_keyed_penalty_block_gamma_priors(&design, priors)
150 .map_err(EstimationError::BasisError)?;
151 let fitted = FittedTermCollection {
152 fit: fit_gamwith_heuristic_lambdas(
153 design.design.clone(),
154 y,
155 weights,
156 effective_offset.view(),
157 &design.penalties,
158 None,
159 family.clone(),
160 &fit_opts,
161 )?,
162 design,
163 adaptive_diagnostics: None,
164 };
165 enforce_term_constraint_feasibility(&fitted.design, &fitted.fit)?;
166 Ok(fitted)
167}
168
169pub fn fit_term_collection_with_coefficient_groups_and_penalty_block_gamma_priors(
170 data: ArrayView2<'_, f64>,
171 y: ArrayView1<'_, f64>,
172 weights: ArrayView1<'_, f64>,
173 offset: ArrayView1<'_, f64>,
174 spec: &TermCollectionSpec,
175 groups: &[CoefficientGroupSpec],
176 priors: &[(String, f64, f64)],
177 family: LikelihoodSpec,
178 options: &FitOptions,
179) -> Result<FittedTermCollection, EstimationError> {
180 if groups.is_empty() {
181 return fit_term_collection_with_penalty_block_gamma_priors(
182 data, y, weights, offset, spec, priors, family, options,
183 );
184 }
185 if priors.is_empty() {
186 return fit_term_collection_with_coefficient_groups(
187 data, y, weights, offset, spec, groups, family, options,
188 );
189 }
190
191 let design = build_term_collection_design_with_policy(data, spec, &options.resource_policy)?;
195 let base_fit_opts = adaptive_fit_options_base(options, &design);
196 let base_rho_prior = realize_keyed_penalty_block_gamma_priors(&design, priors)
197 .map_err(EstimationError::BasisError)?;
198 let realized = design
199 .realize_coefficient_groups(groups, &base_rho_prior)
200 .map_err(EstimationError::BasisError)?;
201 let effective_offset = design
202 .compose_offset(offset, "coefficient-group and penalty-prior fit")
203 .map_err(EstimationError::BasisError)?;
204 let mut grouped_options = base_fit_opts.clone();
205 grouped_options.rho_prior = realized.rho_prior;
206 let fitted = FittedTermCollection {
207 fit: gam_solve::estimate::fit_gam_with_penalty_specs(
208 design.design.clone(),
209 y,
210 weights,
211 effective_offset.view(),
212 realized.penalty_specs,
213 realized.nullspace_dims,
214 family.clone(),
215 &grouped_options,
216 )?,
217 design,
218 adaptive_diagnostics: None,
219 };
220 enforce_term_constraint_feasibility(&fitted.design, &fitted.fit)?;
221 Ok(fitted)
222}
223
224fn fit_term_collection_forspecwith_heuristic_lambdas(
225 data: ArrayView2<'_, f64>,
226 y: ArrayView1<'_, f64>,
227 weights: ArrayView1<'_, f64>,
228 offset: ArrayView1<'_, f64>,
229 spec: &TermCollectionSpec,
230 heuristic_lambdas: Option<&[f64]>,
231 family: LikelihoodSpec,
232 options: &FitOptions,
233) -> Result<FittedTermCollection, EstimationError> {
234 let adaptive_opts = options.adaptive_regularization.clone().unwrap_or_default();
235 let resolved_spec;
236 let design_spec = if adaptive_opts.enabled {
237 resolved_spec = ensure_matern_adaptive_center_resolution(spec, data.nrows());
238 &resolved_spec
239 } else {
240 spec
241 };
242 let base_design =
243 build_term_collection_design_with_policy(data, design_spec, &options.resource_policy)?;
244 fit_term_collection_on_realized_design(
245 y,
246 weights,
247 offset,
248 design_spec,
249 &base_design,
250 heuristic_lambdas,
251 family,
252 options,
253 )
254}
255
256fn ensure_matern_adaptive_center_resolution(
257 spec: &TermCollectionSpec,
258 n_rows: usize,
259) -> TermCollectionSpec {
260 let mut out = spec.clone();
261 for term in &mut out.smooth_terms {
262 let gam_terms::smooth::SmoothBasisSpec::Matern {
263 feature_cols,
264 spec: matern,
265 ..
266 } = &mut term.basis
267 else {
268 continue;
269 };
270 if let gam_terms::basis::CenterStrategy::FarthestPoint { num_centers } =
271 &mut matern.center_strategy
272 {
273 let min_centers = (4 * feature_cols.len()).min(n_rows).max(*num_centers);
286 *num_centers = min_centers;
287 }
288 }
289 out
290}
291
292fn has_bounded_linear_terms(spec: &TermCollectionSpec) -> bool {
293 spec.linear_terms.iter().any(|term| {
294 matches!(
295 term.coefficient_geometry,
296 LinearCoefficientGeometry::Bounded { .. }
297 )
298 })
299}
300
301fn fit_term_collection_on_realized_design(
302 y: ArrayView1<'_, f64>,
303 weights: ArrayView1<'_, f64>,
304 offset: ArrayView1<'_, f64>,
305 spec: &TermCollectionSpec,
306 design: &TermCollectionDesign,
307 heuristic_lambdas: Option<&[f64]>,
308 family: LikelihoodSpec,
309 options: &FitOptions,
310) -> Result<FittedTermCollection, EstimationError> {
311 let effective_offset = design
312 .compose_offset(offset, "term-collection fit")
313 .map_err(EstimationError::BasisError)?;
314 let offset = effective_offset.view();
315 if has_bounded_linear_terms(spec) {
316 return fit_bounded_term_collection_with_design(
317 y,
318 weights,
319 offset,
320 spec,
321 design,
322 heuristic_lambdas,
323 family,
324 options,
325 );
326 }
327 let mut base_fit_opts = adaptive_fit_options_base(options, design);
328 base_fit_opts.rho_prior = relax_smoothing_rho_prior(options, design, y, weights);
335 let fitted = FittedTermCollection {
336 fit: fit_gamwith_heuristic_lambdas(
337 design.design.clone(),
338 y,
339 weights,
340 offset,
341 &design.penalties,
342 heuristic_lambdas,
343 family.clone(),
344 &base_fit_opts,
345 )?,
346 design: design.clone(),
347 adaptive_diagnostics: None,
348 };
349 enforce_term_constraint_feasibility(&fitted.design, &fitted.fit)?;
350
351 let adaptive_opts = options.adaptive_regularization.clone().unwrap_or_default();
352 if !adaptive_opts.enabled {
353 return Ok(fitted);
354 }
355 let runtime_caches = extract_spatial_operator_runtime_caches(spec, &fitted.design)?;
356 if runtime_caches.is_empty() {
357 return Ok(fitted);
358 }
359 fit_term_collectionwith_exact_spatial_adaptive_regularization(
366 fitted,
367 y,
368 weights,
369 offset,
370 family,
371 options,
372 &runtime_caches,
373 )
374}
375
376#[derive(Clone)]
377struct SpatialOperatorRuntimeCache {
378 termname: String,
379 feature_cols: Vec<usize>,
380 coeff_global_range: Range<usize>,
381 mass_penalty_global_idx: usize,
382 tension_penalty_global_idx: usize,
383 stiffness_penalty_global_idx: usize,
384 d0: Array2<f64>,
385 d1: Array2<f64>,
386 d2: Array2<f64>,
387 collocation_points: Array2<f64>,
388 dimension: usize,
389}
390
391#[derive(Clone)]
392struct SpatialAdaptiveWeights {
393 inv_magweight: Array1<f64>,
394 invgradweight: Array1<f64>,
395 inv_lapweight: Array1<f64>,
396}
397
398#[derive(Clone)]
399struct CharbonnierScalarBlockState {
400 signal: Array1<f64>,
401 radius: Array1<f64>,
402 epsilon: f64,
403}
404
405impl CharbonnierScalarBlockState {
406 fn from_signal(signal: Array1<f64>, epsilon: f64) -> Self {
407 let eps = epsilon.max(1e-12);
408 let radius = signal.mapv(|t| (t * t + eps * eps).sqrt());
409 Self {
410 signal,
411 radius,
412 epsilon: eps,
413 }
414 }
415
416 fn absolute_signal(&self) -> Array1<f64> {
417 self.signal.mapv(f64::abs)
418 }
419
420 fn penalty_value(&self) -> f64 {
421 self.radius.iter().map(|r| r - self.epsilon).sum::<f64>()
422 }
423
424 fn betagradient_coeff(&self) -> Array1<f64> {
425 Array1::from_iter(
426 self.signal
427 .iter()
428 .zip(self.radius.iter())
429 .map(|(t, r)| t / r),
430 )
431 }
432
433 fn betahessian_diag(&self) -> Array1<f64> {
434 let eps2 = self.epsilon * self.epsilon;
435 self.radius.mapv(|r| eps2 / r.powi(3))
436 }
437
438 fn log_epsilon_gradient_terms(&self) -> Array1<f64> {
439 let epsilon = self.epsilon;
440 let eps2 = epsilon * epsilon;
441 self.radius.mapv(|r| eps2 / r - epsilon)
442 }
443
444 fn log_epsilon_betagradient_coeff(&self) -> Array1<f64> {
445 let eps2 = self.epsilon * self.epsilon;
446 Array1::from_iter(
447 self.signal
448 .iter()
449 .zip(self.radius.iter())
450 .map(|(t, r)| -eps2 * t / r.powi(3)),
451 )
452 }
453
454 fn log_epsilon_hessian_terms(&self) -> Array1<f64> {
455 let epsilon = self.epsilon;
456 let eps2 = epsilon * epsilon;
457 let eps4 = eps2 * eps2;
458 self.radius
459 .mapv(|r| 2.0 * eps2 / r - eps4 / r.powi(3) - epsilon)
460 }
461
462 fn surrogateweights_posterior_snr(
463 &self,
464 variance: &Array1<f64>,
465 weight_floor: f64,
466 weight_ceiling: f64,
467 ) -> (Array1<f64>, Array1<f64>) {
468 let eps2 = self.epsilon * self.epsilon;
526 let weight = Array1::from_iter(self.signal.iter().zip(variance.iter()).map(|(&t, &v)| {
527 let credible2 = (t * t - v.max(0.0)).max(0.0);
528 let r = (credible2 + eps2).sqrt();
529 (1.0 / r).clamp(weight_floor, weight_ceiling)
530 }));
531 let invweight = weight.mapv(|u| 1.0 / u);
532 (weight, invweight)
533 }
534
535 fn directionalhessian_diag(&self, direction_signal: &Array1<f64>) -> Array1<f64> {
536 let eps2 = self.epsilon * self.epsilon;
551 Array1::from_iter(
552 self.signal
553 .iter()
554 .zip(direction_signal.iter())
555 .zip(self.radius.iter())
556 .map(|((t, q), r)| -3.0 * eps2 * t * q / r.powi(5)),
557 )
558 }
559
560 fn second_directionalhessian_diag(
567 &self,
568 direction1_signal: &Array1<f64>,
569 direction2_signal: &Array1<f64>,
570 ) -> Array1<f64> {
571 let eps2 = self.epsilon * self.epsilon;
572 Array1::from_iter(
573 self.signal
574 .iter()
575 .zip(direction1_signal.iter())
576 .zip(direction2_signal.iter())
577 .zip(self.radius.iter())
578 .map(|(((t, q1), q2), r)| {
579 let r2 = r * r;
580 let psi4 = -3.0 * eps2 / r.powi(5) + 15.0 * eps2 * t * t / (r.powi(5) * r2);
581 psi4 * q1 * q2
582 }),
583 )
584 }
585
586 fn log_epsilon_betahessian_diag(&self) -> Array1<f64> {
587 let eps2 = self.epsilon * self.epsilon;
588 let eps4 = eps2 * eps2;
589 Array1::from_iter(
590 self.signal
591 .iter()
592 .zip(self.radius.iter())
593 .map(|(_, r)| 2.0 * eps2 / r.powi(3) - 3.0 * eps4 / r.powi(5)),
594 )
595 }
596
597 fn log_epsilon_beta_mixed_second_coeff(&self) -> Array1<f64> {
598 let eps2 = self.epsilon * self.epsilon;
599 Array1::from_iter(
600 self.signal
601 .iter()
602 .zip(self.radius.iter())
603 .map(|(t, r)| eps2 * t * (eps2 - 2.0 * t * t) / r.powi(5)),
604 )
605 }
606
607 fn log_epsilon_betahessian_second_diag(&self) -> Array1<f64> {
608 let eps2 = self.epsilon * self.epsilon;
609 let eps4 = eps2 * eps2;
610 let eps6 = eps4 * eps2;
611 Array1::from_iter(
612 self.radius.iter().map(|r| {
613 4.0 * eps2 / r.powi(3) - 18.0 * eps4 / r.powi(5) + 15.0 * eps6 / r.powi(7)
614 }),
615 )
616 }
617
618 fn log_epsilon_betahessian_directional_diag(
619 &self,
620 direction_signal: &Array1<f64>,
621 ) -> Array1<f64> {
622 let eps2 = self.epsilon * self.epsilon;
623 let eps4 = eps2 * eps2;
624 Array1::from_iter(
625 self.signal
626 .iter()
627 .zip(direction_signal.iter())
628 .zip(self.radius.iter())
629 .map(|((t, q), r)| (-6.0 * eps2 * t / r.powi(5) + 15.0 * eps4 * t / r.powi(7)) * q),
630 )
631 }
632}
633
634#[derive(Clone)]
635struct CharbonnierGroupedBlockState {
636 norm: Array1<f64>,
637 radius: Array1<f64>,
638 signal_blocks: Array2<f64>,
639 epsilon: f64,
640}
641
642impl CharbonnierGroupedBlockState {
643 fn from_signal_blocks(signal_blocks: Array2<f64>, epsilon: f64) -> Self {
644 let eps = epsilon.max(1e-12);
645 let norm = Array1::from_iter(
646 signal_blocks
647 .rows()
648 .into_iter()
649 .map(|row| row.iter().map(|v| v * v).sum::<f64>().sqrt()),
650 );
651 let radius = norm.mapv(|g| (g * g + eps * eps).sqrt());
652 Self {
653 norm,
654 radius,
655 signal_blocks,
656 epsilon: eps,
657 }
658 }
659
660 fn penalty_value(&self) -> f64 {
661 self.radius.iter().map(|r| r - self.epsilon).sum::<f64>()
662 }
663
664 fn norm_signal(&self) -> Array1<f64> {
665 self.norm.clone()
666 }
667
668 fn betagradient_blocks(&self) -> Array2<f64> {
669 let mut out = self.signal_blocks.clone();
670 for (k, mut row) in out.rows_mut().into_iter().enumerate() {
671 let scale = 1.0 / self.radius[k];
672 row.mapv_inplace(|v| v * scale);
673 }
674 out
675 }
676
677 fn betahessian_blocks(&self) -> Vec<Array2<f64>> {
678 let mut out = Vec::with_capacity(self.signal_blocks.nrows());
679 for (k, row) in self.signal_blocks.rows().into_iter().enumerate() {
680 let dim = row.len();
681 let mut block = Array2::<f64>::eye(dim);
682 block.mapv_inplace(|v| v / self.radius[k]);
683 for i in 0..dim {
684 for j in 0..dim {
685 block[[i, j]] -= row[i] * row[j] / self.radius[k].powi(3);
686 }
687 }
688 out.push(block);
689 }
690 out
691 }
692
693 fn log_epsilon_gradient_terms(&self) -> Array1<f64> {
694 let epsilon = self.epsilon;
695 let eps2 = epsilon * epsilon;
696 self.radius.mapv(|r| eps2 / r - epsilon)
697 }
698
699 fn log_epsilon_betagradient_blocks(&self) -> Array2<f64> {
700 let mut out = self.signal_blocks.clone();
701 let eps2 = self.epsilon * self.epsilon;
702 for (k, mut row) in out.rows_mut().into_iter().enumerate() {
703 let scale = -eps2 / self.radius[k].powi(3);
704 row.mapv_inplace(|v| v * scale);
705 }
706 out
707 }
708
709 fn log_epsilon_hessian_terms(&self) -> Array1<f64> {
710 let epsilon = self.epsilon;
711 let eps2 = epsilon * epsilon;
712 let eps4 = eps2 * eps2;
713 self.radius
714 .mapv(|r| 2.0 * eps2 / r - eps4 / r.powi(3) - epsilon)
715 }
716
717 fn surrogateweights_posterior_snr(
718 &self,
719 variance: &Array1<f64>,
720 weight_floor: f64,
721 weight_ceiling: f64,
722 ) -> (Array1<f64>, Array1<f64>) {
723 let eps2 = self.epsilon * self.epsilon;
765 let weight = Array1::from_iter(self.norm.iter().zip(variance.iter()).map(|(&g, &v)| {
766 let credible2 = (g * g - v.max(0.0)).max(0.0);
767 let r = (credible2 + eps2).sqrt();
768 (1.0 / r).clamp(weight_floor, weight_ceiling)
769 }));
770 let invweight = weight.mapv(|u| 1.0 / u);
771 (weight, invweight)
772 }
773
774 fn directionalhessian_blocks(&self, direction_blocks: &Array2<f64>) -> Vec<Array2<f64>> {
775 let mut out = Vec::with_capacity(self.signal_blocks.nrows());
800 for (k, (v, q)) in self
801 .signal_blocks
802 .rows()
803 .into_iter()
804 .zip(direction_blocks.rows().into_iter())
805 .enumerate()
806 {
807 let dim = v.len();
808 let dot = v.iter().zip(q.iter()).map(|(a, b)| a * b).sum::<f64>();
809 let r3 = self.radius[k].powi(3);
810 let r5 = self.radius[k].powi(5);
811 let mut block = Array2::<f64>::eye(dim);
812 block.mapv_inplace(|x| -dot * x / r3);
813 for i in 0..dim {
814 for j in 0..dim {
815 block[[i, j]] -= (q[i] * v[j] + v[i] * q[j]) / r3;
816 block[[i, j]] += 3.0 * dot * v[i] * v[j] / r5;
817 }
818 }
819 out.push(block);
820 }
821 out
822 }
823
824 fn second_directionalhessian_blocks(
841 &self,
842 direction1_blocks: &Array2<f64>,
843 direction2_blocks: &Array2<f64>,
844 ) -> Vec<Array2<f64>> {
845 let mut out = Vec::with_capacity(self.signal_blocks.nrows());
846 for ((k, v), (a, b)) in self.signal_blocks.rows().into_iter().enumerate().zip(
847 direction1_blocks
848 .rows()
849 .into_iter()
850 .zip(direction2_blocks.rows().into_iter()),
851 ) {
852 let dim = v.len();
853 let dot = |x: ndarray::ArrayView1<'_, f64>, y: ndarray::ArrayView1<'_, f64>| {
854 x.iter().zip(y.iter()).map(|(p, q)| p * q).sum::<f64>()
855 };
856 let sa = dot(v, a);
857 let sb = dot(v, b);
858 let ab = dot(a, b);
859 let r = self.radius[k];
860 let r3 = r.powi(3);
861 let r5 = r.powi(5);
862 let r7 = r5 * r * r;
863 let diag = -ab / r3 + 3.0 * sa * sb / r5;
864 let mut block = Array2::<f64>::eye(dim);
865 block.mapv_inplace(|x| diag * x);
866 for i in 0..dim {
867 for j in 0..dim {
868 block[[i, j]] -= (a[i] * b[j] + b[i] * a[j]) / r3;
869 block[[i, j]] += 3.0 * sb * (a[i] * v[j] + v[i] * a[j]) / r5;
870 block[[i, j]] += 3.0 * ab * v[i] * v[j] / r5;
871 block[[i, j]] += 3.0 * sa * (b[i] * v[j] + v[i] * b[j]) / r5;
872 block[[i, j]] -= 15.0 * sa * sb * v[i] * v[j] / r7;
873 }
874 }
875 out.push(block);
876 }
877 out
878 }
879
880 fn log_epsilon_betahessian_blocks(&self) -> Vec<Array2<f64>> {
881 let mut out = Vec::with_capacity(self.signal_blocks.nrows());
882 for (k, row) in self.signal_blocks.rows().into_iter().enumerate() {
883 let dim = row.len();
884 let r3 = self.radius[k].powi(3);
885 let r5 = self.radius[k].powi(5);
886 let mut block = Array2::<f64>::eye(dim);
887 let eps2 = self.epsilon * self.epsilon;
888 block.mapv_inplace(|v| -eps2 * v / r3);
889 for i in 0..dim {
890 for j in 0..dim {
891 block[[i, j]] += 3.0 * eps2 * row[i] * row[j] / r5;
892 }
893 }
894 out.push(block);
895 }
896 out
897 }
898
899 fn log_epsilon_beta_mixed_second_blocks(&self) -> Array2<f64> {
900 let mut out = self.signal_blocks.clone();
901 let eps2 = self.epsilon * self.epsilon;
902 for (k, mut row) in out.rows_mut().into_iter().enumerate() {
903 let norm2 = self.norm[k] * self.norm[k];
904 let scale = eps2 * (eps2 - 2.0 * norm2) / self.radius[k].powi(5);
905 row.mapv_inplace(|v| v * scale);
906 }
907 out
908 }
909
910 fn log_epsilon_betahessian_second_blocks(&self) -> Vec<Array2<f64>> {
911 let mut out = Vec::with_capacity(self.signal_blocks.nrows());
912 let eps2 = self.epsilon * self.epsilon;
913 for (k, row) in self.signal_blocks.rows().into_iter().enumerate() {
914 let dim = row.len();
915 let norm2 = self.norm[k] * self.norm[k];
916 let r5 = self.radius[k].powi(5);
917 let r7 = self.radius[k].powi(7);
918 let mut block = Array2::<f64>::eye(dim);
919 block.mapv_inplace(|v| eps2 * (eps2 - 2.0 * norm2) * v / r5);
920 for i in 0..dim {
921 for j in 0..dim {
922 block[[i, j]] += 3.0 * eps2 * (2.0 * norm2 - 3.0 * eps2) * row[i] * row[j] / r7;
923 }
924 }
925 out.push(block);
926 }
927 out
928 }
929
930 fn log_epsilon_betahessian_directional_blocks(
931 &self,
932 direction_blocks: &Array2<f64>,
933 ) -> Vec<Array2<f64>> {
934 let mut out = Vec::with_capacity(self.signal_blocks.nrows());
935 let eps2 = self.epsilon * self.epsilon;
936 for (k, (v, q)) in self
937 .signal_blocks
938 .rows()
939 .into_iter()
940 .zip(direction_blocks.rows().into_iter())
941 .enumerate()
942 {
943 let dim = v.len();
944 let dot = v.iter().zip(q.iter()).map(|(a, b)| a * b).sum::<f64>();
945 let r5 = self.radius[k].powi(5);
946 let r7 = self.radius[k].powi(7);
947 let mut block = Array2::<f64>::eye(dim);
948 block.mapv_inplace(|x| 3.0 * eps2 * dot * x / r5);
949 for i in 0..dim {
950 for j in 0..dim {
951 block[[i, j]] += 3.0 * eps2 * (q[i] * v[j] + v[i] * q[j]) / r5;
952 block[[i, j]] -= 15.0 * eps2 * dot * v[i] * v[j] / r7;
953 }
954 }
955 out.push(block);
956 }
957 out
958 }
959}
960
961fn scalar_operatorgradient(operator: &Array2<f64>, coeff: &Array1<f64>) -> Array1<f64> {
962 operator.t().dot(coeff)
963}
964
965fn scalar_operatorhessian(operator: &Array2<f64>, diag: &Array1<f64>) -> Array2<f64> {
966 let mut weighted = operator.clone();
967 for (k, &w) in diag.iter().enumerate() {
968 weighted.row_mut(k).mapv_inplace(|v| v * w);
969 }
970 let gram = operator.t().dot(&weighted);
971 (&gram + &gram.t().to_owned()) * 0.5
972}
973
974fn grouped_operatorgradient(
975 d1: &Array2<f64>,
976 dimension: usize,
977 blocks: &Array2<f64>,
978) -> Result<Array1<f64>, EstimationError> {
979 if blocks.ncols() != dimension {
980 crate::bail_invalid_estim!(
981 "grouped gradient block dimension mismatch: got {}, expected {dimension}",
982 blocks.ncols()
983 );
984 }
985 if d1.nrows() != blocks.nrows() * dimension {
986 crate::bail_invalid_estim!(
987 "grouped gradient row mismatch: D1 has {} rows, blocks imply {}",
988 d1.nrows(),
989 blocks.nrows() * dimension
990 );
991 }
992 let mut out = Array1::<f64>::zeros(d1.ncols());
993 for k in 0..blocks.nrows() {
994 let gk = d1
995 .slice(s![k * dimension..(k + 1) * dimension, ..])
996 .to_owned();
997 out += &gk.t().dot(&blocks.row(k));
998 }
999 Ok(out)
1000}
1001
1002fn grouped_operatorhessian(
1003 d1: &Array2<f64>,
1004 dimension: usize,
1005 blocks: &[Array2<f64>],
1006) -> Result<Array2<f64>, EstimationError> {
1007 if d1.nrows() != blocks.len() * dimension {
1008 crate::bail_invalid_estim!(
1009 "grouped Hessian row mismatch: D1 has {} rows, blocks imply {}",
1010 d1.nrows(),
1011 blocks.len() * dimension
1012 );
1013 }
1014 let p = d1.ncols();
1015 let mut out = Array2::<f64>::zeros((p, p));
1016 for (k, block) in blocks.iter().enumerate() {
1017 if block.nrows() != dimension || block.ncols() != dimension {
1018 crate::bail_invalid_estim!(
1019 "grouped Hessian block {k} has shape {}x{}, expected {}x{}",
1020 block.nrows(),
1021 block.ncols(),
1022 dimension,
1023 dimension
1024 );
1025 }
1026 let gk = d1
1027 .slice(s![k * dimension..(k + 1) * dimension, ..])
1028 .to_owned();
1029 out += &gk.t().dot(&block.dot(&gk));
1030 }
1031 Ok((&out + &out.t().to_owned()) * 0.5)
1032}
1033
1034#[derive(Clone)]
1035struct SpatialPenaltyExactState {
1036 magnitude: CharbonnierScalarBlockState,
1037 gradient: CharbonnierGroupedBlockState,
1038 curvature: CharbonnierGroupedBlockState,
1039}
1040
1041fn collocationgradient_blocks(
1042 gradrows: &Array1<f64>,
1043 dimension: usize,
1044) -> Result<Array2<f64>, EstimationError> {
1045 if dimension == 0 || !gradrows.len().is_multiple_of(dimension) {
1046 crate::bail_invalid_estim!(
1047 "invalid collocation gradient layout: rows={}, dimension={dimension}",
1048 gradrows.len()
1049 );
1050 }
1051 let p = gradrows.len() / dimension;
1052 let mut out = Array2::<f64>::zeros((p, dimension));
1053 for k in 0..p {
1054 for axis in 0..dimension {
1055 out[[k, axis]] = gradrows[k * dimension + axis];
1056 }
1057 }
1058 Ok(out)
1059}
1060
1061fn collocationhessian_blocks(
1062 hessianrows: &Array1<f64>,
1063 dimension: usize,
1064) -> Result<Array2<f64>, EstimationError> {
1065 let block_dim = dimension.checked_mul(dimension).ok_or_else(|| {
1066 EstimationError::InvalidInput("invalid collocation Hessian dimension overflow".to_string())
1067 })?;
1068 if block_dim == 0 || !hessianrows.len().is_multiple_of(block_dim) {
1069 crate::bail_invalid_estim!(
1070 "invalid collocation Hessian layout: rows={}, dimension={dimension}",
1071 hessianrows.len()
1072 );
1073 }
1074 let p = hessianrows.len() / block_dim;
1075 let mut out = Array2::<f64>::zeros((p, block_dim));
1076 for k in 0..p {
1077 for idx in 0..block_dim {
1078 out[[k, idx]] = hessianrows[k * block_dim + idx];
1079 }
1080 }
1081 Ok(out)
1082}
1083
1084impl SpatialPenaltyExactState {
1085 fn from_beta_local(
1086 beta_local: ArrayView1<'_, f64>,
1087 cache: &SpatialOperatorRuntimeCache,
1088 epsilons: [f64; 3],
1089 ) -> Result<Self, EstimationError> {
1090 let gradientrows = cache.d1.dot(&beta_local);
1120 let hessianrows = cache.d2.dot(&beta_local);
1121 Ok(Self {
1122 magnitude: CharbonnierScalarBlockState::from_signal(
1123 cache.d0.dot(&beta_local),
1124 epsilons[0],
1125 ),
1126 gradient: CharbonnierGroupedBlockState::from_signal_blocks(
1127 collocationgradient_blocks(&gradientrows, cache.dimension)?,
1128 epsilons[1],
1129 ),
1130 curvature: CharbonnierGroupedBlockState::from_signal_blocks(
1131 collocationhessian_blocks(&hessianrows, cache.dimension)?,
1132 epsilons[2],
1133 ),
1134 })
1135 }
1136
1137 fn absolute_collocation_magnitudes(&self) -> (Array1<f64>, Array1<f64>, Array1<f64>) {
1138 (
1139 self.magnitude.absolute_signal(),
1140 self.gradient.norm_signal(),
1141 self.curvature.norm_signal(),
1142 )
1143 }
1144}
1145
1146fn robust_epsilon_from_samples(values: &[f64], min_epsilon_cfg: f64) -> f64 {
1147 if values.is_empty() {
1148 return min_epsilon_cfg.max(1e-12);
1149 }
1150 let mut clean = values
1151 .iter()
1152 .copied()
1153 .filter(|v| v.is_finite() && *v >= 0.0)
1154 .collect::<Vec<_>>();
1155 if clean.is_empty() {
1156 return min_epsilon_cfg.max(1e-12);
1157 }
1158 clean.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
1159
1160 let n = clean.len();
1161 let median = quantile_from_sorted(&clean, 0.5);
1162 let q75 = quantile_from_sorted(&clean, 0.75);
1163 let q95 = quantile_from_sorted(&clean, 0.95);
1164
1165 let mut abs_dev = clean
1166 .iter()
1167 .map(|v| (v - median).abs())
1168 .filter(|v| v.is_finite())
1169 .collect::<Vec<_>>();
1170 abs_dev.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
1171 let mad = 1.4826 * quantile_from_sorted(&abs_dev, 0.5);
1172
1173 let mut scale = median.max(mad).max(q75);
1183
1184 let delta = (f64::EPSILON.sqrt() * q95.max(1.0))
1186 .max(min_epsilon_cfg)
1187 .max(1e-12);
1188 let s_min = min_epsilon_cfg.max(1e-12);
1189
1190 if scale <= delta {
1192 let rms = (clean.iter().map(|v| v * v).sum::<f64>() / n as f64).sqrt();
1193 scale = q95.max(rms);
1194 }
1195 if scale <= delta {
1196 scale = s_min;
1197 }
1198
1199 let kappa = 1.0_f64;
1202 (kappa * scale).max(s_min)
1203}
1204
1205fn extract_spatial_operator_runtime_caches(
1206 spec: &TermCollectionSpec,
1207 design: &TermCollectionDesign,
1208) -> Result<Vec<SpatialOperatorRuntimeCache>, EstimationError> {
1209 let smooth_start = design
1210 .design
1211 .ncols()
1212 .saturating_sub(design.smooth.total_smooth_cols());
1213 let mut out = Vec::<SpatialOperatorRuntimeCache>::new();
1214 for (term_idx, (termspec, term_fit)) in spec
1215 .smooth_terms
1216 .iter()
1217 .zip(design.smooth.terms.iter())
1218 .enumerate()
1219 {
1220 let Some(global_range) = design
1221 .smooth_term_penalty_range(term_idx)
1222 .map_err(EstimationError::InvalidInput)?
1223 else {
1224 continue;
1225 };
1226 let global_base_idx = global_range.start;
1227 let mut mass_local_idx = None;
1228 let mut tension_local_idx = None;
1229 let mut stiffness_local_idx = None;
1230 let mut mass_norm = None;
1231 let mut tension_norm = None;
1232 let mut stiffness_norm = None;
1233 for (active_local_idx, penalty) in term_fit.active_penalties.iter().enumerate() {
1234 match penalty.info.source {
1235 PenaltySource::OperatorMass => {
1236 mass_local_idx = Some(active_local_idx);
1237 mass_norm = Some(penalty.info.normalization_scale);
1238 }
1239 PenaltySource::OperatorTension => {
1240 tension_local_idx = Some(active_local_idx);
1241 tension_norm = Some(penalty.info.normalization_scale);
1242 }
1243 PenaltySource::OperatorStiffness => {
1244 stiffness_local_idx = Some(active_local_idx);
1245 stiffness_norm = Some(penalty.info.normalization_scale);
1246 }
1247 _ => {}
1248 }
1249 }
1250 let (
1263 Some(mass_local),
1264 Some(tension_local),
1265 Some(stiffness_local),
1266 Some(mass_scale),
1267 Some(tension_scale),
1268 Some(stiffness_scale),
1269 ) = (
1270 mass_local_idx,
1271 tension_local_idx,
1272 stiffness_local_idx,
1273 mass_norm,
1274 tension_norm,
1275 stiffness_norm,
1276 )
1277 else {
1278 continue;
1279 };
1280 let mass_global_idx = global_base_idx + mass_local;
1281 let tension_global_idx = global_base_idx + tension_local;
1282 let stiffness_global_idx = global_base_idx + stiffness_local;
1283
1284 let (feature_cols, mut d0, mut d1, mut d2, collocation_points, dim, center_mass_rows) =
1285 match (&termspec.basis, &term_fit.metadata) {
1286 (
1287 SmoothBasisSpec::Matern { feature_cols, .. },
1288 BasisMetadata::Matern {
1289 centers,
1290 length_scale,
1291 nu,
1292 include_intercept,
1293 identifiability_transform,
1294 aniso_log_scales,
1295 input_scale,
1296 ..
1297 },
1298 ) => {
1299 let collocation_length_scale =
1305 input_scale.to_standardized_units(*length_scale);
1306 let ops = build_matern_collocation_operator_matrices(
1307 centers.view(),
1308 None,
1309 collocation_length_scale,
1310 *nu,
1311 *include_intercept,
1312 identifiability_transform.as_ref().map(|z| z.view()),
1313 aniso_log_scales.as_deref(),
1314 )?;
1315 (
1316 feature_cols.clone(),
1317 ops.d0,
1318 ops.d1,
1319 ops.d2,
1320 ops.collocation_points,
1321 centers.ncols(),
1322 false,
1323 )
1324 }
1325 (
1326 SmoothBasisSpec::Duchon { feature_cols, .. },
1327 BasisMetadata::Duchon {
1328 centers,
1329 length_scale,
1330 power,
1331 nullspace_order,
1332 identifiability_transform,
1333 input_scale,
1334 aniso_log_scales,
1335 operator_collocation_points: Some(collocation_points),
1336 radial_reparam,
1337 ..
1338 },
1339 ) => {
1340 let collocation_length_scale = (*length_scale)
1341 .map(|length| input_scale.to_standardized_units(length));
1342 let ops =
1343 gam_terms::basis::build_duchon_collocation_operator_matriceswithworkspace(
1344 centers.view(),
1345 collocation_points.view(),
1346 None,
1347 collocation_length_scale,
1348 *power,
1349 *nullspace_order,
1350 aniso_log_scales.as_deref(),
1351 identifiability_transform.as_ref().map(|z| z.view()),
1352 2,
1353 radial_reparam.as_ref().map(|v| v.view()),
1354 &mut BasisWorkspace::default(),
1355 )?;
1356 (
1357 feature_cols.clone(),
1358 ops.d0,
1359 ops.d1,
1360 ops.d2,
1361 ops.collocation_points,
1362 centers.ncols(),
1363 true,
1364 )
1365 }
1366 _ => continue,
1367 };
1368 if center_mass_rows && d0.nrows() > 0 && d0.ncols() > 0 {
1369 let means = d0.sum_axis(Axis(0)).mapv(|v| v / d0.nrows() as f64);
1370 for mut row in d0.rows_mut() {
1371 row -= &means;
1372 }
1373 }
1374
1375 let mass_scale = mass_scale.max(1e-12).sqrt();
1393 let tension_scale = tension_scale.max(1e-12).sqrt();
1394 let stiffness_scale = stiffness_scale.max(1e-12).sqrt();
1395 d0.mapv_inplace(|v| v / mass_scale);
1396 d1.mapv_inplace(|v| v / tension_scale);
1397 d2.mapv_inplace(|v| v / stiffness_scale);
1398
1399 let coeff_global_range =
1400 (smooth_start + term_fit.coeff_range.start)..(smooth_start + term_fit.coeff_range.end);
1401 if d0.ncols() != coeff_global_range.len()
1402 || d1.ncols() != coeff_global_range.len()
1403 || d2.ncols() != coeff_global_range.len()
1404 {
1405 crate::bail_invalid_estim!(
1406 "spatial operator dimension mismatch for term '{}': D0 cols={}, D1 cols={}, D2 cols={}, coeffs={}",
1407 term_fit.name,
1408 d0.ncols(),
1409 d1.ncols(),
1410 d2.ncols(),
1411 coeff_global_range.len()
1412 );
1413 }
1414 out.push(SpatialOperatorRuntimeCache {
1415 termname: term_fit.name.clone(),
1416 feature_cols,
1417 coeff_global_range,
1418 mass_penalty_global_idx: mass_global_idx,
1419 tension_penalty_global_idx: tension_global_idx,
1420 stiffness_penalty_global_idx: stiffness_global_idx,
1421 d0,
1422 d1,
1423 d2,
1424 collocation_points,
1425 dimension: dim,
1426 });
1427 }
1428 Ok(out)
1429}
1430
1431fn scalar_operator_response_variance(
1443 operator: &Array2<f64>,
1444 cov_local: &Array2<f64>,
1445) -> Array1<f64> {
1446 Array1::from_iter(operator.rows().into_iter().map(|row| {
1447 let s = cov_local.dot(&row);
1448 row.dot(&s).max(0.0)
1449 }))
1450}
1451
1452fn grouped_operator_response_variance(
1463 operator: &Array2<f64>,
1464 block_dim: usize,
1465 cov_local: &Array2<f64>,
1466) -> Result<Array1<f64>, EstimationError> {
1467 if block_dim == 0 || !operator.nrows().is_multiple_of(block_dim) {
1468 crate::bail_invalid_estim!(
1469 "grouped variance row layout invalid: rows={}, block_dim={block_dim}",
1470 operator.nrows()
1471 );
1472 }
1473 let p = operator.nrows() / block_dim;
1474 let mut out = Array1::<f64>::zeros(p);
1475 for k in 0..p {
1476 let mut acc = 0.0;
1477 for axis in 0..block_dim {
1478 let row = operator.row(k * block_dim + axis);
1479 let s = cov_local.dot(&row);
1480 acc += row.dot(&s);
1481 }
1482 out[k] = acc.max(0.0);
1483 }
1484 Ok(out)
1485}
1486
1487fn compute_spatial_adaptiveweights_for_beta(
1488 beta: &Array1<f64>,
1489 caches: &[SpatialOperatorRuntimeCache],
1490 epsilon_0: f64,
1491 epsilon_g: f64,
1492 epsilon_c: f64,
1493 weight_floor: f64,
1494 weight_ceiling: f64,
1495 beta_covariance: Option<&Array2<f64>>,
1496) -> Result<Vec<SpatialAdaptiveWeights>, EstimationError> {
1497 caches
1529 .iter()
1530 .map(|cache| {
1531 let beta_local = beta.slice(s![cache.coeff_global_range.clone()]);
1532 let exact = SpatialPenaltyExactState::from_beta_local(
1533 beta_local,
1534 cache,
1535 [epsilon_0, epsilon_g, epsilon_c],
1536 )?;
1537 let cov_local = beta_covariance.map(|cov| {
1538 cov.slice(s![
1539 cache.coeff_global_range.clone(),
1540 cache.coeff_global_range.clone()
1541 ])
1542 .to_owned()
1543 });
1544 let dim = cache.dimension;
1545 let (var_0, var_g, var_c) = match cov_local.as_ref() {
1546 Some(cov) => (
1547 scalar_operator_response_variance(&cache.d0, cov),
1548 grouped_operator_response_variance(&cache.d1, dim, cov)?,
1549 grouped_operator_response_variance(&cache.d2, dim * dim, cov)?,
1550 ),
1551 None => (
1552 Array1::<f64>::zeros(exact.magnitude.signal.len()),
1553 Array1::<f64>::zeros(exact.gradient.norm.len()),
1554 Array1::<f64>::zeros(exact.curvature.norm.len()),
1555 ),
1556 };
1557 let (_, inv_0) = exact.magnitude.surrogateweights_posterior_snr(
1558 &var_0,
1559 weight_floor,
1560 weight_ceiling,
1561 );
1562 let (_, inv_g) =
1563 exact
1564 .gradient
1565 .surrogateweights_posterior_snr(&var_g, weight_floor, weight_ceiling);
1566 let (_, inv_c) = exact.curvature.surrogateweights_posterior_snr(
1567 &var_c,
1568 weight_floor,
1569 weight_ceiling,
1570 );
1571 Ok(SpatialAdaptiveWeights {
1572 inv_magweight: inv_0,
1573 invgradweight: inv_g,
1574 inv_lapweight: inv_c,
1575 })
1576 })
1577 .collect()
1578}
1579
1580fn compute_initial_epsilons(
1581 beta: &Array1<f64>,
1582 caches: &[SpatialOperatorRuntimeCache],
1583 min_epsilon: f64,
1584) -> Result<(f64, f64, f64), EstimationError> {
1585 let mut fvals = Vec::<f64>::new();
1586 let mut gvals = Vec::<f64>::new();
1587 let mut cvals = Vec::<f64>::new();
1588 for cache in caches {
1589 let beta_local = beta.slice(s![cache.coeff_global_range.clone()]);
1590 let exact = SpatialPenaltyExactState::from_beta_local(
1591 beta_local,
1592 cache,
1593 [min_epsilon, min_epsilon, min_epsilon],
1594 )?;
1595 let (f, g, c) = exact.absolute_collocation_magnitudes();
1596 fvals.extend(f.iter().copied());
1597 gvals.extend(g.iter().copied());
1598 cvals.extend(c.iter().copied());
1599 }
1600 let eps_0 = robust_epsilon_from_samples(&fvals, min_epsilon);
1606 let eps_g = robust_epsilon_from_samples(&gvals, min_epsilon);
1607 let eps_c = robust_epsilon_from_samples(&cvals, min_epsilon);
1608 Ok((eps_0, eps_g, eps_c))
1609}
1610
1611fn exact_spatial_adaptive_penalty_index_set(
1612 caches: &[SpatialOperatorRuntimeCache],
1613) -> BTreeSet<usize> {
1614 let mut out = BTreeSet::new();
1615 for cache in caches {
1616 out.insert(cache.mass_penalty_global_idx);
1617 out.insert(cache.tension_penalty_global_idx);
1618 out.insert(cache.stiffness_penalty_global_idx);
1619 }
1620 out
1621}
1622
1623fn checked_fit_log_lambdas(
1624 lambdas: &Array1<f64>,
1625 context: &str,
1626) -> Result<Array1<f64>, EstimationError> {
1627 let values = lambdas
1628 .iter()
1629 .copied()
1630 .enumerate()
1631 .map(|(coordinate, lambda)| {
1632 gam_problem::checked_log_strength(lambda).map_err(|error| {
1633 EstimationError::InvalidInput(format!(
1634 "{context} lambda coordinate {coordinate} is outside the canonical physical-strength domain: {error}"
1635 ))
1636 })
1637 })
1638 .collect::<Result<Vec<_>, _>>()?;
1639 Ok(Array1::from_vec(values))
1640}
1641
1642fn build_spatial_adaptive_hyperspecs(cache_count: usize) -> Vec<SpatialAdaptiveHyperSpec> {
1643 let mut out = Vec::with_capacity(cache_count * 3 + 3);
1644 for cache_index in 0..cache_count {
1645 out.push(SpatialAdaptiveHyperSpec {
1646 cache_index,
1647 kind: SpatialAdaptiveHyperKind::LogLambdaMagnitude,
1648 });
1649 out.push(SpatialAdaptiveHyperSpec {
1650 cache_index,
1651 kind: SpatialAdaptiveHyperKind::LogLambdaGradient,
1652 });
1653 out.push(SpatialAdaptiveHyperSpec {
1654 cache_index,
1655 kind: SpatialAdaptiveHyperKind::LogLambdaCurvature,
1656 });
1657 }
1658 out.push(SpatialAdaptiveHyperSpec {
1659 cache_index: 0,
1660 kind: SpatialAdaptiveHyperKind::LogEpsilonMagnitude,
1661 });
1662 out.push(SpatialAdaptiveHyperSpec {
1663 cache_index: 0,
1664 kind: SpatialAdaptiveHyperKind::LogEpsilonGradient,
1665 });
1666 out.push(SpatialAdaptiveHyperSpec {
1667 cache_index: 0,
1668 kind: SpatialAdaptiveHyperKind::LogEpsilonCurvature,
1669 });
1670 out
1671}
1672
1673fn penalty_matrixwith_local_block(
1674 total_dim: usize,
1675 coeff_range: Range<usize>,
1676 local: &Array2<f64>,
1677) -> Array2<f64> {
1678 let mut out = Array2::<f64>::zeros((total_dim, total_dim));
1679 out.slice_mut(s![coeff_range.clone(), coeff_range])
1680 .assign(local);
1681 out
1682}
1683
1684fn fit_term_collectionwith_exact_spatial_adaptive_regularization(
1685 baseline: FittedTermCollection,
1686 y: ArrayView1<'_, f64>,
1687 weights: ArrayView1<'_, f64>,
1688 offset: ArrayView1<'_, f64>,
1689 family: LikelihoodSpec,
1690 options: &FitOptions,
1691 runtime_caches: &[SpatialOperatorRuntimeCache],
1692) -> Result<FittedTermCollection, EstimationError> {
1693 let adaptive_opts = options.adaptive_regularization.clone().unwrap_or_default();
1722 let adaptive_penalty_indices = exact_spatial_adaptive_penalty_index_set(runtime_caches);
1723 let p_total = baseline.design.design.ncols();
1724 if baseline.fit.lambdas.len() != baseline.design.penalties.len() {
1725 crate::bail_invalid_estim!(
1726 "exact spatial adaptive fit received {} baseline lambdas for {} penalties",
1727 baseline.fit.lambdas.len(),
1728 baseline.design.penalties.len(),
1729 );
1730 }
1731 let baseline_log_lambdas =
1732 checked_fit_log_lambdas(&baseline.fit.lambdas, "exact spatial adaptive baseline")?;
1733 for (cache_idx, cache) in runtime_caches.iter().enumerate() {
1734 for (operator, penalty_idx) in [
1735 ("mass", cache.mass_penalty_global_idx),
1736 ("tension", cache.tension_penalty_global_idx),
1737 ("stiffness", cache.stiffness_penalty_global_idx),
1738 ] {
1739 if penalty_idx >= baseline.fit.lambdas.len() {
1740 crate::bail_invalid_estim!(
1741 "exact spatial adaptive cache {cache_idx} {operator} penalty index {penalty_idx} is out of bounds for {} baseline lambdas",
1742 baseline.fit.lambdas.len(),
1743 );
1744 }
1745 }
1746 }
1747 struct RetainedPenaltySetup {
1748 global_idx: usize,
1749 global_penalty: Array2<f64>,
1750 nullspace_dim: usize,
1751 log_lambda: f64,
1752 }
1753 use rayon::iter::{IndexedParallelIterator, IntoParallelRefIterator, ParallelIterator};
1754 let retained_setups = baseline
1755 .design
1756 .penalties
1757 .par_iter()
1758 .enumerate()
1759 .map(|(idx, bp)| {
1760 if adaptive_penalty_indices.contains(&idx) {
1761 return None;
1762 }
1763 Some(RetainedPenaltySetup {
1764 global_idx: idx,
1765 global_penalty: bp.to_global(p_total),
1766 nullspace_dim: baseline
1767 .design
1768 .nullspace_dims
1769 .get(idx)
1770 .copied()
1771 .unwrap_or(0),
1772 log_lambda: baseline_log_lambdas[idx],
1773 })
1774 })
1775 .collect::<Vec<_>>();
1776 let retained_count = retained_setups
1777 .iter()
1778 .filter(|setup| setup.is_some())
1779 .count();
1780 let mut retained_penalties = Vec::<Array2<f64>>::with_capacity(retained_count);
1781 let mut retained_nullspace_dims = Vec::<usize>::with_capacity(retained_count);
1782 let mut retained_log_lambdas = Vec::<f64>::with_capacity(retained_count);
1783 let mut retained_global_indices = Vec::<usize>::with_capacity(retained_count);
1784 for setup in retained_setups.into_iter().flatten() {
1785 retained_penalties.push(setup.global_penalty);
1786 retained_nullspace_dims.push(setup.nullspace_dim);
1787 retained_log_lambdas.push(setup.log_lambda);
1788 retained_global_indices.push(setup.global_idx);
1789 }
1790
1791 let (eps_0_init, eps_g_init, eps_c_init) = compute_initial_epsilons(
1792 &baseline.fit.beta,
1793 runtime_caches,
1794 adaptive_opts.min_epsilon,
1795 )?;
1796 let mut initial_theta =
1797 Array1::<f64>::zeros(retained_penalties.len() + runtime_caches.len() * 3 + 3);
1798 for (idx, value) in retained_log_lambdas.iter().enumerate() {
1799 initial_theta[idx] = *value;
1800 }
1801 let adaptive_log_lambda_components = runtime_caches
1802 .par_iter()
1803 .map(|cache| {
1804 [
1805 baseline_log_lambdas[cache.mass_penalty_global_idx],
1806 baseline_log_lambdas[cache.tension_penalty_global_idx],
1807 baseline_log_lambdas[cache.stiffness_penalty_global_idx],
1808 ]
1809 })
1810 .collect::<Vec<_>>();
1811 let mut at = retained_penalties.len();
1812 for logs in &adaptive_log_lambda_components {
1813 initial_theta[at] = logs[0];
1814 initial_theta[at + 1] = logs[1];
1815 initial_theta[at + 2] = logs[2];
1816 at += 3;
1817 }
1818 let minimum_log_epsilon = gam_problem::checked_log_strength(adaptive_opts.min_epsilon)
1819 .map_err(|error| {
1820 EstimationError::InvalidInput(format!(
1821 "adaptive minimum epsilon is outside the canonical positive-strength domain: {error}"
1822 ))
1823 })?;
1824 for (slot, epsilon) in [eps_0_init, eps_g_init, eps_c_init].into_iter().enumerate() {
1825 initial_theta[at + slot] =
1826 gam_problem::checked_log_strength(epsilon.max(adaptive_opts.min_epsilon)).map_err(
1827 |error| {
1828 EstimationError::InvalidInput(format!(
1829 "adaptive initial epsilon coordinate {slot} is outside the canonical positive-strength domain: {error}"
1830 ))
1831 },
1832 )?;
1833 }
1834
1835 let hyperspecs = build_spatial_adaptive_hyperspecs(runtime_caches.len());
1836 let zero_psi_op: std::sync::Arc<dyn gam_custom_family::CustomFamilyPsiDerivativeOperator> =
1837 std::sync::Arc::new(gam_custom_family::ZeroPsiDerivativeOperator::new(
1838 baseline.design.design.nrows(),
1839 baseline.design.design.ncols(),
1840 ));
1841 let derivative_blocks = vec![
1842 hyperspecs
1843 .par_iter()
1844 .map(|_| CustomFamilyBlockPsiDerivative {
1845 penalty_index: None,
1846 x_psi: Array2::<f64>::zeros((0, 0)),
1847 s_psi: Array2::<f64>::zeros((0, 0)),
1848 s_psi_components: None,
1849 s_psi_penalty_components: None,
1850 x_psi_psi: None,
1851 s_psi_psi: None,
1852 s_psi_psi_components: None,
1853 s_psi_psi_penalty_components: None,
1854 implicit_operator: Some(std::sync::Arc::clone(&zero_psi_op)),
1855 implicit_axis: 0,
1856 implicit_group_id: None,
1857 })
1858 .collect::<Vec<_>>(),
1859 ];
1860
1861 let mixture_link_state = options
1862 .mixture_link
1863 .clone()
1864 .as_ref()
1865 .map(state_fromspec)
1866 .transpose()
1867 .map_err(EstimationError::InvalidInput)?;
1868 let sas_link_state = options
1869 .sas_link
1870 .map(|spec| {
1871 if family.is_binomial_beta_logistic() {
1872 state_from_beta_logisticspec(spec)
1873 } else {
1874 state_from_sasspec(spec)
1875 }
1876 })
1877 .transpose()
1878 .map_err(EstimationError::InvalidInput)?;
1879 let latent_cloglog_state = options.latent_cloglog;
1880 let shared_y = Arc::new(y.to_owned());
1881 let sharedweights = Arc::new(weights.to_owned());
1882 let shared_design = baseline
1883 .design
1884 .design
1885 .try_to_dense_arc("spatial adaptive exact hyperfit design")
1886 .map_err(EstimationError::InvalidInput)?;
1887 let shared_offset = Arc::new(offset.to_owned());
1888 let shared_runtime_caches = Arc::new(runtime_caches.to_vec());
1889 let shared_hyperspecs = Arc::new(hyperspecs.clone());
1890 let zero_quadratic = ValidatedFixedQuadraticHessian::zero(
1891 baseline.design.design.ncols(),
1892 )
1893 .map_err(EstimationError::InvalidInput)?;
1894 let base_family = SpatialAdaptiveExactFamily {
1895 family: family.clone(),
1896 latent_cloglog_state,
1897 mixture_link_state: mixture_link_state.clone(),
1898 sas_link_state,
1899 y: shared_y.clone(),
1900 weights: sharedweights.clone(),
1901 design: shared_design.clone(),
1902 offset: shared_offset.clone(),
1903 linear_constraints: baseline.design.linear_constraints.clone(),
1904 runtime_caches: shared_runtime_caches.clone(),
1905 adaptive_params: Vec::new(),
1906 fixed_quadratic_hessian: zero_quadratic.clone(),
1907 hyperspecs: shared_hyperspecs.clone(),
1908 exact_eval_cache: Arc::new(Mutex::new(None)),
1909 };
1910
1911 let rho_dim = retained_penalties.len();
1912 let operator_slots_end = rho_dim + runtime_caches.len() * 3;
1913 const UNIFIED_LOG_WINDOW: f64 = 6.0;
1923 const RETAINED_LAMBDA_LOG_LOWER_FLOOR: f64 = -30.0;
1924 const RETAINED_LAMBDA_LOG_UPPER_CAP: f64 = 30.0;
1925 const OPERATOR_LAMBDA_LOG_LOWER_FLOOR: f64 = -10.0;
1926 const OPERATOR_LAMBDA_LOG_UPPER_CAP: f64 = 30.0;
1927 let epsilon_floor_log = minimum_log_epsilon;
1928 let anchored_bound = |idx: usize, sign: f64| -> f64 {
1929 let raw = initial_theta[idx] + sign * UNIFIED_LOG_WINDOW;
1930 if idx < rho_dim {
1931 raw.clamp(
1932 RETAINED_LAMBDA_LOG_LOWER_FLOOR,
1933 RETAINED_LAMBDA_LOG_UPPER_CAP,
1934 )
1935 } else if idx < operator_slots_end {
1936 raw.clamp(
1937 OPERATOR_LAMBDA_LOG_LOWER_FLOOR,
1938 OPERATOR_LAMBDA_LOG_UPPER_CAP,
1939 )
1940 } else {
1941 raw.clamp(epsilon_floor_log, gam_problem::LOG_STRENGTH_MAX)
1942 }
1943 };
1944 let eps_lower =
1945 Array1::from_iter((0..initial_theta.len()).map(|idx| anchored_bound(idx, -1.0)));
1946 let eps_upper = Array1::from_iter((0..initial_theta.len()).map(|idx| anchored_bound(idx, 1.0)));
1947 let blockspec = ParameterBlockSpec {
1948 name: "eta".to_string(),
1949 design: baseline.design.design.clone(),
1950 offset: offset.to_owned(),
1951 penalties: retained_penalties
1952 .iter()
1953 .cloned()
1954 .map(PenaltyMatrix::Dense)
1955 .collect(),
1956 nullspace_dims: retained_nullspace_dims.clone(),
1957 initial_log_lambdas: Array1::from_vec(retained_log_lambdas.clone()),
1958 initial_beta: Some(baseline.fit.beta.clone()),
1959 gauge_priority: 100,
1960 jacobian_callback: None,
1961 stacked_design: None,
1962 stacked_offset: None,
1963 };
1964 let screening_cap = Arc::new(AtomicUsize::new(0));
1965 let outer_opts = BlockwiseFitOptions {
1966 inner_max_cycles: options.max_iter,
1967 inner_tol: options.tol,
1968 outer_max_iter: options.max_iter,
1969 outer_tol: options.tol,
1970 compute_covariance: false,
1971 screening_max_inner_iterations: Some(Arc::clone(&screening_cap)),
1972 ..BlockwiseFitOptions::default()
1973 };
1974
1975 use gam_problem::{DeclaredHessianForm, Derivative, HessianValue, OuterEval};
1976 use gam_solve::rho_optimizer::OuterProblem;
1977
1978 struct SpatialAdaptiveOuterState {
1979 warm_cache: Option<CustomFamilyWarmStart>,
1980 terminal_mode: Option<(Array1<f64>, f64, CustomFamilyOwnedMode)>,
1981 last_eval: Option<(
1982 Array1<f64>,
1983 f64,
1984 Array1<f64>,
1985 HessianValue,
1986 CustomFamilyWarmStart,
1987 )>,
1988 }
1989
1990 struct DecodedSpatialAdaptiveTheta {
1991 rho: Array1<f64>,
1992 retained_lambdas: Array1<f64>,
1993 adaptive_params: Vec<SpatialAdaptiveTermHyperParams>,
1994 epsilon: [f64; 3],
1995 }
1996
1997 let n_theta = initial_theta.len();
1998
1999 let theta_bounds = Some((eps_lower.clone(), eps_upper.clone()));
2000
2001 let decode_theta =
2002 |theta: &Array1<f64>| -> Result<DecodedSpatialAdaptiveTheta, EstimationError> {
2003 let physical = gam_problem::checked_exp_log_strengths(theta.iter().copied()).map_err(
2004 |error| {
2005 EstimationError::InvalidInput(format!(
2006 "exact spatial adaptive outer coordinate is outside the canonical log-strength domain: {error}"
2007 ))
2008 },
2009 )?;
2010 let rho = theta.slice(s![..rho_dim]).to_owned();
2011 let retained_lambdas = Array1::from_vec(physical[..rho_dim].to_vec());
2012 let adaptive_lambda_start = rho_dim;
2013 let adaptive_lambda_end = adaptive_lambda_start + runtime_caches.len() * 3;
2014 let eps = [
2015 physical[adaptive_lambda_end],
2016 physical[adaptive_lambda_end + 1],
2017 physical[adaptive_lambda_end + 2],
2018 ];
2019 let adaptive_params = runtime_caches
2020 .iter()
2021 .enumerate()
2022 .map(|(cache_idx, _)| SpatialAdaptiveTermHyperParams {
2023 lambda: [
2024 physical[adaptive_lambda_start + cache_idx * 3],
2025 physical[adaptive_lambda_start + cache_idx * 3 + 1],
2026 physical[adaptive_lambda_start + cache_idx * 3 + 2],
2027 ],
2028 epsilon: eps,
2029 })
2030 .collect::<Vec<_>>();
2031 Ok(DecodedSpatialAdaptiveTheta {
2032 rho,
2033 retained_lambdas,
2034 adaptive_params,
2035 epsilon: eps,
2036 })
2037 };
2038 let clamp_theta = |theta: &Array1<f64>| -> Array1<f64> {
2044 Array1::from_shape_fn(theta.len(), |i| theta[i].clamp(eps_lower[i], eps_upper[i]))
2045 };
2046 let realize_hyper_layout = |theta: &Array1<f64>| {
2047 gam_custom_family::CustomFamilyHyperLayout::new(
2048 derivative_blocks.clone(),
2049 Vec::new(),
2050 theta.slice(s![rho_dim..]).to_owned(),
2051 )
2052 .map_err(EstimationError::InvalidInput)
2053 };
2054 let analytic_outer_hessian_available =
2055 gam_custom_family::joint_exact_analytic_outer_hessian_available()
2056 && base_family
2057 .exact_outer_derivative_order(std::slice::from_ref(&blockspec), &outer_opts)
2058 .has_hessian()
2059 && gam_custom_family::exact_newton_outer_geometry_supports_second_order_solver(
2060 &base_family,
2061 );
2062 let problem = OuterProblem::new(n_theta)
2068 .with_gradient(Derivative::Analytic)
2069 .with_hessian(if analytic_outer_hessian_available {
2070 DeclaredHessianForm::Either
2071 } else {
2072 DeclaredHessianForm::Unavailable
2073 })
2074 .with_prefer_gradient_only(true)
2075 .with_fallback_policy(gam_solve::rho_optimizer::FallbackPolicy::Disabled)
2076 .with_psi_dim(n_theta.saturating_sub(rho_dim))
2077 .with_tolerance(options.tol)
2078 .with_rel_cost_tolerance(Some(options.tol))
2083 .with_max_iter(options.max_iter)
2084 .with_seed_config(gam_problem::SeedConfig::default())
2085 .with_screening_cap(Arc::clone(&screening_cap))
2086 .with_initial_rho(initial_theta.clone());
2087 let problem = if let Some((lo, hi)) = theta_bounds {
2088 problem.with_bounds(lo, hi)
2089 } else {
2090 problem
2091 };
2092
2093 let eval_outer = |st: &mut SpatialAdaptiveOuterState,
2094 theta: &Array1<f64>,
2095 order: gam_solve::rho_optimizer::OuterEvalOrder|
2096 -> Result<OuterEval, EstimationError> {
2097 let decoded = decode_theta(theta)?;
2098
2099 if let Some((cached_theta, cached_cost, cached_grad, cached_hess, cached_warm)) =
2100 &st.last_eval
2101 && cached_theta.len() == theta.len()
2102 && cached_theta
2103 .iter()
2104 .zip(theta.iter())
2105 .all(|(&a, &b)| a.to_bits() == b.to_bits())
2106 && st
2107 .terminal_mode
2108 .as_ref()
2109 .is_some_and(|(mode_theta, mode_objective, _)| {
2110 mode_theta.len() == theta.len()
2111 && mode_theta
2112 .iter()
2113 .zip(theta.iter())
2114 .all(|(&a, &b)| a.to_bits() == b.to_bits())
2115 && mode_objective.to_bits() == cached_cost.to_bits()
2116 })
2117 && (!matches!(
2118 order,
2119 gam_solve::rho_optimizer::OuterEvalOrder::ValueGradientHessian
2120 ) || analytic_outer_hessian_available)
2121 {
2122 st.warm_cache = Some(cached_warm.clone());
2123 return Ok(OuterEval {
2124 cost: *cached_cost,
2125 gradient: cached_grad.clone(),
2126 hessian: if matches!(
2127 order,
2128 gam_solve::rho_optimizer::OuterEvalOrder::ValueGradientHessian
2129 ) && analytic_outer_hessian_available
2130 {
2131 cached_hess.clone()
2132 } else {
2133 HessianValue::Unavailable
2134 },
2135 inner_beta_hint: None,
2136 });
2137 }
2138
2139 let family_eval =
2140 base_family.with_adaptive_params(decoded.adaptive_params, zero_quadratic.clone());
2141 let hyper_layout = realize_hyper_layout(theta)?;
2142 let need_hessian = matches!(
2143 order,
2144 gam_solve::rho_optimizer::OuterEvalOrder::ValueGradientHessian
2145 ) && analytic_outer_hessian_available;
2146 let owned = evaluate_custom_family_joint_hyper_owned(
2147 &family_eval,
2148 std::slice::from_ref(&blockspec),
2149 &outer_opts,
2150 &decoded.rho,
2151 &hyper_layout,
2152 st.warm_cache.as_ref(),
2153 if need_hessian {
2154 gam_solve::estimate::reml::reml_outer_engine::EvalMode::ValueGradientHessian
2155 } else {
2156 gam_solve::estimate::reml::reml_outer_engine::EvalMode::ValueAndGradient
2157 },
2158 )
2159 .map_err(|e| {
2160 EstimationError::RemlOptimizationFailed(format!("spatial adaptive eval failed: {e}"))
2161 })?;
2162 if !owned.result.inner_converged {
2163 st.warm_cache = Some(owned.result.warm_start.clone());
2164 return Err(EstimationError::RemlOptimizationFailed(
2165 "exact spatial adaptive inner solve did not converge".to_string(),
2166 ));
2167 }
2168 if !owned.result.objective.is_finite()
2169 || owned.result.gradient.iter().any(|v| !v.is_finite())
2170 {
2171 return Err(EstimationError::RemlOptimizationFailed(
2172 "exact spatial adaptive objective returned non-finite values".to_string(),
2173 ));
2174 }
2175 let hessian_result = if need_hessian {
2176 if !owned.result.outer_hessian.is_analytic() {
2177 return Err(EstimationError::RemlOptimizationFailed(
2178 "exact spatial adaptive objective did not return an exact outer Hessian"
2179 .to_string(),
2180 ));
2181 }
2182 match owned.result.outer_hessian.dim() {
2183 Some(dim) if dim == theta.len() => {}
2184 Some(dim) => {
2185 return Err(EstimationError::RemlOptimizationFailed(format!(
2186 "exact spatial adaptive outer Hessian dimension mismatch: got {dim}, expected {}",
2187 theta.len(),
2188 )));
2189 }
2190 None => {
2191 return Err(EstimationError::RemlOptimizationFailed(
2192 "exact spatial adaptive objective did not report an outer Hessian dimension"
2193 .to_string(),
2194 ));
2195 }
2196 }
2197 st.last_eval = Some((
2198 theta.to_owned(),
2199 owned.result.objective,
2200 owned.result.gradient.clone(),
2201 owned.result.outer_hessian.clone(),
2202 owned.result.warm_start.clone(),
2203 ));
2204 owned.result.outer_hessian
2205 } else {
2206 HessianValue::Unavailable
2207 };
2208 let objective = owned.result.objective;
2209 let gradient = owned.result.gradient;
2210 st.warm_cache = Some(owned.result.warm_start);
2211 st.terminal_mode = Some((theta.to_owned(), objective, owned.mode));
2212 Ok(OuterEval {
2213 cost: objective,
2214 gradient,
2215 hessian: hessian_result,
2216 inner_beta_hint: None,
2217 })
2218 };
2219
2220 let mut obj = problem.build_objective_with_screening_proxy(
2221 SpatialAdaptiveOuterState {
2222 warm_cache: None,
2223 terminal_mode: None,
2224 last_eval: None,
2225 },
2226 |st: &mut SpatialAdaptiveOuterState, theta: &Array1<f64>| {
2227 let theta = clamp_theta(theta);
2228 let DecodedSpatialAdaptiveTheta {
2229 rho,
2230 adaptive_params,
2231 ..
2232 } = decode_theta(&theta)?;
2233 let family_eval =
2234 base_family.with_adaptive_params(adaptive_params, zero_quadratic.clone());
2235 let hyper_layout = realize_hyper_layout(&theta)?;
2236 let owned = evaluate_custom_family_joint_hyper_owned(
2237 &family_eval,
2238 std::slice::from_ref(&blockspec),
2239 &outer_opts,
2240 &rho,
2241 &hyper_layout,
2242 st.warm_cache.as_ref(),
2243 gam_solve::estimate::reml::reml_outer_engine::EvalMode::ValueOnly,
2244 )
2245 .map_err(|e| {
2246 EstimationError::RemlOptimizationFailed(format!(
2247 "spatial adaptive cost eval failed: {e}"
2248 ))
2249 })?;
2250 if !owned.result.inner_converged {
2251 st.warm_cache = Some(owned.result.warm_start);
2252 return Err(EstimationError::RemlOptimizationFailed(
2253 "exact spatial adaptive cost inner solve did not converge".to_string(),
2254 ));
2255 }
2256 let objective = owned.result.objective;
2257 st.warm_cache = Some(owned.result.warm_start);
2258 st.terminal_mode = Some((theta, objective, owned.mode));
2259 Ok(objective)
2260 },
2261 |st: &mut SpatialAdaptiveOuterState, theta: &Array1<f64>| {
2262 eval_outer(
2263 st,
2264 theta,
2265 if analytic_outer_hessian_available {
2266 gam_solve::rho_optimizer::OuterEvalOrder::ValueGradientHessian
2267 } else {
2268 gam_solve::rho_optimizer::OuterEvalOrder::ValueAndGradient
2269 },
2270 )
2271 },
2272 |st: &mut SpatialAdaptiveOuterState,
2273 theta: &Array1<f64>,
2274 order: gam_solve::rho_optimizer::OuterEvalOrder| { eval_outer(st, theta, order) },
2275 Some(|st: &mut SpatialAdaptiveOuterState| {
2276 st.warm_cache = None;
2277 st.terminal_mode = None;
2278 st.last_eval = None;
2279 }),
2280 Some(|st: &mut SpatialAdaptiveOuterState, theta: &Array1<f64>| {
2281 let theta = clamp_theta(theta);
2282 let DecodedSpatialAdaptiveTheta {
2283 rho,
2284 adaptive_params,
2285 ..
2286 } = decode_theta(&theta)?;
2287 let family_eval =
2288 base_family.with_adaptive_params(adaptive_params, zero_quadratic.clone());
2289 let hyper_layout = realize_hyper_layout(&theta)?;
2290 let owned = evaluate_custom_family_joint_hyper_efs_owned(
2291 &family_eval,
2292 std::slice::from_ref(&blockspec),
2293 &outer_opts,
2294 &rho,
2295 &hyper_layout,
2296 st.warm_cache.as_ref(),
2297 )
2298 .map_err(|e| {
2299 EstimationError::RemlOptimizationFailed(format!(
2300 "spatial adaptive EFS eval failed: {e}"
2301 ))
2302 })?;
2303 if !owned.result.inner_converged {
2304 st.warm_cache = Some(owned.result.warm_start);
2305 return Err(EstimationError::RemlOptimizationFailed(
2306 "exact spatial adaptive EFS inner solve did not converge".to_string(),
2307 ));
2308 }
2309 let objective = owned.result.efs_eval.cost;
2310 st.warm_cache = Some(owned.result.warm_start);
2311 st.terminal_mode = Some((theta, objective, owned.mode));
2312 Ok(owned.result.efs_eval)
2313 }),
2314 |st: &mut SpatialAdaptiveOuterState, theta: &Array1<f64>| {
2326 let theta = clamp_theta(theta);
2327 let DecodedSpatialAdaptiveTheta {
2328 rho,
2329 adaptive_params,
2330 ..
2331 } = decode_theta(&theta)?;
2332 let family_eval =
2333 base_family.with_adaptive_params(adaptive_params, zero_quadratic.clone());
2334 let hyper_layout = realize_hyper_layout(&theta)?;
2335 let owned = evaluate_custom_family_joint_hyper_owned(
2336 &family_eval,
2337 std::slice::from_ref(&blockspec),
2338 &outer_opts,
2339 &rho,
2340 &hyper_layout,
2341 st.warm_cache.as_ref(),
2342 gam_solve::estimate::reml::reml_outer_engine::EvalMode::ValueOnly,
2343 )
2344 .map_err(|e| {
2345 EstimationError::RemlOptimizationFailed(format!(
2346 "spatial adaptive screening eval failed: {e}"
2347 ))
2348 })?;
2349 st.warm_cache = Some(owned.result.warm_start);
2350 Ok(owned.result.objective)
2351 },
2352 );
2353
2354 let certified_outer = problem
2355 .run_certified(&mut obj, "exact spatial adaptive regularization")
2356 .map_err(|e| {
2357 EstimationError::InvalidInput(format!(
2358 "exact spatial adaptive outer optimization failed: {e}"
2359 ))
2360 })?;
2361 let outer_iterations = certified_outer.iterations();
2362 let outer_grad_norm = certified_outer.final_grad_norm();
2363 let theta_star = certified_outer.rho().clone();
2364 let (mode_theta, mode_objective, terminal_mode) =
2365 obj.state.terminal_mode.take().ok_or_else(|| {
2366 EstimationError::InvalidInput(
2367 "exact spatial adaptive optimization certified without retaining its terminal coefficient mode"
2368 .to_string(),
2369 )
2370 })?;
2371 if mode_theta.len() != theta_star.len()
2372 || mode_theta
2373 .iter()
2374 .zip(theta_star.iter())
2375 .any(|(mode, certified)| mode.to_bits() != certified.to_bits())
2376 {
2377 return Err(EstimationError::InvalidInput(
2378 "exact spatial adaptive terminal coefficient mode does not bitwise match the certified hyperparameter vector"
2379 .to_string(),
2380 ));
2381 }
2382 if mode_objective.to_bits() != certified_outer.final_value().to_bits() {
2383 return Err(EstimationError::InvalidInput(format!(
2384 "exact spatial adaptive terminal coefficient mode objective does not bitwise match the certified objective: mode={mode_objective:.17e}, certified={:.17e}",
2385 certified_outer.final_value(),
2386 )));
2387 }
2388 let DecodedSpatialAdaptiveTheta {
2389 rho: _,
2390 retained_lambdas,
2391 adaptive_params,
2392 epsilon: eps_star,
2393 } = decode_theta(&theta_star)?;
2394 let mut fixed_total = Array2::<f64>::zeros((
2395 baseline.design.design.ncols(),
2396 baseline.design.design.ncols(),
2397 ));
2398 for (idx, penalty) in retained_penalties.iter().enumerate() {
2399 fixed_total.scaled_add(retained_lambdas[idx], penalty);
2400 }
2401 let certified_final_family = base_family.with_adaptive_params(
2407 adaptive_params.clone(),
2408 zero_quadratic.clone(),
2409 );
2410 let fixed_total = ValidatedFixedQuadraticHessian::try_from_dense(
2411 fixed_total,
2412 baseline.design.design.ncols(),
2413 )
2414 .map_err(|error| {
2415 EstimationError::InvalidInput(format!(
2416 "optimized spatial adaptive fixed quadratic Hessian is invalid: {error}"
2417 ))
2418 })?;
2419 let final_family =
2420 base_family.with_adaptive_params(adaptive_params.clone(), fixed_total.clone());
2421 let final_blockspec = ParameterBlockSpec {
2422 name: "eta".to_string(),
2423 design: baseline.design.design.clone(),
2424 offset: offset.to_owned(),
2425 penalties: retained_penalties
2426 .iter()
2427 .cloned()
2428 .map(PenaltyMatrix::Dense)
2429 .collect(),
2430 nullspace_dims: retained_nullspace_dims.clone(),
2431 initial_log_lambdas: theta_star.slice(s![..rho_dim]).to_owned(),
2432 initial_beta: Some(baseline.fit.beta.clone()),
2433 gauge_priority: 100,
2434 jacobian_callback: None,
2435 stacked_design: None,
2436 stacked_offset: None,
2437 };
2438 let final_fit = fit_custom_family_fixed_log_lambdas_from_owned_mode(
2439 &certified_final_family,
2440 &[final_blockspec],
2441 &BlockwiseFitOptions {
2442 inner_max_cycles: options.max_iter,
2443 inner_tol: options.tol,
2444 outer_max_iter: 1,
2445 outer_tol: options.tol,
2446 compute_covariance: true,
2447 ..BlockwiseFitOptions::default()
2448 },
2449 terminal_mode,
2450 &theta_star,
2451 &certified_outer,
2452 )
2453 .map_err(EstimationError::CustomFamily)?;
2454 let beta = final_fit.block_states[0].beta.clone();
2455 let final_eval = final_family
2456 .exact_evaluation(&beta)
2457 .map_err(EstimationError::InvalidInput)?;
2458 let penalized_hessian = final_eval
2459 .totalobjectivehessian(&final_family.design)
2460 .map_err(EstimationError::InvalidInput)?;
2461 let beta_covariance = final_fit.covariance_conditional.clone();
2462 let beta_standard_errors = beta_covariance
2463 .as_ref()
2464 .map(|cov| Array1::from_iter((0..cov.nrows()).map(|i| cov[[i, i]].max(0.0).sqrt())));
2465
2466 let mut full_lambdas = baseline.fit.lambdas.clone();
2467 for (idx, &global_idx) in retained_global_indices.iter().enumerate() {
2468 full_lambdas[global_idx] = retained_lambdas[idx];
2469 }
2470 for (cache_idx, cache) in runtime_caches.iter().enumerate() {
2471 full_lambdas[cache.mass_penalty_global_idx] = adaptive_params[cache_idx].lambda[0];
2472 full_lambdas[cache.tension_penalty_global_idx] = adaptive_params[cache_idx].lambda[1];
2473 full_lambdas[cache.stiffness_penalty_global_idx] = adaptive_params[cache_idx].lambda[2];
2474 }
2475
2476 let deviance = -2.0 * final_eval.obs.log_likelihood;
2477 let mut local_penalty_blocks =
2478 Vec::<PenaltySpec>::with_capacity(baseline.design.penalties.len());
2479 for (global_idx, bp) in baseline.design.penalties.iter().enumerate() {
2480 if adaptive_penalty_indices.contains(&global_idx) {
2481 let cache = runtime_caches
2482 .iter()
2483 .find(|cache| {
2484 cache.mass_penalty_global_idx == global_idx
2485 || cache.tension_penalty_global_idx == global_idx
2486 || cache.stiffness_penalty_global_idx == global_idx
2487 })
2488 .ok_or_else(|| {
2489 EstimationError::InvalidInput(format!(
2490 "missing runtime cache for adaptive penalty index {global_idx}"
2491 ))
2492 })?;
2493 let cache_idx = runtime_caches
2494 .iter()
2495 .position(|c| {
2496 c.mass_penalty_global_idx == global_idx
2497 || c.tension_penalty_global_idx == global_idx
2498 || c.stiffness_penalty_global_idx == global_idx
2499 })
2500 .ok_or_else(|| {
2501 EstimationError::InvalidInput(format!(
2502 "missing adaptive cache position for penalty index {global_idx}"
2503 ))
2504 })?;
2505 let state = &final_eval.adaptive_states[cache_idx];
2506 let local = if cache.mass_penalty_global_idx == global_idx {
2507 scalar_operatorhessian(&cache.d0, &state.magnitude.betahessian_diag())
2508 .mapv(|v| adaptive_params[cache_idx].lambda[0] * v)
2509 } else if cache.tension_penalty_global_idx == global_idx {
2510 grouped_operatorhessian(
2511 &cache.d1,
2512 cache.dimension,
2513 &state.gradient.betahessian_blocks(),
2514 )?
2515 .mapv(|v| adaptive_params[cache_idx].lambda[1] * v)
2516 } else {
2517 grouped_operatorhessian(
2518 &cache.d2,
2519 cache.dimension * cache.dimension,
2520 &state.curvature.betahessian_blocks(),
2521 )?
2522 .mapv(|v| adaptive_params[cache_idx].lambda[2] * v)
2523 };
2524 local_penalty_blocks.push(PenaltySpec::Dense(penalty_matrixwith_local_block(
2526 baseline.design.design.ncols(),
2527 cache.coeff_global_range.clone(),
2528 &local,
2529 )));
2530 } else {
2531 local_penalty_blocks.push(PenaltySpec::Dense(
2532 bp.to_global(p_total).mapv(|v| v * full_lambdas[global_idx]),
2533 ));
2534 }
2535 }
2536 let (edf_by_block, penalty_block_trace, edf_total) = if let Some(cov) = beta_covariance.as_ref()
2537 {
2538 exact_bounded_edf(
2539 &local_penalty_blocks,
2540 &Array1::from_elem(local_penalty_blocks.len(), 1.0),
2541 cov,
2542 )?
2543 } else {
2544 (
2545 vec![0.0; local_penalty_blocks.len()],
2546 vec![0.0; local_penalty_blocks.len()],
2547 0.0,
2548 )
2549 };
2550 let stable_penalty_term = 2.0 * final_eval.adaptive_penalty_value
2551 + beta.dot(&fixed_total.as_dense().dot(&beta));
2552 let standard_deviation = if family.is_gaussian_identity() {
2553 let denom = (y.len() as f64 - edf_total).max(1.0);
2554 (deviance / denom).sqrt()
2555 } else {
2556 1.0
2557 };
2558 let maps = compute_spatial_adaptiveweights_for_beta(
2559 &beta,
2560 runtime_caches,
2561 eps_star[0],
2562 eps_star[1],
2563 eps_star[2],
2564 adaptive_opts.weight_floor,
2565 adaptive_opts.weight_ceiling,
2566 beta_covariance.as_ref(),
2570 )?
2571 .into_iter()
2572 .zip(runtime_caches.iter())
2573 .map(|(w, cache)| AdaptiveSpatialMap {
2574 termname: cache.termname.clone(),
2575 feature_cols: cache.feature_cols.clone(),
2576 collocation_points: cache.collocation_points.clone(),
2577 inv_magweight: w.inv_magweight,
2578 invgradweight: w.invgradweight,
2579 inv_lapweight: w.inv_lapweight,
2580 })
2581 .collect::<Vec<_>>();
2582 let fitted_link = if family.is_latent_cloglog() {
2583 FittedLinkState::LatentCLogLog {
2584 state: latent_cloglog_state
2585 .expect("BinomialLatentCLogLog requires an explicit latent-cloglog state"),
2586 }
2587 } else if family.is_binomial_mixture() {
2588 mixture_link_state
2589 .clone()
2590 .map(|state| FittedLinkState::Mixture {
2591 state,
2592 covariance: None,
2593 })
2594 .unwrap_or(FittedLinkState::Standard(None))
2595 } else if family.is_binomial_sas() {
2596 sas_link_state
2597 .map(|state| FittedLinkState::Sas {
2598 state,
2599 covariance: None,
2600 })
2601 .unwrap_or(FittedLinkState::Standard(None))
2602 } else if family.is_binomial_beta_logistic() {
2603 sas_link_state
2604 .map(|state| FittedLinkState::BetaLogistic {
2605 state,
2606 covariance: None,
2607 })
2608 .unwrap_or(FittedLinkState::Standard(None))
2609 } else {
2610 FittedLinkState::Standard(None)
2611 };
2612 let max_abs_eta = final_eval
2613 .obs
2614 .eta
2615 .iter()
2616 .fold(0.0_f64, |acc, &v| acc.max(v.abs()));
2617 let fitted = FittedTermCollection {
2618 fit: {
2619 let log_lambdas =
2620 checked_fit_log_lambdas(&full_lambdas, "final exact spatial adaptive fit")?;
2621 let working = gam_solve::estimate::WorkingGeometry {
2622 weights: final_eval.obs.fisherweight.clone(),
2623 response: exact_standard_working_response(&final_eval.obs)?,
2624 };
2625 let inf = FitInference {
2626 edf_by_block,
2627 penalty_block_trace,
2628 edf_total,
2629 smoothing_correction: None,
2630 smoothing_correction_method: None,
2631 smoothing_correction_first_order: None,
2632 smoothing_correction_method_first_order: None,
2633 penalized_hessian: penalized_hessian.clone().into(),
2636 reparam_qs: None,
2637 dispersion: gam_solve::estimate::Dispersion::UNIT,
2638 beta_covariance: beta_covariance
2639 .clone()
2640 .map(gam_problem::dispersion_cov::PhiScaledCovariance::from),
2641 beta_standard_errors,
2642 beta_covariance_corrected: None,
2643 beta_standard_errors_corrected: None,
2644 beta_covariance_frequentist: None,
2645 coefficient_influence: None,
2646 weighted_gram: None,
2647 bias_correction_beta: None,
2648 bias_correction_jacobian: None,
2649 };
2650 let geometry = Some(gam_solve::estimate::FitGeometry {
2651 coefficient_gauge: gam_problem::gauge::Gauge::identity(&[beta.len()]),
2652 penalized_hessian: penalized_hessian.into(),
2653 constrained_posterior: None,
2654 working: Some(working),
2655 });
2656 let covariance_conditional = beta_covariance;
2657 let convergence = final_fit.convergence_evidence();
2658 let pirls_status_val = convergence.inner_status();
2659 let certified_outer_present = convergence.outer_certificate().is_some();
2660 UnifiedFitResult::try_from_parts(UnifiedFitResultParts {
2661 blocks: vec![gam_solve::estimate::FittedBlock {
2662 beta: beta.clone(),
2663 role: gam_problem::BlockRole::Mean,
2664 edf: edf_total,
2665 lambdas: full_lambdas.clone(),
2666 }],
2667 log_lambdas,
2668 lambdas: full_lambdas,
2669 likelihood_scale: family.default_scale_metadata(),
2670 likelihood_family: Some(family),
2671 log_likelihood_normalization: gam_spec::LogLikelihoodNormalization::UserProvided,
2672 log_likelihood: final_eval.obs.log_likelihood,
2673 deviance,
2674 reml_score: final_fit.penalized_objective,
2675 stable_penalty_term,
2676 penalized_objective: final_fit.penalized_objective,
2677 used_device: false,
2678 outer_iterations,
2679 outer_converged: certified_outer_present,
2680 outer_gradient_norm: outer_grad_norm,
2681 standard_deviation,
2682 covariance_conditional,
2683 covariance_corrected: None,
2684 inference: Some(inf),
2685 fitted_link,
2686 geometry,
2687 block_states: Vec::new(),
2688 pirls_status: pirls_status_val,
2689 max_abs_eta,
2690 constraint_kkt: None,
2691 artifacts: gam_solve::estimate::FitArtifacts {
2692 pirls: None,
2693 criterion_certificate: final_fit.artifacts.criterion_certificate.clone(),
2694 ..Default::default()
2695 },
2696 inner_cycles: 0,
2697 })?
2698 },
2699 design: baseline.design,
2700 adaptive_diagnostics: Some(AdaptiveRegularizationDiagnostics {
2701 epsilon_0: eps_star[0],
2702 epsilon_g: eps_star[1],
2703 epsilon_c: eps_star[2],
2704 epsilon_outer_iterations: outer_iterations,
2705 mm_iterations: 0,
2706 converged: true,
2707 maps,
2708 }),
2709 };
2710 enforce_term_constraint_feasibility(&fitted.design, &fitted.fit)?;
2711 Ok(fitted)
2712}
2713
2714fn relax_smoothing_rho_prior(
2754 options: &FitOptions,
2755 design: &TermCollectionDesign,
2756 y: ArrayView1<'_, f64>,
2757 weights: ArrayView1<'_, f64>,
2758) -> gam_spec::RhoPrior {
2759 use gam_terms::basis::BasisMetadata;
2760 let base = &options.rho_prior;
2761 if !base.is_unset() {
2789 return base.clone();
2790 }
2791 let has_link_aux = options.sas_link.is_some()
2811 || options.optimize_sas
2812 || options.mixture_link.is_some()
2813 || options.optimize_mixture;
2814 let has_moving_kappa = design.smooth.terms.iter().any(|t| {
2815 if let BasisMetadata::Duchon {
2828 length_scale,
2829 aniso_log_scales,
2830 ..
2831 } = &t.metadata
2832 {
2833 return length_scale.is_some() || aniso_log_scales.is_some();
2834 }
2835 matches!(
2836 t.metadata,
2837 BasisMetadata::Matern { .. }
2838 | BasisMetadata::Sphere { .. }
2839 | BasisMetadata::SphereHarmonics { .. }
2840 | BasisMetadata::ConstantCurvature { .. }
2841 | BasisMetadata::MeasureJet { .. }
2842 )
2843 });
2844 let length_safe = !has_link_aux && !has_moving_kappa;
2851 if !length_safe {
2852 return base.clone();
2853 }
2854 let coords = &design.penaltyinfo;
2855 if coords.is_empty() {
2856 return base.clone();
2857 }
2858 let n_obs = design.design.nrows();
2869 let p_total = design.design.ncols();
2870 let underdetermined = n_obs < 2 * p_total;
2901 let relaxable_terms: std::collections::HashSet<&str> = design
2913 .smooth
2914 .terms
2915 .iter()
2916 .filter(|t| {
2917 (matches!(
2918 t.metadata,
2919 BasisMetadata::BSpline1D { .. }
2920 | BasisMetadata::ThinPlate { .. }
2921 | BasisMetadata::TensorBSpline { .. }
2922 )
2923 || matches!(
2935 t.metadata,
2936 BasisMetadata::Duchon {
2937 length_scale: None,
2938 aniso_log_scales: None,
2939 ..
2940 }
2941 ))
2942 && matches!(t.shape, gam_terms::smooth::ShapeConstraint::None)
2956 })
2957 .map(|t| t.name.as_str())
2958 .collect();
2959 let any_relaxed = coords.iter().any(|info| {
2960 info.termname
2961 .as_deref()
2962 .is_some_and(|name| relaxable_terms.contains(name))
2963 });
2964 if !any_relaxed {
2965 return base.clone();
2966 }
2967 let relaxed_prior = if underdetermined {
2972 gam_spec::RhoPrior::Normal {
2973 mean: 0.0,
2974 sd: RELAX_UNDERDETERMINED_RHO_SD,
2975 }
2976 } else {
2977 gam_spec::RhoPrior::Flat
2978 };
2979 let nullspace_select_prior = gam_spec::RhoPrior::PenalizedComplexity {
3006 upper: NULLSPACE_SELECT_PC_UPPER,
3007 tail_prob: NULLSPACE_SELECT_PC_TAIL_PROB,
3008 };
3009 let nullspace_degeneracy_prior = gam_spec::RhoPrior::Normal {
3036 mean: 0.0,
3037 sd: NULLSPACE_WELLDET_DEGENERACY_RHO_SD,
3038 };
3039 let per_coord = coords
3040 .iter()
3041 .enumerate()
3042 .map(|(coord_idx, info)| {
3043 let relax = info
3044 .termname
3045 .as_deref()
3046 .is_some_and(|name| relaxable_terms.contains(name));
3047 if !relax {
3048 return base.clone();
3049 }
3050 let is_nullspace = matches!(info.penalty.source, PenaltySource::DoublePenaltyNullspace);
3051 if is_nullspace {
3090 if underdetermined
3111 && !nullspace_directions_are_supported(design, coord_idx, y, weights)
3112 {
3113 nullspace_select_prior.clone()
3114 } else {
3115 nullspace_degeneracy_prior.clone()
3116 }
3117 } else {
3118 relaxed_prior.clone()
3119 }
3120 })
3121 .collect::<Vec<_>>();
3122 gam_spec::RhoPrior::Independent(per_coord)
3123}
3124
3125const NULLSPACE_SUPPORT_FRACTION_THRESHOLD: f64 = 0.5;
3132
3133fn nullspace_directions_are_supported(
3155 design: &TermCollectionDesign,
3156 penalty_idx: usize,
3157 y: ArrayView1<'_, f64>,
3158 weights: ArrayView1<'_, f64>,
3159) -> bool {
3160 use gam_linalg::faer_ndarray::FaerEigh;
3161
3162 let Some(pen) = design.penalties.get(penalty_idx) else {
3163 return false;
3164 };
3165 let col_range = pen.col_range.clone();
3166 if col_range.is_empty() {
3167 return false;
3168 }
3169 let x = design.design.to_dense();
3170 let n = x.nrows();
3171 if n == 0 || y.len() != n || weights.len() != n || x.ncols() < col_range.end {
3172 return false;
3173 }
3174 let mut resp = y.to_owned();
3178 if design.affine_offset.len() == n {
3179 resp -= &design.affine_offset;
3180 }
3181
3182 let Ok((evals, evecs)) = pen.local.eigh(faer::Side::Lower) else {
3185 return false;
3186 };
3187 let max_eig = evals.iter().cloned().fold(0.0_f64, |m, v| m.max(v));
3188 if !(max_eig > 0.0) {
3189 return false;
3190 }
3191 let tol = 1.0e-9 * max_eig;
3192 let pos_cols: Vec<usize> = (0..evals.len()).filter(|&j| evals[j] > tol).collect();
3193 if pos_cols.is_empty() {
3194 return false;
3195 }
3196 let xblock = x.slice(s![.., col_range.clone()]);
3197 let mut z = Array2::<f64>::zeros((n, pos_cols.len()));
3198 for (out_j, &j) in pos_cols.iter().enumerate() {
3199 let v = evecs.column(j);
3200 if v.len() != xblock.ncols() {
3202 return false;
3203 }
3204 z.column_mut(out_j).assign(&xblock.dot(&v));
3205 }
3206
3207 let mut control: Vec<usize> = design.intercept_range.clone().collect();
3212 for (_, r) in &design.linear_ranges {
3213 control.extend(r.clone());
3214 }
3215 control.retain(|&c| c < x.ncols());
3216 let mut cmat = Array2::<f64>::zeros((n, control.len().max(1)));
3217 if control.is_empty() {
3218 cmat.column_mut(0).fill(1.0);
3220 } else {
3221 for (out_j, &c) in control.iter().enumerate() {
3222 cmat.column_mut(out_j).assign(&x.column(c));
3223 }
3224 }
3225
3226 let rss_c = weighted_regression_rss(cmat.view(), resp.view(), weights);
3227 let Some(rss_c) = rss_c else { return false };
3228 let base_scale = weighted_total_ss(resp.view(), weights);
3231 if !(rss_c > 1.0e-12 * base_scale.max(f64::MIN_POSITIVE)) {
3232 return false;
3233 }
3234 let mut cz = Array2::<f64>::zeros((n, cmat.ncols() + z.ncols()));
3235 cz.slice_mut(s![.., ..cmat.ncols()]).assign(&cmat);
3236 cz.slice_mut(s![.., cmat.ncols()..]).assign(&z);
3237 let Some(rss_cz) = weighted_regression_rss(cz.view(), resp.view(), weights) else {
3238 return false;
3239 };
3240
3241 let support = (rss_c - rss_cz) / rss_c;
3242 support.is_finite() && support > NULLSPACE_SUPPORT_FRACTION_THRESHOLD
3243}
3244
3245fn weighted_total_ss(y: ArrayView1<'_, f64>, w: ArrayView1<'_, f64>) -> f64 {
3247 let mut sw = 0.0;
3248 let mut swy = 0.0;
3249 for (&yi, &wi) in y.iter().zip(w.iter()) {
3250 if wi > 0.0 && yi.is_finite() {
3251 sw += wi;
3252 swy += wi * yi;
3253 }
3254 }
3255 if sw <= 0.0 {
3256 return 0.0;
3257 }
3258 let mean = swy / sw;
3259 let mut ss = 0.0;
3260 for (&yi, &wi) in y.iter().zip(w.iter()) {
3261 if wi > 0.0 && yi.is_finite() {
3262 ss += wi * (yi - mean) * (yi - mean);
3263 }
3264 }
3265 ss
3266}
3267
3268fn weighted_regression_rss(
3275 d: ArrayView2<'_, f64>,
3276 y: ArrayView1<'_, f64>,
3277 w: ArrayView1<'_, f64>,
3278) -> Option<f64> {
3279 use gam_linalg::faer_ndarray::FaerCholesky;
3280
3281 let m = d.ncols();
3282 if m == 0 {
3283 return Some(weighted_total_ss(y, w));
3284 }
3285 let mut gram = Array2::<f64>::zeros((m, m));
3286 let mut rhs = Array1::<f64>::zeros(m);
3287 for row in 0..d.nrows() {
3288 let wi = w[row];
3289 if !(wi > 0.0) || !y[row].is_finite() {
3290 continue;
3291 }
3292 let dr = d.row(row);
3293 for a in 0..m {
3294 let wda = wi * dr[a];
3295 rhs[a] += wda * y[row];
3296 for b in a..m {
3297 gram[[a, b]] += wda * dr[b];
3298 }
3299 }
3300 }
3301 for a in 0..m {
3302 for b in (a + 1)..m {
3303 gram[[b, a]] = gram[[a, b]];
3304 }
3305 }
3306 let trace = (0..m).map(|i| gram[[i, i]]).sum::<f64>();
3307 if !(trace > 0.0) {
3308 return Some(weighted_total_ss(y, w));
3309 }
3310 let ridge = 1.0e-10 * trace / (m as f64);
3311 for i in 0..m {
3312 gram[[i, i]] += ridge;
3313 }
3314 let chol = gram.cholesky(faer::Side::Lower).ok()?;
3315 let beta = chol.solvevec(&rhs);
3316 let mut rss = 0.0;
3317 for row in 0..d.nrows() {
3318 let wi = w[row];
3319 if !(wi > 0.0) || !y[row].is_finite() {
3320 continue;
3321 }
3322 let fitted = d.row(row).dot(&beta);
3323 let resid = y[row] - fitted;
3324 rss += wi * resid * resid;
3325 }
3326 Some(rss)
3327}
3328
3329const RELAX_UNDERDETERMINED_RHO_SD: f64 = 15.0;
3342
3343const NULLSPACE_SELECT_PC_UPPER: f64 = 0.05;
3361
3362const NULLSPACE_SELECT_PC_TAIL_PROB: f64 = 0.01;
3372
3373fn adaptive_fit_options_base(options: &FitOptions, design: &TermCollectionDesign) -> FitOptions {
3374 FitOptions {
3375 resource_policy: options.resource_policy.clone(),
3376 latent_cloglog: options.latent_cloglog,
3377 mixture_link: options.mixture_link.clone(),
3378 optimize_mixture: options.optimize_mixture,
3379 sas_link: options.sas_link,
3380 optimize_sas: options.optimize_sas,
3381 compute_inference: options.compute_inference,
3382 skip_rho_posterior_inference: options.skip_rho_posterior_inference,
3383 max_iter: options.max_iter,
3384 tol: options.tol,
3385 nullspace_dims: design.nullspace_dims.clone(),
3386 linear_constraints: design.linear_constraints.clone(),
3387 firth_bias_reduction: options.firth_bias_reduction,
3388 adaptive_regularization: None,
3389 penalty_shrinkage_floor: options.penalty_shrinkage_floor,
3390 rho_prior: options.rho_prior.clone(),
3393 kronecker_penalty_system: design.kronecker_penalty_system(),
3394 kronecker_factored: design
3395 .smooth
3396 .terms
3397 .iter()
3398 .find_map(|t| t.kronecker_factored.clone()),
3399 persist_warm_start_disk: options.persist_warm_start_disk,
3400 }
3401}
3402
3403fn superseded_fit_options(options: &FitOptions) -> FitOptions {
3404 let mut fit_options = options.clone();
3405 fit_options.skip_rho_posterior_inference = true;
3406 fit_options
3407}
3408
3409#[derive(Clone)]
3410struct BoundedLinearTermMeta {
3411 col_idx: usize,
3412 min: f64,
3413 max: f64,
3414 prior: BoundedCoefficientPriorSpec,
3415}
3416
3417struct BoundedEffectiveJacobian {
3441 design: Array2<f64>,
3442 bounded_terms: Vec<BoundedLinearTermMeta>,
3443}
3444
3445impl BlockEffectiveJacobian for BoundedEffectiveJacobian {
3446 fn effective_jacobian_rows(
3447 &self,
3448 state: &FamilyLinearizationState<'_>,
3449 rows: std::ops::Range<usize>,
3450 ) -> Result<Array2<f64>, String> {
3451 let p = self.design.ncols();
3452 let n = self.design.nrows();
3453 let rows = rows.start.min(n)..rows.end.min(n);
3454 if !state.beta.is_empty() {
3455 if state.beta.len() != p {
3456 return Err(format!(
3457 "BoundedEffectiveJacobian::effective_jacobian_at: beta length {} != design \
3458 ncols {p}",
3459 state.beta.len(),
3460 ));
3461 }
3462 if state.beta.iter().any(|v| !v.is_finite()) {
3463 return Err(
3464 "BoundedEffectiveJacobian::effective_jacobian_at: beta contains a non-finite value"
3465 .to_string(),
3466 );
3467 }
3468 }
3469 let mut jac = self
3470 .design
3471 .slice(ndarray::s![rows.start..rows.end, ..])
3472 .to_owned();
3473 for term in &self.bounded_terms {
3474 if term.col_idx >= p {
3475 return Err(format!(
3476 "BoundedEffectiveJacobian::effective_jacobian_at: bounded column {} is outside {p} columns",
3477 term.col_idx
3478 ));
3479 }
3480 let theta = if state.beta.is_empty() {
3481 0.0
3482 } else {
3483 state.beta[term.col_idx]
3484 };
3485 let (_, _, db_dtheta, _, _) = bounded_latent_derivatives(theta, term.min, term.max);
3486 if !(db_dtheta.is_finite() && db_dtheta > 0.0) {
3487 return Err(format!(
3488 "BoundedEffectiveJacobian::effective_jacobian_at: bounded column {} has unrepresentable derivative {db_dtheta} at theta={theta}",
3489 term.col_idx
3490 ));
3491 }
3492 jac.column_mut(term.col_idx).mapv_inplace(|v| v * db_dtheta);
3493 }
3494 Ok(jac)
3495 }
3496}
3497
3498#[derive(Clone)]
3499struct BoundedLinearFamily {
3500 likelihood: gam_spec::GlmLikelihoodSpec,
3501 latent_cloglog_state: Option<LatentCLogLogState>,
3502 mixture_link_state: Option<MixtureLinkState>,
3503 sas_link_state: Option<SasLinkState>,
3504 y: Array1<f64>,
3505 weights: Array1<f64>,
3506 design: Array2<f64>,
3507 designzeroed: Array2<f64>,
3508 offset: Array1<f64>,
3509 bounded_terms: Vec<BoundedLinearTermMeta>,
3510}
3511
3512#[derive(Clone, Debug)]
3513struct StandardFamilyObservationState {
3514 eta: Array1<f64>,
3515 score: Array1<f64>,
3516 fisherweight: Array1<f64>,
3517 neghessian_eta: Array1<f64>,
3518 neghessian_eta_derivative: Array1<f64>,
3519 log_likelihood: f64,
3520}
3521
3522fn bounded_latent_to_user(theta: f64, min: f64, max: f64) -> (f64, f64, f64) {
3523 let jet = logit_inverse_link_jet5(theta);
3524 let z = jet.mu;
3525 let width = max - min;
3526 let beta = min + width * z;
3527 let db_dtheta = width * jet.d1;
3528 (beta, z, db_dtheta)
3529}
3530
3531fn bounded_user_to_latent(beta: f64, min: f64, max: f64) -> f64 {
3540 (beta - min).ln() - (max - beta).ln()
3541}
3542
3543#[derive(Debug, Clone, Copy)]
3547pub struct BoundedSampleColumn {
3548 pub col_idx: usize,
3550 pub min: f64,
3552 pub max: f64,
3554}
3555
3556pub fn sample_bounded_latent_posterior_internal(
3594 beta_user: &Array1<f64>,
3595 user_hessian: &Array2<f64>,
3596 bounded_columns: &[BoundedSampleColumn],
3597 n_draws: usize,
3598 sqrt_cov_scale: f64,
3599 base_seed: u64,
3600) -> Result<Array2<f64>, EstimationError> {
3601 let p = beta_user.len();
3602 if user_hessian.nrows() != p || user_hessian.ncols() != p {
3603 crate::bail_invalid_estim!(
3604 "bounded posterior sampling dimension mismatch: mode has {p} entries, user Hessian is {}x{}",
3605 user_hessian.nrows(),
3606 user_hessian.ncols()
3607 );
3608 }
3609 if beta_user.iter().any(|value| !value.is_finite()) {
3610 crate::bail_invalid_estim!("bounded posterior sampling requires a finite mode");
3611 }
3612 if user_hessian.iter().any(|value| !value.is_finite()) {
3613 crate::bail_invalid_estim!("bounded posterior sampling requires a finite Hessian");
3614 }
3615 if !(sqrt_cov_scale.is_finite() && sqrt_cov_scale >= 0.0) {
3616 crate::bail_invalid_estim!(
3617 "bounded posterior sampling covariance scale must be finite and non-negative, got {sqrt_cov_scale}"
3618 );
3619 }
3620
3621 let mut theta_mode = beta_user.clone();
3623 let mut jac_diag = Array1::<f64>::ones(p);
3624 for bc in bounded_columns {
3625 if bc.col_idx >= p {
3626 crate::bail_invalid_estim!(
3627 "bounded posterior sampling: bounded column index {} out of range for {p} coefficients",
3628 bc.col_idx
3629 );
3630 }
3631 if !(bc.min.is_finite()
3632 && bc.max.is_finite()
3633 && (bc.max - bc.min).is_finite()
3634 && bc.min < beta_user[bc.col_idx]
3635 && beta_user[bc.col_idx] < bc.max)
3636 {
3637 crate::bail_invalid_estim!(
3638 "bounded posterior sampling column {} requires finite bounds with a finite width and a mode strictly inside ({}, {}); got {}",
3639 bc.col_idx,
3640 bc.min,
3641 bc.max,
3642 beta_user[bc.col_idx]
3643 );
3644 }
3645 let theta_i = bounded_user_to_latent(beta_user[bc.col_idx], bc.min, bc.max);
3646 let (_, _, db_dtheta) = bounded_latent_to_user(theta_i, bc.min, bc.max);
3647 if !(theta_i.is_finite() && db_dtheta.is_finite() && db_dtheta > 0.0) {
3648 crate::bail_invalid_estim!(
3649 "bounded posterior sampling column {} has unrepresentable latent geometry: theta={theta_i}, d_beta/d_theta={db_dtheta}",
3650 bc.col_idx
3651 );
3652 }
3653 theta_mode[bc.col_idx] = theta_i;
3654 jac_diag[bc.col_idx] = db_dtheta;
3655 }
3656
3657 let mut h_latent = user_hessian.clone();
3660 for i in 0..p {
3661 let ji = jac_diag[i];
3662 if ji != 1.0 {
3663 h_latent.row_mut(i).mapv_inplace(|v| v * ji);
3664 h_latent.column_mut(i).mapv_inplace(|v| v * ji);
3665 }
3666 }
3667
3668 use gam_linalg::faer_ndarray::FaerCholesky as _;
3671 use rand::SeedableRng as _;
3672 let chol = h_latent.cholesky(faer::Side::Lower).map_err(|err| {
3673 EstimationError::InvalidInput(format!(
3674 "bounded posterior sampling: Cholesky of the latent penalized Hessian failed: {err:?}"
3675 ))
3676 })?;
3677 let l = chol.lower_triangular();
3678
3679 let mut draws = Array2::<f64>::zeros((n_draws, p));
3680 let mut eps = Array1::<f64>::zeros(p);
3681 let mut delta = Array1::<f64>::zeros(p);
3682 let mut rng = rand::rngs::StdRng::seed_from_u64(base_seed);
3683 for k in 0..n_draws {
3684 for e in eps.iter_mut() {
3685 *e = standard_normal_draw(&mut rng);
3686 }
3687 solve_lower_transpose_into(&l, &eps, &mut delta)?;
3688 for i in 0..p {
3689 draws[(k, i)] = theta_mode[i] + sqrt_cov_scale * delta[i];
3692 }
3693 for bc in bounded_columns {
3698 let (beta_draw, _, _) = bounded_latent_to_user(draws[(k, bc.col_idx)], bc.min, bc.max);
3699 draws[(k, bc.col_idx)] = beta_draw;
3700 }
3701 }
3702
3703 Ok(draws)
3704}
3705
3706#[inline]
3709fn standard_normal_draw<R: rand::Rng + ?Sized>(rng: &mut R) -> f64 {
3710 use rand::RngExt as _;
3711 let u1 = loop {
3712 let candidate = rng.random::<f64>();
3713 if candidate > 0.0 {
3714 break candidate;
3715 }
3716 };
3717 let u2 = rng.random::<f64>();
3718 (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
3719}
3720
3721fn solve_lower_transpose_into(
3725 l: &Array2<f64>,
3726 b: &Array1<f64>,
3727 out: &mut Array1<f64>,
3728) -> Result<(), EstimationError> {
3729 let p = l.nrows();
3730 if l.ncols() != p || b.len() != p || out.len() != p {
3731 crate::bail_invalid_estim!(
3732 "bounded triangular solve dimension mismatch: L={}x{}, b={}, out={}",
3733 l.nrows(),
3734 l.ncols(),
3735 b.len(),
3736 out.len()
3737 );
3738 }
3739 for i in (0..p).rev() {
3740 let mut acc = b[i];
3741 for j in (i + 1)..p {
3742 acc -= l[(j, i)] * out[j];
3743 }
3744 let diag = l[(i, i)];
3745 if !(diag.is_finite() && diag > 0.0 && acc.is_finite()) {
3746 crate::bail_invalid_estim!(
3747 "bounded triangular solve has invalid row {i}: diagonal={diag}, residual={acc}"
3748 );
3749 }
3750 let value = acc / diag;
3751 if !value.is_finite() {
3752 crate::bail_invalid_estim!(
3753 "bounded triangular solve produced a non-finite value at row {i}: {acc}/{diag}"
3754 );
3755 }
3756 out[i] = value;
3757 }
3758 Ok(())
3759}
3760
3761fn bounded_latent_derivatives(theta: f64, min: f64, max: f64) -> (f64, f64, f64, f64, f64) {
3762 let jet = logit_inverse_link_jet5(theta);
3763 let z = jet.mu;
3764 let width = max - min;
3765 let beta = min + width * z;
3766 let db_dtheta = width * jet.d1;
3767 let d2b_dtheta2 = width * jet.d2;
3768 let d3b_dtheta3 = width * jet.d3;
3769 (beta, z, db_dtheta, d2b_dtheta2, d3b_dtheta3)
3770}
3771
3772fn bounded_prior_terms(
3773 theta: f64,
3774 prior: &BoundedCoefficientPriorSpec,
3775) -> Result<(f64, f64, f64, f64), String> {
3776 if !theta.is_finite() {
3777 return Err(format!(
3778 "bounded coefficient prior requires a finite latent coordinate, got {theta}"
3779 ));
3780 }
3781 let (a, b) = match prior {
3782 BoundedCoefficientPriorSpec::None => return Ok((0.0, 0.0, 0.0, 0.0)),
3784 BoundedCoefficientPriorSpec::Uniform => (1.0, 1.0),
3787 BoundedCoefficientPriorSpec::Beta { a, b } => (*a, *b),
3788 };
3789 if !(a.is_finite() && a > 0.0 && b.is_finite() && b > 0.0) {
3790 return Err(format!(
3791 "bounded coefficient Beta prior requires finite positive shapes, got ({a}, {b})"
3792 ));
3793 }
3794 let jet = logit_inverse_link_jet5(theta);
3795 let z = jet.mu;
3796 let logp = -a * gam_linalg::utils::stable_softplus(-theta)
3801 - b * gam_linalg::utils::stable_softplus(theta);
3802 let grad = a - (a + b) * z;
3803 let neghess = (a + b) * jet.d1;
3804 let neghess_derivative = (a + b) * jet.d2;
3805 let terms = (logp, grad, neghess, neghess_derivative);
3806 if [terms.0, terms.1, terms.2, terms.3]
3807 .iter()
3808 .any(|value| !value.is_finite())
3809 {
3810 return Err(format!(
3811 "bounded coefficient prior geometry is not representable at theta={theta}: {terms:?}"
3812 ));
3813 }
3814 Ok(terms)
3815}
3816
3817#[derive(Clone, Copy)]
3818struct ExactStandardObservationRow {
3819 mu: f64,
3820 score: f64,
3821 fisherweight: f64,
3822 neghessian_eta: f64,
3823 neghessian_eta_derivative: f64,
3824 log_likelihood: f64,
3825}
3826
3827impl ExactStandardObservationRow {
3828 #[inline]
3829 fn zero_weight(mu: f64) -> Self {
3830 Self {
3831 mu,
3832 score: 0.0,
3833 fisherweight: 0.0,
3834 neghessian_eta: 0.0,
3835 neghessian_eta_derivative: 0.0,
3836 log_likelihood: 0.0,
3837 }
3838 }
3839}
3840
3841#[inline]
3842fn bounded_row_error(row: usize, quantity: &'static str, eta: f64, value: f64) -> EstimationError {
3843 EstimationError::PirlsRowGeometryUnrepresentable {
3844 row,
3845 quantity,
3846 eta,
3847 value,
3848 }
3849}
3850
3851#[inline]
3852fn certify_bounded_row(
3853 row: usize,
3854 eta: f64,
3855 state: ExactStandardObservationRow,
3856) -> Result<ExactStandardObservationRow, EstimationError> {
3857 for (quantity, value) in [
3858 ("bounded-family mean", state.mu),
3859 ("bounded-family score", state.score),
3860 ("bounded-family Fisher weight", state.fisherweight),
3861 ("bounded-family observed Hessian", state.neghessian_eta),
3862 (
3863 "bounded-family observed Hessian derivative",
3864 state.neghessian_eta_derivative,
3865 ),
3866 ("bounded-family log likelihood", state.log_likelihood),
3867 ] {
3868 if !value.is_finite() {
3869 return Err(bounded_row_error(row, quantity, eta, value));
3870 }
3871 }
3872 if state.fisherweight < 0.0 {
3873 return Err(bounded_row_error(
3874 row,
3875 "bounded-family Fisher weight",
3876 eta,
3877 state.fisherweight,
3878 ));
3879 }
3880 Ok(state)
3881}
3882
3883#[inline]
3884fn weighted_positive_from_log(weight: f64, log_value: f64) -> f64 {
3885 if weight == 0.0 {
3886 return 0.0;
3887 }
3888 (weight.ln() + log_value).exp()
3889}
3890
3891#[inline]
3892fn weighted_product3(a: f64, b: f64, c: f64) -> f64 {
3893 crate::gamlss::scaled_signed_product3(a, b, c)
3894}
3895
3896#[inline]
3897fn convex_combination(y: f64, left: f64, right: f64) -> f64 {
3898 if y == 0.0 {
3899 right
3900 } else if y == 1.0 {
3901 left
3902 } else {
3903 y.mul_add(left, (1.0 - y) * right)
3904 }
3905}
3906
3907#[derive(Clone, Copy)]
3913struct BernoulliNaturalJet {
3914 mu: f64,
3915 log_mu: [f64; 4],
3916 log_one_minus_mu: [f64; 4],
3917 log_fisher: f64,
3918}
3919
3920#[inline]
3921fn probit_natural_jet(eta: f64) -> BernoulliNaturalJet {
3922 let left = gam_math::probability::normal_logcdf_derivatives(eta);
3923 let right_at_neg_eta = gam_math::probability::normal_logcdf_derivatives(-eta);
3924 let log_pdf = if eta.abs() <= f64::MAX.sqrt() {
3925 -0.5 * eta * eta - 0.5 * (2.0 * std::f64::consts::PI).ln()
3926 } else {
3927 f64::NEG_INFINITY
3928 };
3929 BernoulliNaturalJet {
3930 mu: left[0].exp(),
3931 log_mu: [left[0], left[1], left[2], left[3]],
3932 log_one_minus_mu: [
3933 right_at_neg_eta[0],
3934 -right_at_neg_eta[1],
3935 right_at_neg_eta[2],
3936 -right_at_neg_eta[3],
3937 ],
3938 log_fisher: 2.0 * log_pdf - left[0] - right_at_neg_eta[0],
3939 }
3940}
3941
3942#[inline]
3943fn cloglog_natural_jet(eta: f64) -> BernoulliNaturalJet {
3944 let x = eta.exp();
3945 if x == f64::INFINITY {
3946 return BernoulliNaturalJet {
3947 mu: 1.0,
3948 log_mu: [0.0; 4],
3949 log_one_minus_mu: [f64::NEG_INFINITY; 4],
3950 log_fisher: f64::NEG_INFINITY,
3951 };
3952 }
3953 if x == 0.0 {
3954 return BernoulliNaturalJet {
3955 mu: 0.0,
3956 log_mu: [eta, 1.0, 0.0, 0.0],
3957 log_one_minus_mu: [0.0; 4],
3958 log_fisher: eta,
3959 };
3960 }
3961 let mu = -(-x).exp_m1();
3962 let log_mu = if x < 0.5 {
3963 eta + (mu / x).ln()
3964 } else {
3965 mu.ln()
3966 };
3967 let h = if x < 1.0 {
3968 x / x.exp_m1()
3969 } else {
3970 let exp_neg_x = (-x).exp();
3971 x * exp_neg_x / (1.0 - exp_neg_x)
3972 };
3973 let a = 1.0 - x - h;
3974 let d2_log_mu = h * a;
3975 let d3_log_mu = h * (a * a - x - h * a);
3976 BernoulliNaturalJet {
3977 mu,
3978 log_mu: [log_mu, h, d2_log_mu, d3_log_mu],
3979 log_one_minus_mu: [-x, -x, -x, -x],
3980 log_fisher: 2.0 * eta - x - log_mu,
3981 }
3982}
3983
3984#[inline]
3985fn loglog_natural_jet(eta: f64) -> BernoulliNaturalJet {
3986 let mirrored = cloglog_natural_jet(-eta);
3987 BernoulliNaturalJet {
3988 mu: mirrored.log_one_minus_mu[0].exp(),
3989 log_mu: [
3990 mirrored.log_one_minus_mu[0],
3991 -mirrored.log_one_minus_mu[1],
3992 mirrored.log_one_minus_mu[2],
3993 -mirrored.log_one_minus_mu[3],
3994 ],
3995 log_one_minus_mu: [
3996 mirrored.log_mu[0],
3997 -mirrored.log_mu[1],
3998 mirrored.log_mu[2],
3999 -mirrored.log_mu[3],
4000 ],
4001 log_fisher: mirrored.log_fisher,
4002 }
4003}
4004
4005#[inline]
4006fn cauchit_natural_jet(eta: f64) -> BernoulliNaturalJet {
4007 let (mu, one_minus_mu) = if eta > 0.0 {
4008 let q = (eta.recip()).atan() / std::f64::consts::PI;
4009 (1.0 - q, q)
4010 } else if eta < 0.0 {
4011 let p = (-eta.recip()).atan() / std::f64::consts::PI;
4012 (p, 1.0 - p)
4013 } else {
4014 (0.5, 0.5)
4015 };
4016 let abs_eta = eta.abs();
4017 let log_one_plus_eta_sq = if abs_eta <= f64::MAX.sqrt() {
4018 (eta * eta).ln_1p()
4019 } else {
4020 2.0 * abs_eta.ln() + eta.recip().powi(2).ln_1p()
4021 };
4022 let log_d1 = -std::f64::consts::PI.ln() - log_one_plus_eta_sq;
4023 let ratio = if abs_eta <= 1.0 {
4024 eta / (1.0 + eta * eta)
4025 } else {
4026 1.0 / (eta + eta.recip())
4027 };
4028 let d2_over_d1 = -2.0 * ratio;
4029 let inv_one_plus_sq = if abs_eta <= 1.0 {
4030 1.0 / (1.0 + eta * eta)
4031 } else {
4032 let inv = eta.recip();
4033 inv * inv / (1.0 + inv * inv)
4034 };
4035 let d3_over_d1 = inv_one_plus_sq * (6.0 * (eta * ratio) - 2.0 * inv_one_plus_sq);
4036 let d1_over_mu = (log_d1 - mu.ln()).exp();
4037 let d1_over_q = (log_d1 - one_minus_mu.ln()).exp();
4038 let left_d2_ratio = d2_over_d1 * d1_over_mu;
4039 let right_d2_ratio = d2_over_d1 * d1_over_q;
4040 BernoulliNaturalJet {
4041 mu,
4042 log_mu: [
4043 mu.ln(),
4044 d1_over_mu,
4045 left_d2_ratio - d1_over_mu * d1_over_mu,
4046 d3_over_d1 * d1_over_mu - 3.0 * d1_over_mu * left_d2_ratio + 2.0 * d1_over_mu.powi(3),
4047 ],
4048 log_one_minus_mu: [
4049 one_minus_mu.ln(),
4050 -d1_over_q,
4051 -right_d2_ratio - d1_over_q * d1_over_q,
4052 -d3_over_d1 * d1_over_q - 3.0 * d1_over_q * right_d2_ratio - 2.0 * d1_over_q.powi(3),
4053 ],
4054 log_fisher: 2.0 * log_d1 - mu.ln() - one_minus_mu.ln(),
4055 }
4056}
4057
4058#[inline]
4059fn generic_bernoulli_natural_jet(
4060 row: usize,
4061 eta: f64,
4062 link: &InverseLink,
4063) -> Result<BernoulliNaturalJet, EstimationError> {
4064 let jet = inverse_link_jet_for_inverse_link(link, eta)?;
4065 if !(jet.mu.is_finite()
4066 && jet.mu > 0.0
4067 && jet.mu < 1.0
4068 && jet.d1.is_finite()
4069 && jet.d1 > 0.0
4070 && jet.d2.is_finite()
4071 && jet.d3.is_finite())
4072 {
4073 return Err(bounded_row_error(
4074 row,
4075 "bounded-family inverse-link jet",
4076 eta,
4077 jet.mu,
4078 ));
4079 }
4080 let mu = jet.mu;
4081 let q = 1.0 - mu;
4082 let r1 = jet.d1 / mu;
4083 let r2 = jet.d2 / mu;
4084 let r3 = jet.d3 / mu;
4085 let s1 = jet.d1 / q;
4086 let s2 = jet.d2 / q;
4087 let s3 = jet.d3 / q;
4088 Ok(BernoulliNaturalJet {
4089 mu,
4090 log_mu: [
4091 mu.ln(),
4092 r1,
4093 r2 - r1 * r1,
4094 r3 - 3.0 * r1 * r2 + 2.0 * r1.powi(3),
4095 ],
4096 log_one_minus_mu: [
4097 (-mu).ln_1p(),
4098 -s1,
4099 -s2 - s1 * s1,
4100 -s3 - 3.0 * s1 * s2 - 2.0 * s1.powi(3),
4101 ],
4102 log_fisher: 2.0 * jet.d1.ln() - mu.ln() - q.ln(),
4103 })
4104}
4105
4106fn resolved_bounded_binomial_link(
4107 family: &LikelihoodSpec,
4108 latent_cloglog_state: Option<&LatentCLogLogState>,
4109 mixture_link_state: Option<&MixtureLinkState>,
4110 sas_link_state: Option<&SasLinkState>,
4111) -> InverseLink {
4112 match &family.link {
4113 InverseLink::LatentCLogLog(_) => latent_cloglog_state
4114 .copied()
4115 .map(InverseLink::LatentCLogLog)
4116 .unwrap_or_else(|| family.link.clone()),
4117 InverseLink::Mixture(_) => mixture_link_state
4118 .cloned()
4119 .map(InverseLink::Mixture)
4120 .unwrap_or_else(|| family.link.clone()),
4121 InverseLink::Sas(_) => sas_link_state
4122 .copied()
4123 .map(InverseLink::Sas)
4124 .unwrap_or_else(|| family.link.clone()),
4125 InverseLink::BetaLogistic(_) => sas_link_state
4126 .copied()
4127 .map(InverseLink::BetaLogistic)
4128 .unwrap_or_else(|| family.link.clone()),
4129 InverseLink::Standard(_) => family.link.clone(),
4130 }
4131}
4132
4133fn binomial_natural_jet(
4134 row: usize,
4135 eta: f64,
4136 link: &InverseLink,
4137) -> Result<BernoulliNaturalJet, EstimationError> {
4138 match link {
4139 InverseLink::Standard(StandardLink::Probit) => Ok(probit_natural_jet(eta)),
4140 InverseLink::Standard(StandardLink::CLogLog) => Ok(cloglog_natural_jet(eta)),
4141 InverseLink::Standard(StandardLink::LogLog) => Ok(loglog_natural_jet(eta)),
4142 InverseLink::Standard(StandardLink::Cauchit) => Ok(cauchit_natural_jet(eta)),
4143 _ => generic_bernoulli_natural_jet(row, eta, link),
4144 }
4145}
4146
4147fn exact_logit_observation_row(
4148 row: usize,
4149 y: f64,
4150 weight: f64,
4151 eta: f64,
4152) -> Result<ExactStandardObservationRow, EstimationError> {
4153 let tail = (-eta.abs()).exp();
4154 let (mu, one_minus_mu) = if eta >= 0.0 {
4155 let q = tail / (1.0 + tail);
4156 (1.0 - q, q)
4157 } else {
4158 let p = tail / (1.0 + tail);
4159 (p, 1.0 - p)
4160 };
4161 if weight == 0.0 {
4162 return Ok(ExactStandardObservationRow::zero_weight(mu));
4163 }
4164 let log_fisher =
4165 -gam_linalg::utils::stable_softplus(eta) - gam_linalg::utils::stable_softplus(-eta);
4166 let fisherweight = weighted_positive_from_log(weight, log_fisher);
4167 if !(fisherweight.is_finite() && fisherweight > 0.0) {
4168 return Err(bounded_row_error(
4169 row,
4170 "bounded logit Fisher weight",
4171 eta,
4172 fisherweight,
4173 ));
4174 }
4175 let residual = if eta >= 0.0 {
4176 if y == 1.0 {
4177 one_minus_mu
4178 } else {
4179 (y - 1.0) + one_minus_mu
4180 }
4181 } else {
4182 y - mu
4183 };
4184 let log_likelihood_unit = if eta >= 0.0 {
4185 -(1.0 - y) * eta - gam_linalg::utils::stable_softplus(-eta)
4186 } else {
4187 y * eta - gam_linalg::utils::stable_softplus(eta)
4188 };
4189 certify_bounded_row(
4190 row,
4191 eta,
4192 ExactStandardObservationRow {
4193 mu,
4194 score: weight * residual,
4195 fisherweight,
4196 neghessian_eta: fisherweight,
4197 neghessian_eta_derivative: fisherweight * (one_minus_mu - mu),
4198 log_likelihood: weight * log_likelihood_unit,
4199 },
4200 )
4201}
4202
4203fn exact_noncanonical_binomial_observation_row(
4204 row: usize,
4205 y: f64,
4206 weight: f64,
4207 eta: f64,
4208 link: &InverseLink,
4209) -> Result<ExactStandardObservationRow, EstimationError> {
4210 let jet = binomial_natural_jet(row, eta, link)?;
4211 if weight == 0.0 {
4212 return Ok(ExactStandardObservationRow::zero_weight(jet.mu));
4213 }
4214 let fisherweight = weighted_positive_from_log(weight, jet.log_fisher);
4215 if !(fisherweight.is_finite() && fisherweight > 0.0) {
4216 return Err(bounded_row_error(
4217 row,
4218 "bounded binomial Fisher weight",
4219 eta,
4220 fisherweight,
4221 ));
4222 }
4223 let log_likelihood = weight * convex_combination(y, jet.log_mu[0], jet.log_one_minus_mu[0]);
4224 let score = weight * convex_combination(y, jet.log_mu[1], jet.log_one_minus_mu[1]);
4225 let neghessian_eta = -weight * convex_combination(y, jet.log_mu[2], jet.log_one_minus_mu[2]);
4226 let neghessian_eta_derivative =
4227 -weight * convex_combination(y, jet.log_mu[3], jet.log_one_minus_mu[3]);
4228 certify_bounded_row(
4229 row,
4230 eta,
4231 ExactStandardObservationRow {
4232 mu: jet.mu,
4233 score,
4234 fisherweight,
4235 neghessian_eta,
4236 neghessian_eta_derivative,
4237 log_likelihood,
4238 },
4239 )
4240}
4241
4242#[inline]
4243fn eta_exprel(rate: f64, eta: f64) -> f64 {
4244 (rate * eta).exp_m1() / rate
4245}
4246
4247fn validate_bounded_observation_inputs(
4248 likelihood: &gam_spec::GlmLikelihoodSpec,
4249 y: &Array1<f64>,
4250 weights: &Array1<f64>,
4251 eta: &Array1<f64>,
4252) -> Result<gam_spec::ResolvedLikelihoodScale, EstimationError> {
4253 let family = &likelihood.spec;
4254 if weights.len() != y.len() || eta.len() != y.len() {
4255 crate::bail_invalid_estim!(
4256 "bounded family observation size mismatch: y={}, weights={}, eta={}",
4257 y.len(),
4258 weights.len(),
4259 eta.len()
4260 );
4261 }
4262 if !LikelihoodSpec::is_legal_cell(&family.response, &family.link) {
4263 crate::bail_invalid_estim!(
4264 "bounded family received illegal likelihood cell response={} link={}",
4265 family.response.name(),
4266 family.link.link_function().name()
4267 );
4268 }
4269 let resolved_scale = likelihood
4270 .resolved_scale()
4271 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
4272 match &family.response {
4273 ResponseFamily::Tweedie { p } if !(p.is_finite() && *p > 1.0 && *p < 2.0) => {
4274 crate::bail_invalid_estim!(
4275 "bounded Tweedie power must be finite and strictly inside (1, 2), got {p}"
4276 );
4277 }
4278 ResponseFamily::NegativeBinomial { theta, .. } if !(theta.is_finite() && *theta > 0.0) => {
4279 crate::bail_invalid_estim!(
4280 "bounded negative-binomial theta must be finite and positive, got {theta}"
4281 );
4282 }
4283 _ => {}
4284 }
4285 for (i, &wi) in weights.iter().enumerate() {
4288 if !(wi.is_finite() && wi >= 0.0) {
4289 return Err(EstimationError::InvalidInput(format!(
4290 "bounded-family row {} has invalid prior weight {wi:?}; expected finite weight >= 0",
4291 i + 1
4292 )));
4293 }
4294 }
4295 for i in 0..y.len() {
4296 let wi = weights[i];
4297 if wi == 0.0 {
4298 continue;
4299 }
4300 if !eta[i].is_finite() {
4301 return Err(bounded_row_error(i, "linear predictor", eta[i], eta[i]));
4302 }
4303 if !y[i].is_finite() {
4304 return Err(bounded_row_error(
4305 i,
4306 "bounded-family response",
4307 eta[i],
4308 y[i],
4309 ));
4310 }
4311 let yi = y[i];
4312 let valid = match &family.response {
4313 ResponseFamily::Gaussian => yi.is_finite(),
4314 ResponseFamily::Binomial => yi.is_finite() && (0.0..=1.0).contains(&yi),
4315 ResponseFamily::Poisson | ResponseFamily::NegativeBinomial { .. } => {
4316 yi.is_finite() && yi >= 0.0 && (yi - yi.round()).abs() <= 1e-9
4317 }
4318 ResponseFamily::Tweedie { .. } => yi.is_finite() && yi >= 0.0,
4319 ResponseFamily::Gamma => yi.is_finite() && yi > 0.0,
4320 ResponseFamily::Beta { .. } | ResponseFamily::RoystonParmar => false,
4321 };
4322 if !valid {
4323 return Err(bounded_row_error(i, "bounded-family response", eta[i], yi));
4324 }
4325 }
4326 Ok(resolved_scale)
4327}
4328
4329fn exact_standard_observation_row(
4330 likelihood: &gam_spec::GlmLikelihoodSpec,
4331 resolved_scale: gam_spec::ResolvedLikelihoodScale,
4332 binomial_link: &InverseLink,
4333 row: usize,
4334 y: f64,
4335 weight: f64,
4336 eta: f64,
4337) -> Result<ExactStandardObservationRow, EstimationError> {
4338 if weight == 0.0 {
4339 return Ok(ExactStandardObservationRow::zero_weight(0.0));
4340 }
4341 let family = &likelihood.spec;
4342 match &family.response {
4343 ResponseFamily::Gaussian => {
4344 let scaled_weight = match resolved_scale {
4345 gam_spec::ResolvedLikelihoodScale::ProfiledGaussian => weight,
4346 gam_spec::ResolvedLikelihoodScale::FixedGaussian { phi } => {
4347 crate::gamlss::scaled_positive_product_quotient(weight, 1.0, 1.0, phi.value())
4348 }
4349 _ => {
4350 crate::bail_invalid_estim!(
4351 "bounded Gaussian received a non-Gaussian resolved scale"
4352 );
4353 }
4354 };
4355 if !(scaled_weight.is_finite() && scaled_weight > 0.0) {
4356 return Err(bounded_row_error(
4357 row,
4358 "bounded Gaussian dispersion-scaled weight",
4359 eta,
4360 scaled_weight,
4361 ));
4362 }
4363 let residual = y - eta;
4364 let loss = if residual == 0.0 {
4365 0.0
4366 } else {
4367 crate::gamlss::scaled_positive_product_quotient(
4368 scaled_weight,
4369 residual.abs(),
4370 residual.abs(),
4371 2.0,
4372 )
4373 };
4374 certify_bounded_row(
4375 row,
4376 eta,
4377 ExactStandardObservationRow {
4378 mu: eta,
4379 score: scaled_weight * residual,
4380 fisherweight: scaled_weight,
4381 neghessian_eta: scaled_weight,
4382 neghessian_eta_derivative: 0.0,
4383 log_likelihood: -loss,
4384 },
4385 )
4386 }
4387 ResponseFamily::Binomial
4388 if matches!(binomial_link, InverseLink::Standard(StandardLink::Logit)) =>
4389 {
4390 exact_logit_observation_row(row, y, weight, eta)
4391 }
4392 ResponseFamily::Binomial => {
4393 exact_noncanonical_binomial_observation_row(row, y, weight, eta, binomial_link)
4394 }
4395 ResponseFamily::Poisson => {
4396 let mu = inverse_link_jet_for_inverse_link(&family.link, eta)?.mu;
4397 let fisherweight = weight * mu;
4398 let score = weight * (y - mu);
4399 let raw_log_likelihood = y.mul_add(eta, -mu);
4400 let log_likelihood = if raw_log_likelihood.is_finite() {
4401 weight * raw_log_likelihood
4402 } else {
4403 weighted_product3(weight, y, eta) - weight * mu
4404 };
4405 if !(fisherweight.is_finite() && fisherweight > 0.0) {
4406 return Err(bounded_row_error(
4407 row,
4408 "bounded Poisson Fisher weight",
4409 eta,
4410 fisherweight,
4411 ));
4412 }
4413 certify_bounded_row(
4414 row,
4415 eta,
4416 ExactStandardObservationRow {
4417 mu,
4418 score,
4419 fisherweight,
4420 neghessian_eta: fisherweight,
4421 neghessian_eta_derivative: fisherweight,
4422 log_likelihood,
4423 },
4424 )
4425 }
4426 ResponseFamily::Gamma => {
4427 let mu = inverse_link_jet_for_inverse_link(&family.link, eta)?.mu;
4428 let shape = resolved_scale
4429 .gamma_shape()
4430 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
4431 let weighted_shape = weight * shape;
4432 if !(weighted_shape.is_finite() && weighted_shape > 0.0) {
4433 return Err(bounded_row_error(
4434 row,
4435 "bounded Gamma shape-scaled weight",
4436 eta,
4437 weighted_shape,
4438 ));
4439 }
4440 let weighted_ratio =
4441 crate::gamlss::scaled_positive_product_quotient(weight, y, shape, mu);
4442 if !(weighted_ratio.is_finite() && weighted_ratio > 0.0) {
4443 return Err(bounded_row_error(
4444 row,
4445 "bounded Gamma observed Hessian",
4446 eta,
4447 weighted_ratio,
4448 ));
4449 }
4450 certify_bounded_row(
4451 row,
4452 eta,
4453 ExactStandardObservationRow {
4454 mu,
4455 score: weighted_ratio - weighted_shape,
4456 fisherweight: weighted_shape,
4457 neghessian_eta: weighted_ratio,
4458 neghessian_eta_derivative: -weighted_ratio,
4459 log_likelihood: -weighted_ratio - weighted_shape * eta,
4460 },
4461 )
4462 }
4463 ResponseFamily::Tweedie { p } => {
4464 let p = *p;
4465 let mu = inverse_link_jet_for_inverse_link(&family.link, eta)?.mu;
4466 let phi = resolved_scale
4467 .tweedie_phi()
4468 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
4469 let weight = crate::gamlss::scaled_positive_product_quotient(weight, 1.0, 1.0, phi);
4470 if !(weight.is_finite() && weight > 0.0) {
4471 return Err(bounded_row_error(
4472 row,
4473 "bounded Tweedie dispersion-scaled weight",
4474 eta,
4475 weight,
4476 ));
4477 }
4478 let a = ((1.0 - p) * eta).exp();
4479 let b = ((2.0 - p) * eta).exp();
4480 let score_unit = y.mul_add(a, -b);
4481 let score = if score_unit.is_finite() {
4482 weight * score_unit
4483 } else {
4484 weighted_product3(weight, y, a) - weight * b
4485 };
4486 let fisherweight = weight * b;
4487 let observed_unit = (p - 1.0) * y * a + (2.0 - p) * b;
4488 let neghessian_eta = if observed_unit.is_finite() {
4489 weight * observed_unit
4490 } else {
4491 weighted_product3(weight * (p - 1.0), y, a) + weight * (2.0 - p) * b
4492 };
4493 let observed_derivative_unit = -(p - 1.0).powi(2) * y * a + (2.0 - p).powi(2) * b;
4494 let neghessian_eta_derivative = if observed_derivative_unit.is_finite() {
4495 weight * observed_derivative_unit
4496 } else {
4497 -weighted_product3(weight * (p - 1.0).powi(2), y, a)
4498 + weight * (2.0 - p).powi(2) * b
4499 };
4500 let q_left = eta_exprel(1.0 - p, eta);
4503 let q_right = eta_exprel(2.0 - p, eta);
4504 let q = y.mul_add(q_left, -q_right);
4505 let log_likelihood = if q.is_finite() {
4506 weight * q
4507 } else {
4508 weighted_product3(weight, y, q_left) - weight * q_right
4509 };
4510 if !(fisherweight.is_finite() && fisherweight > 0.0) {
4511 return Err(bounded_row_error(
4512 row,
4513 "bounded Tweedie Fisher weight",
4514 eta,
4515 fisherweight,
4516 ));
4517 }
4518 certify_bounded_row(
4519 row,
4520 eta,
4521 ExactStandardObservationRow {
4522 mu,
4523 score,
4524 fisherweight,
4525 neghessian_eta,
4526 neghessian_eta_derivative,
4527 log_likelihood,
4528 },
4529 )
4530 }
4531 ResponseFamily::NegativeBinomial { .. } => {
4532 let theta = resolved_scale
4533 .negative_binomial_theta()
4534 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
4535 let mu = inverse_link_jet_for_inverse_link(&family.link, eta)?.mu;
4536 let log_theta = theta.ln();
4537 let delta = eta - log_theta;
4538 let log_q = -gam_linalg::utils::stable_softplus(-delta);
4539 let log_r = -gam_linalg::utils::stable_softplus(delta);
4540 let q = log_q.exp();
4541 let r = log_r.exp();
4542 let y_r = if y == 0.0 {
4543 0.0
4544 } else {
4545 (y.ln() + log_r).exp()
4546 };
4547 let theta_q = (log_theta + log_q).exp();
4548 let score = weight * (y_r - theta_q);
4549 let fisherweight = weighted_positive_from_log(weight, log_theta + log_q);
4550 let log_qr = log_q + log_r;
4551 let observed_y = if y == 0.0 {
4552 0.0
4553 } else {
4554 weighted_positive_from_log(weight, y.ln() + log_qr)
4555 };
4556 let observed_theta = weighted_positive_from_log(weight, log_theta + log_qr);
4557 let neghessian_eta = observed_y + observed_theta;
4558 let neghessian_eta_derivative = neghessian_eta * (r - q);
4559 let softplus_tail = if delta >= 0.0 {
4560 gam_linalg::utils::stable_softplus(-delta)
4561 } else {
4562 gam_linalg::utils::stable_softplus(delta)
4563 };
4564 let log_likelihood = if delta >= 0.0 {
4565 -weighted_product3(weight, theta, delta)
4566 - weighted_product3(weight, y, softplus_tail)
4567 - weighted_product3(weight, theta, softplus_tail)
4568 } else {
4569 weighted_product3(weight, y, delta)
4570 - weighted_product3(weight, y, softplus_tail)
4571 - weighted_product3(weight, theta, softplus_tail)
4572 };
4573 if !(fisherweight.is_finite() && fisherweight > 0.0) {
4574 return Err(bounded_row_error(
4575 row,
4576 "bounded negative-binomial Fisher weight",
4577 eta,
4578 fisherweight,
4579 ));
4580 }
4581 certify_bounded_row(
4582 row,
4583 eta,
4584 ExactStandardObservationRow {
4585 mu,
4586 score,
4587 fisherweight,
4588 neghessian_eta,
4589 neghessian_eta_derivative,
4590 log_likelihood,
4591 },
4592 )
4593 }
4594 ResponseFamily::Beta { .. } => {
4595 crate::bail_invalid_estim!("bounded linear terms are not supported for BetaLogit fits");
4596 }
4597 ResponseFamily::RoystonParmar => {
4598 crate::bail_invalid_estim!(
4599 "bounded linear terms are not supported for survival model fits"
4600 );
4601 }
4602 }
4603}
4604
4605fn evaluate_resolved_standard_family_observations(
4606 likelihood: &gam_spec::GlmLikelihoodSpec,
4607 latent_cloglog_state: Option<&LatentCLogLogState>,
4608 mixture_link_state: Option<&MixtureLinkState>,
4609 sas_link_state: Option<&SasLinkState>,
4610 y: &Array1<f64>,
4611 weights: &Array1<f64>,
4612 eta: &Array1<f64>,
4613) -> Result<StandardFamilyObservationState, EstimationError> {
4614 let n = y.len();
4615 let resolved_scale = validate_bounded_observation_inputs(likelihood, y, weights, eta)?;
4616 let family = &likelihood.spec;
4617 let binomial_link = resolved_bounded_binomial_link(
4618 &family,
4619 latent_cloglog_state,
4620 mixture_link_state,
4621 sas_link_state,
4622 );
4623
4624 let mut score = Array1::<f64>::zeros(n);
4625 let mut fisherweight = Array1::<f64>::zeros(n);
4626 let mut neghessian_eta = Array1::<f64>::zeros(n);
4627 let mut neghessian_eta_derivative = Array1::<f64>::zeros(n);
4628 let mut log_likelihood = 0.0;
4629 let mut log_likelihood_compensation = 0.0;
4630
4631 for i in 0..n {
4632 let row = exact_standard_observation_row(
4633 likelihood,
4634 resolved_scale,
4635 &binomial_link,
4636 i,
4637 y[i],
4638 weights[i],
4639 eta[i],
4640 )?;
4641 score[i] = row.score;
4642 fisherweight[i] = row.fisherweight;
4643 neghessian_eta[i] = row.neghessian_eta;
4644 neghessian_eta_derivative[i] = row.neghessian_eta_derivative;
4645 let adjusted = row.log_likelihood - log_likelihood_compensation;
4646 let updated = log_likelihood + adjusted;
4647 log_likelihood_compensation = (updated - log_likelihood) - adjusted;
4648 log_likelihood = updated;
4649 if !log_likelihood.is_finite() {
4650 return Err(bounded_row_error(
4651 i,
4652 "bounded-family cumulative log likelihood",
4653 eta[i],
4654 log_likelihood,
4655 ));
4656 }
4657 }
4658
4659 Ok(StandardFamilyObservationState {
4660 eta: eta.clone(),
4661 score,
4662 fisherweight,
4663 neghessian_eta,
4664 neghessian_eta_derivative,
4665 log_likelihood,
4666 })
4667}
4668
4669fn evaluate_standard_familyobservations(
4674 family: LikelihoodSpec,
4675 latent_cloglog_state: Option<&LatentCLogLogState>,
4676 mixture_link_state: Option<&MixtureLinkState>,
4677 sas_link_state: Option<&SasLinkState>,
4678 y: &Array1<f64>,
4679 weights: &Array1<f64>,
4680 eta: &Array1<f64>,
4681) -> Result<StandardFamilyObservationState, EstimationError> {
4682 let likelihood = gam_spec::GlmLikelihoodSpec::canonical(family);
4683 evaluate_resolved_standard_family_observations(
4684 &likelihood,
4685 latent_cloglog_state,
4686 mixture_link_state,
4687 sas_link_state,
4688 y,
4689 weights,
4690 eta,
4691 )
4692}
4693
4694fn exact_standard_working_response(
4695 state: &StandardFamilyObservationState,
4696) -> Result<Array1<f64>, EstimationError> {
4697 let mut out = state.eta.clone();
4698 for i in 0..out.len() {
4699 let weight = state.fisherweight[i];
4700 let score = state.score[i];
4701 if weight == 0.0 {
4702 if score != 0.0 {
4703 return Err(bounded_row_error(
4704 i,
4705 "zero-Fisher row with nonzero score",
4706 state.eta[i],
4707 score,
4708 ));
4709 }
4710 continue;
4711 }
4712 let increment = score / weight;
4713 let value = out[i] + increment;
4714 if !increment.is_finite() || !value.is_finite() {
4715 return Err(bounded_row_error(
4716 i,
4717 "bounded-family working response",
4718 state.eta[i],
4719 value,
4720 ));
4721 }
4722 out[i] = value;
4723 }
4724 Ok(out)
4725}
4726
4727#[derive(Clone, Copy, Debug, PartialEq, Eq)]
4728enum SpatialAdaptiveHyperKind {
4729 LogLambdaMagnitude,
4730 LogLambdaGradient,
4731 LogLambdaCurvature,
4732 LogEpsilonMagnitude,
4733 LogEpsilonGradient,
4734 LogEpsilonCurvature,
4735}
4736
4737impl SpatialAdaptiveHyperKind {
4738 fn component_index(self) -> usize {
4739 match self {
4740 SpatialAdaptiveHyperKind::LogLambdaMagnitude
4741 | SpatialAdaptiveHyperKind::LogEpsilonMagnitude => 0,
4742 SpatialAdaptiveHyperKind::LogLambdaGradient
4743 | SpatialAdaptiveHyperKind::LogEpsilonGradient => 1,
4744 SpatialAdaptiveHyperKind::LogLambdaCurvature
4745 | SpatialAdaptiveHyperKind::LogEpsilonCurvature => 2,
4746 }
4747 }
4748
4749 fn is_log_lambda(self) -> bool {
4750 matches!(
4751 self,
4752 SpatialAdaptiveHyperKind::LogLambdaMagnitude
4753 | SpatialAdaptiveHyperKind::LogLambdaGradient
4754 | SpatialAdaptiveHyperKind::LogLambdaCurvature
4755 )
4756 }
4757
4758 fn is_log_epsilon(self) -> bool {
4759 matches!(
4760 self,
4761 SpatialAdaptiveHyperKind::LogEpsilonMagnitude
4762 | SpatialAdaptiveHyperKind::LogEpsilonGradient
4763 | SpatialAdaptiveHyperKind::LogEpsilonCurvature
4764 )
4765 }
4766}
4767
4768#[derive(Clone, Copy, Debug)]
4769struct SpatialAdaptiveHyperSpec {
4770 cache_index: usize,
4771 kind: SpatialAdaptiveHyperKind,
4772}
4773
4774#[derive(Clone, Copy, Debug, PartialEq, Eq)]
4775enum SpatialAdaptiveExplicitSecondOrderKind {
4776 StructuralZero,
4777 LocalAlphaAlpha,
4778 LocalAlphaEta,
4779 SharedEtaEta,
4780}
4781
4782#[derive(Clone, Copy, Debug, PartialEq, Eq)]
4787enum AdaptiveComponent {
4788 Magnitude,
4789 Gradient,
4790 Curvature,
4791}
4792
4793impl AdaptiveComponent {
4794 fn from_index(index: usize) -> Result<Self, String> {
4795 match index {
4796 0 => Ok(AdaptiveComponent::Magnitude),
4797 1 => Ok(AdaptiveComponent::Gradient),
4798 2 => Ok(AdaptiveComponent::Curvature),
4799 other => Err(SmoothError::invalid_index(format!(
4800 "invalid adaptive component index {}",
4801 other
4802 ))
4803 .into()),
4804 }
4805 }
4806}
4807
4808#[derive(Clone, Copy, Debug, PartialEq, Eq)]
4813enum HyperDerivativeKind {
4814 Rho,
4816 LogEpsilonFirst,
4818 LogEpsilonSecond,
4820}
4821
4822#[derive(Clone, Copy, Debug, PartialEq, Eq)]
4827enum HyperDriftKind {
4828 Rho,
4829 LogEpsilon,
4830}
4831
4832impl SpatialAdaptiveHyperSpec {
4833 fn component_index(self) -> usize {
4834 self.kind.component_index()
4835 }
4836
4837 fn explicit_second_order_kind(self, other: Self) -> SpatialAdaptiveExplicitSecondOrderKind {
4838 if self.component_index() != other.component_index() {
4839 return SpatialAdaptiveExplicitSecondOrderKind::StructuralZero;
4840 }
4841 match (
4842 self.kind.is_log_lambda(),
4843 other.kind.is_log_lambda(),
4844 self.kind.is_log_epsilon(),
4845 other.kind.is_log_epsilon(),
4846 ) {
4847 (true, true, false, false) if self.cache_index == other.cache_index => {
4848 SpatialAdaptiveExplicitSecondOrderKind::LocalAlphaAlpha
4849 }
4850 (true, false, false, true) | (false, true, true, false) => {
4851 SpatialAdaptiveExplicitSecondOrderKind::LocalAlphaEta
4852 }
4853 (false, false, true, true) => SpatialAdaptiveExplicitSecondOrderKind::SharedEtaEta,
4854 _ => SpatialAdaptiveExplicitSecondOrderKind::StructuralZero,
4855 }
4856 }
4857}
4858
4859#[derive(Clone, Debug)]
4860struct SpatialAdaptiveTermHyperParams {
4861 lambda: [f64; 3],
4862 epsilon: [f64; 3],
4863}
4864
4865#[derive(Clone, Debug)]
4874struct ValidatedFixedQuadraticHessian {
4875 dense: Arc<Array2<f64>>,
4876}
4877
4878impl ValidatedFixedQuadraticHessian {
4879 fn try_from_dense(dense: Array2<f64>, coefficient_dim: usize) -> Result<Self, String> {
4880 gam_linalg::utils::validate_finite_symmetric_matrix(
4881 &dense,
4882 "spatial adaptive fixed quadratic Hessian",
4883 )
4884 .map_err(|error| error.to_string())?;
4885 PenaltyMatrix::Dense(dense.clone())
4886 .validate(coefficient_dim)
4887 .map_err(|error| {
4888 format!(
4889 "spatial adaptive fixed quadratic Hessian failed quadratic-form validation: {error}"
4890 )
4891 })?;
4892 Ok(Self {
4893 dense: Arc::new(dense),
4894 })
4895 }
4896
4897 fn zero(coefficient_dim: usize) -> Result<Self, String> {
4898 Self::try_from_dense(
4899 Array2::<f64>::zeros((coefficient_dim, coefficient_dim)),
4900 coefficient_dim,
4901 )
4902 }
4903
4904 fn as_dense(&self) -> &Array2<f64> {
4905 self.dense.as_ref()
4906 }
4907
4908 fn quadratic_terms(&self, beta: &Array1<f64>) -> Result<(f64, Array1<f64>), String> {
4909 if beta.len() != self.dense.ncols() {
4910 return Err(format!(
4911 "spatial adaptive fixed quadratic beta length {} does not match validated Hessian dimension {}",
4912 beta.len(),
4913 self.dense.ncols()
4914 ));
4915 }
4916 let gradient = self.dense.dot(beta);
4917 let value = 0.5 * beta.dot(&gradient);
4918 Ok((value, gradient))
4919 }
4920}
4921
4922#[derive(Clone)]
4923struct SpatialAdaptiveExactEvaluation {
4924 obs: StandardFamilyObservationState,
4925 adaptive_states: Vec<SpatialPenaltyExactState>,
4926 adaptive_penalty_value: f64,
4927 adaptive_penaltygradient: Array1<f64>,
4928 adaptive_penaltyhessian: Array2<f64>,
4929 fixed_quadraticvalue: f64,
4930 fixed_quadraticgradient: Array1<f64>,
4931 fixed_quadratic_hessian: ValidatedFixedQuadraticHessian,
4932}
4933
4934#[derive(Clone)]
4935struct CachedSpatialAdaptiveExactEvaluation {
4936 beta: Array1<f64>,
4937 eval: Arc<SpatialAdaptiveExactEvaluation>,
4938}
4939
4940impl SpatialAdaptiveExactEvaluation {
4941 fn total_penalty_value(&self) -> f64 {
4942 self.adaptive_penalty_value + self.fixed_quadraticvalue
4943 }
4944
4945 fn total_penaltygradient(&self) -> Array1<f64> {
4946 &self.adaptive_penaltygradient + &self.fixed_quadraticgradient
4947 }
4948
4949 fn total_penaltyhessian(&self) -> Array2<f64> {
4950 &self.adaptive_penaltyhessian + self.fixed_quadratic_hessian.as_dense()
4951 }
4952
4953 fn totalobjectivehessian(&self, design: &Array2<f64>) -> Result<Array2<f64>, String> {
4954 let mut out = xt_diag_x_dense(design.view(), self.obs.neghessian_eta.view())?;
4955 out += &self.total_penaltyhessian();
4956 Ok(out)
4957 }
4958}
4959
4960#[derive(Clone)]
4961struct SpatialAdaptiveExactFamily {
4962 family: LikelihoodSpec,
4963 latent_cloglog_state: Option<LatentCLogLogState>,
4964 mixture_link_state: Option<MixtureLinkState>,
4965 sas_link_state: Option<SasLinkState>,
4966 y: Arc<Array1<f64>>,
4967 weights: Arc<Array1<f64>>,
4968 design: Arc<Array2<f64>>,
4969 offset: Arc<Array1<f64>>,
4970 linear_constraints: Option<LinearInequalityConstraints>,
4971 runtime_caches: Arc<Vec<SpatialOperatorRuntimeCache>>,
4972 adaptive_params: Vec<SpatialAdaptiveTermHyperParams>,
4973 fixed_quadratic_hessian: ValidatedFixedQuadraticHessian,
4974 hyperspecs: Arc<Vec<SpatialAdaptiveHyperSpec>>,
4975 exact_eval_cache: Arc<Mutex<Option<CachedSpatialAdaptiveExactEvaluation>>>,
4976}
4977
4978impl SpatialAdaptiveExactFamily {
4979 fn with_adaptive_params(
4980 &self,
4981 adaptive_params: Vec<SpatialAdaptiveTermHyperParams>,
4982 fixed_quadratic_hessian: ValidatedFixedQuadraticHessian,
4983 ) -> Self {
4984 Self {
4985 family: self.family.clone(),
4986 latent_cloglog_state: self.latent_cloglog_state,
4987 mixture_link_state: self.mixture_link_state.clone(),
4988 sas_link_state: self.sas_link_state,
4989 y: self.y.clone(),
4990 weights: self.weights.clone(),
4991 design: self.design.clone(),
4992 offset: self.offset.clone(),
4993 linear_constraints: self.linear_constraints.clone(),
4994 runtime_caches: self.runtime_caches.clone(),
4995 adaptive_params,
4996 fixed_quadratic_hessian,
4997 hyperspecs: self.hyperspecs.clone(),
4998 exact_eval_cache: Arc::new(Mutex::new(None)),
4999 }
5000 }
5001
5002 fn total_eta(&self, beta: &Array1<f64>) -> Array1<f64> {
5003 gam_linalg::faer_ndarray::fast_av(self.design.as_ref(), beta) + self.offset.as_ref()
5004 }
5005
5006 fn fixed_quadratic_terms(
5007 &self,
5008 beta: &Array1<f64>,
5009 ) -> Result<(f64, Array1<f64>), String> {
5010 self.fixed_quadratic_hessian.quadratic_terms(beta)
5011 }
5012
5013 fn adaptive_penalty_value_only(&self, beta: &Array1<f64>) -> Result<f64, String> {
5014 let mut penalty_value = 0.0;
5015 for (cache_idx, cache) in self.runtime_caches.iter().enumerate() {
5016 let params = self.adaptive_params.get(cache_idx).ok_or_else(|| {
5017 format!(
5018 "missing adaptive parameter block for cache {}",
5019 cache.termname
5020 )
5021 })?;
5022 let beta_local = beta.slice(s![cache.coeff_global_range.clone()]);
5023 let state =
5024 SpatialPenaltyExactState::from_beta_local(beta_local, cache, params.epsilon)
5025 .map_err(|e| e.to_string())?;
5026 penalty_value += params.lambda[0] * state.magnitude.penalty_value();
5027 penalty_value += params.lambda[1] * state.gradient.penalty_value();
5028 penalty_value += params.lambda[2] * state.curvature.penalty_value();
5029 }
5030 Ok(penalty_value)
5031 }
5032
5033 fn zero_hyper_parts(&self) -> (Array1<f64>, Array2<f64>) {
5034 let total_dim = self.design.ncols();
5035 (
5036 Array1::<f64>::zeros(total_dim),
5037 Array2::<f64>::zeros((total_dim, total_dim)),
5038 )
5039 }
5040
5041 fn embed_local_hyper_parts(
5042 &self,
5043 coeff_range: &Range<usize>,
5044 local_grad: &Array1<f64>,
5045 local_hess: &Array2<f64>,
5046 ) -> (Array1<f64>, Array2<f64>) {
5047 let (mut beta_mixed, mut betahessian) = self.zero_hyper_parts();
5048 beta_mixed
5049 .slice_mut(s![coeff_range.clone()])
5050 .assign(local_grad);
5051 betahessian
5052 .slice_mut(s![coeff_range.clone(), coeff_range.clone()])
5053 .assign(local_hess);
5054 (beta_mixed, betahessian)
5055 }
5056
5057 fn embed_local_hyper_hessian(
5058 &self,
5059 coeff_range: &Range<usize>,
5060 local_hess: &Array2<f64>,
5061 ) -> Array2<f64> {
5062 let total_dim = self.design.ncols();
5063 let mut out = Array2::<f64>::zeros((total_dim, total_dim));
5064 out.slice_mut(s![coeff_range.clone(), coeff_range.clone()])
5065 .assign(local_hess);
5066 out
5067 }
5068
5069 fn adaptive_block_eval(
5078 &self,
5079 eval: &SpatialAdaptiveExactEvaluation,
5080 cache_idx: usize,
5081 component: AdaptiveComponent,
5082 derivative: HyperDerivativeKind,
5083 ) -> Result<(f64, Array1<f64>, Array2<f64>), String> {
5084 let cache = self
5085 .runtime_caches
5086 .get(cache_idx)
5087 .ok_or_else(|| format!("adaptive cache index {} out of bounds", cache_idx))?;
5088 let params = self
5089 .adaptive_params
5090 .get(cache_idx)
5091 .ok_or_else(|| format!("adaptive hyperparameter block {} out of bounds", cache_idx))?;
5092 let state = eval
5093 .adaptive_states
5094 .get(cache_idx)
5095 .ok_or_else(|| format!("adaptive exact state index {} out of bounds", cache_idx))?;
5096
5097 let (objective_local, beta_mixed_local, betahessian_local) = match component {
5098 AdaptiveComponent::Magnitude => {
5099 let lambda = params.lambda[0];
5100 let mag = &state.magnitude;
5101 let (objective, gradient_coeff, hessian_diag) = match derivative {
5102 HyperDerivativeKind::Rho => (
5103 mag.penalty_value(),
5104 mag.betagradient_coeff(),
5105 mag.betahessian_diag(),
5106 ),
5107 HyperDerivativeKind::LogEpsilonFirst => (
5108 mag.log_epsilon_gradient_terms().sum(),
5109 mag.log_epsilon_betagradient_coeff(),
5110 mag.log_epsilon_betahessian_diag(),
5111 ),
5112 HyperDerivativeKind::LogEpsilonSecond => (
5113 mag.log_epsilon_hessian_terms().sum(),
5114 mag.log_epsilon_beta_mixed_second_coeff(),
5115 mag.log_epsilon_betahessian_second_diag(),
5116 ),
5117 };
5118 (
5119 lambda * objective,
5120 lambda * scalar_operatorgradient(&cache.d0, &gradient_coeff),
5121 lambda * scalar_operatorhessian(&cache.d0, &hessian_diag),
5122 )
5123 }
5124 AdaptiveComponent::Gradient => {
5125 let lambda = params.lambda[1];
5126 let grad = &state.gradient;
5127 let (objective, gradient_blocks, hessian_blocks) = match derivative {
5128 HyperDerivativeKind::Rho => (
5129 grad.penalty_value(),
5130 grad.betagradient_blocks(),
5131 grad.betahessian_blocks(),
5132 ),
5133 HyperDerivativeKind::LogEpsilonFirst => (
5134 grad.log_epsilon_gradient_terms().sum(),
5135 grad.log_epsilon_betagradient_blocks(),
5136 grad.log_epsilon_betahessian_blocks(),
5137 ),
5138 HyperDerivativeKind::LogEpsilonSecond => (
5139 grad.log_epsilon_hessian_terms().sum(),
5140 grad.log_epsilon_beta_mixed_second_blocks(),
5141 grad.log_epsilon_betahessian_second_blocks(),
5142 ),
5143 };
5144 (
5145 lambda * objective,
5146 lambda
5147 * grouped_operatorgradient(&cache.d1, cache.dimension, &gradient_blocks)
5148 .map_err(|e| e.to_string())?,
5149 lambda
5150 * grouped_operatorhessian(&cache.d1, cache.dimension, &hessian_blocks)
5151 .map_err(|e| e.to_string())?,
5152 )
5153 }
5154 AdaptiveComponent::Curvature => {
5155 let lambda = params.lambda[2];
5156 let group = cache.dimension * cache.dimension;
5157 let curv = &state.curvature;
5158 let (objective, gradient_blocks, hessian_blocks) = match derivative {
5159 HyperDerivativeKind::Rho => (
5160 curv.penalty_value(),
5161 curv.betagradient_blocks(),
5162 curv.betahessian_blocks(),
5163 ),
5164 HyperDerivativeKind::LogEpsilonFirst => (
5165 curv.log_epsilon_gradient_terms().sum(),
5166 curv.log_epsilon_betagradient_blocks(),
5167 curv.log_epsilon_betahessian_blocks(),
5168 ),
5169 HyperDerivativeKind::LogEpsilonSecond => (
5170 curv.log_epsilon_hessian_terms().sum(),
5171 curv.log_epsilon_beta_mixed_second_blocks(),
5172 curv.log_epsilon_betahessian_second_blocks(),
5173 ),
5174 };
5175 (
5176 lambda * objective,
5177 lambda
5178 * grouped_operatorgradient(&cache.d2, group, &gradient_blocks)
5179 .map_err(|e| e.to_string())?,
5180 lambda
5181 * grouped_operatorhessian(&cache.d2, group, &hessian_blocks)
5182 .map_err(|e| e.to_string())?,
5183 )
5184 }
5185 };
5186
5187 let (beta_mixed, betahessian) = self.embed_local_hyper_parts(
5188 &cache.coeff_global_range,
5189 &beta_mixed_local,
5190 &betahessian_local,
5191 );
5192 Ok((objective_local, beta_mixed, betahessian))
5193 }
5194
5195 fn adaptive_shared_log_epsilon_parts(
5196 &self,
5197 eval: &SpatialAdaptiveExactEvaluation,
5198 component: usize,
5199 ) -> Result<(f64, Array1<f64>, Array2<f64>), String> {
5200 self.adaptive_shared_block_eval(eval, component, HyperDerivativeKind::LogEpsilonFirst)
5206 }
5207
5208 fn adaptive_shared_log_epsilon_second_parts(
5209 &self,
5210 eval: &SpatialAdaptiveExactEvaluation,
5211 component: usize,
5212 ) -> Result<(f64, Array1<f64>, Array2<f64>), String> {
5213 self.adaptive_shared_block_eval(eval, component, HyperDerivativeKind::LogEpsilonSecond)
5219 }
5220
5221 fn adaptive_shared_block_eval(
5226 &self,
5227 eval: &SpatialAdaptiveExactEvaluation,
5228 component: usize,
5229 derivative: HyperDerivativeKind,
5230 ) -> Result<(f64, Array1<f64>, Array2<f64>), String> {
5231 let component = AdaptiveComponent::from_index(component)?;
5232 let (mut score, mut hessian) = self.zero_hyper_parts();
5233 let mut objective = 0.0;
5234 for cache_idx in 0..self.runtime_caches.len() {
5235 let (local_objective, local_score, local_hessian) =
5236 self.adaptive_block_eval(eval, cache_idx, component, derivative)?;
5237 objective += local_objective;
5238 score += &local_score;
5239 hessian += &local_hessian;
5240 }
5241 Ok((objective, score, hessian))
5242 }
5243
5244 fn adaptive_shared_log_epsilon_drift(
5245 &self,
5246 eval: &SpatialAdaptiveExactEvaluation,
5247 component: usize,
5248 direction: &Array1<f64>,
5249 ) -> Result<Array2<f64>, String> {
5250 let component = AdaptiveComponent::from_index(component)?;
5254 let total_dim = self.design.ncols();
5255 let mut total = Array2::<f64>::zeros((total_dim, total_dim));
5256 for cache_idx in 0..self.runtime_caches.len() {
5257 total += &self.adaptive_block_drift_eval(
5258 eval,
5259 cache_idx,
5260 component,
5261 HyperDriftKind::LogEpsilon,
5262 direction,
5263 )?;
5264 }
5265 Ok(total)
5266 }
5267
5268 fn adaptive_explicit_second_order_parts(
5269 &self,
5270 eval: &SpatialAdaptiveExactEvaluation,
5271 left: SpatialAdaptiveHyperSpec,
5272 right: SpatialAdaptiveHyperSpec,
5273 ) -> Result<(f64, Array1<f64>, Array2<f64>), String> {
5274 match left.explicit_second_order_kind(right) {
5283 SpatialAdaptiveExplicitSecondOrderKind::StructuralZero => {
5284 let (score, hessian) = self.zero_hyper_parts();
5285 Ok((0.0, score, hessian))
5286 }
5287 SpatialAdaptiveExplicitSecondOrderKind::LocalAlphaAlpha => self.adaptive_block_eval(
5288 eval,
5289 left.cache_index,
5290 AdaptiveComponent::from_index(left.component_index())?,
5291 HyperDerivativeKind::Rho,
5292 ),
5293 SpatialAdaptiveExplicitSecondOrderKind::LocalAlphaEta => {
5294 let local_alpha = if left.kind.is_log_lambda() {
5295 left
5296 } else {
5297 right
5298 };
5299 self.adaptive_block_eval(
5300 eval,
5301 local_alpha.cache_index,
5302 AdaptiveComponent::from_index(local_alpha.component_index())?,
5303 HyperDerivativeKind::LogEpsilonFirst,
5304 )
5305 }
5306 SpatialAdaptiveExplicitSecondOrderKind::SharedEtaEta => {
5307 self.adaptive_shared_log_epsilon_second_parts(eval, left.component_index())
5308 }
5309 }
5310 }
5311
5312 fn adaptive_block_drift_eval(
5320 &self,
5321 eval: &SpatialAdaptiveExactEvaluation,
5322 cache_idx: usize,
5323 component: AdaptiveComponent,
5324 drift: HyperDriftKind,
5325 direction: &Array1<f64>,
5326 ) -> Result<Array2<f64>, String> {
5327 let cache = self
5328 .runtime_caches
5329 .get(cache_idx)
5330 .ok_or_else(|| format!("adaptive cache index {} out of bounds", cache_idx))?;
5331 let params = self
5332 .adaptive_params
5333 .get(cache_idx)
5334 .ok_or_else(|| format!("adaptive hyperparameter block {} out of bounds", cache_idx))?;
5335 let state = eval
5336 .adaptive_states
5337 .get(cache_idx)
5338 .ok_or_else(|| format!("adaptive exact state index {} out of bounds", cache_idx))?;
5339 let direction_local = direction.slice(s![cache.coeff_global_range.clone()]);
5340
5341 let local_hessian = match component {
5342 AdaptiveComponent::Magnitude => {
5343 let d0_u = cache.d0.dot(&direction_local);
5344 let mag = &state.magnitude;
5345 let diag = match drift {
5346 HyperDriftKind::Rho => mag.directionalhessian_diag(&d0_u),
5347 HyperDriftKind::LogEpsilon => {
5348 mag.log_epsilon_betahessian_directional_diag(&d0_u)
5349 }
5350 };
5351 params.lambda[0] * scalar_operatorhessian(&cache.d0, &diag)
5352 }
5353 AdaptiveComponent::Gradient => {
5354 let d1_u = cache.d1.dot(&direction_local);
5355 let direction_blocks = collocationgradient_blocks(&d1_u, cache.dimension)
5356 .map_err(|e| e.to_string())?;
5357 let grad = &state.gradient;
5358 let blocks = match drift {
5359 HyperDriftKind::Rho => grad.directionalhessian_blocks(&direction_blocks),
5360 HyperDriftKind::LogEpsilon => {
5361 grad.log_epsilon_betahessian_directional_blocks(&direction_blocks)
5362 }
5363 };
5364 params.lambda[1]
5365 * grouped_operatorhessian(&cache.d1, cache.dimension, &blocks)
5366 .map_err(|e| e.to_string())?
5367 }
5368 AdaptiveComponent::Curvature => {
5369 let group = cache.dimension * cache.dimension;
5370 let d2_u = cache.d2.dot(&direction_local);
5371 let direction_blocks =
5372 collocationhessian_blocks(&d2_u, cache.dimension).map_err(|e| e.to_string())?;
5373 let curv = &state.curvature;
5374 let blocks = match drift {
5375 HyperDriftKind::Rho => curv.directionalhessian_blocks(&direction_blocks),
5376 HyperDriftKind::LogEpsilon => {
5377 curv.log_epsilon_betahessian_directional_blocks(&direction_blocks)
5378 }
5379 };
5380 params.lambda[2]
5381 * grouped_operatorhessian(&cache.d2, group, &blocks)
5382 .map_err(|e| e.to_string())?
5383 }
5384 };
5385
5386 Ok(self.embed_local_hyper_hessian(&cache.coeff_global_range, &local_hessian))
5387 }
5388
5389 fn adaptive_hyper_parts(
5390 &self,
5391 eval: &SpatialAdaptiveExactEvaluation,
5392 hyper: SpatialAdaptiveHyperSpec,
5393 ) -> Result<(f64, Array1<f64>, Array2<f64>), String> {
5394 match hyper.kind {
5395 SpatialAdaptiveHyperKind::LogLambdaMagnitude
5398 | SpatialAdaptiveHyperKind::LogLambdaGradient
5399 | SpatialAdaptiveHyperKind::LogLambdaCurvature => self.adaptive_block_eval(
5400 eval,
5401 hyper.cache_index,
5402 AdaptiveComponent::from_index(hyper.component_index())?,
5403 HyperDerivativeKind::Rho,
5404 ),
5405 SpatialAdaptiveHyperKind::LogEpsilonMagnitude
5407 | SpatialAdaptiveHyperKind::LogEpsilonGradient
5408 | SpatialAdaptiveHyperKind::LogEpsilonCurvature => {
5409 self.adaptive_shared_log_epsilon_parts(eval, hyper.component_index())
5410 }
5411 }
5412 }
5413
5414 fn exact_evaluation_uncached(
5415 &self,
5416 beta: &Array1<f64>,
5417 ) -> Result<SpatialAdaptiveExactEvaluation, String> {
5418 let eta = self.total_eta(beta);
5419 let obs = evaluate_standard_familyobservations(
5420 self.family.clone(),
5421 self.latent_cloglog_state.as_ref(),
5422 self.mixture_link_state.as_ref(),
5423 self.sas_link_state.as_ref(),
5424 &self.y,
5425 &self.weights,
5426 &eta,
5427 )
5428 .map_err(|e| e.to_string())?;
5429 let p = beta.len();
5430 let mut penalty_value = 0.0;
5431 let mut penaltygradient = Array1::<f64>::zeros(p);
5432 let mut penaltyhessian = Array2::<f64>::zeros((p, p));
5433 let mut adaptive_states = Vec::with_capacity(self.runtime_caches.len());
5434
5435 for (cache_idx, cache) in self.runtime_caches.iter().enumerate() {
5436 let params = self.adaptive_params.get(cache_idx).ok_or_else(|| {
5437 format!(
5438 "missing adaptive parameter block for cache {}",
5439 cache.termname
5440 )
5441 })?;
5442 let beta_local = beta.slice(s![cache.coeff_global_range.clone()]);
5443 let state =
5444 SpatialPenaltyExactState::from_beta_local(beta_local, cache, params.epsilon)
5445 .map_err(|e| e.to_string())?;
5446
5447 let g0 = scalar_operatorgradient(&cache.d0, &state.magnitude.betagradient_coeff());
5448 let gg = grouped_operatorgradient(
5449 &cache.d1,
5450 cache.dimension,
5451 &state.gradient.betagradient_blocks(),
5452 )
5453 .map_err(|e| e.to_string())?;
5454 let gc = grouped_operatorgradient(
5455 &cache.d2,
5456 cache.dimension * cache.dimension,
5457 &state.curvature.betagradient_blocks(),
5458 )
5459 .map_err(|e| e.to_string())?;
5460 let h0 = scalar_operatorhessian(&cache.d0, &state.magnitude.betahessian_diag());
5461 let hg = grouped_operatorhessian(
5462 &cache.d1,
5463 cache.dimension,
5464 &state.gradient.betahessian_blocks(),
5465 )
5466 .map_err(|e| e.to_string())?;
5467 let hc = grouped_operatorhessian(
5468 &cache.d2,
5469 cache.dimension * cache.dimension,
5470 &state.curvature.betahessian_blocks(),
5471 )
5472 .map_err(|e| e.to_string())?;
5473
5474 let lambda0 = params.lambda[0];
5475 let lambdag = params.lambda[1];
5476 let lambdac = params.lambda[2];
5477
5478 penalty_value += lambda0 * state.magnitude.penalty_value();
5479 penalty_value += lambdag * state.gradient.penalty_value();
5480 penalty_value += lambdac * state.curvature.penalty_value();
5481
5482 let range = cache.coeff_global_range.clone();
5483 {
5484 let mut grad_local = penaltygradient.slice_mut(s![range.clone()]);
5485 grad_local += &(g0.mapv(|v| lambda0 * v));
5486 grad_local += &(gg.mapv(|v| lambdag * v));
5487 grad_local += &(gc.mapv(|v| lambdac * v));
5488 }
5489 {
5490 let mut h_local = penaltyhessian.slice_mut(s![range.clone(), range]);
5491 h_local += &h0.mapv(|v| lambda0 * v);
5492 h_local += &hg.mapv(|v| lambdag * v);
5493 h_local += &hc.mapv(|v| lambdac * v);
5494 }
5495
5496 adaptive_states.push(state);
5497 }
5498
5499 let (fixed_quadraticvalue, fixed_quadraticgradient) =
5500 self.fixed_quadratic_terms(beta)?;
5501 Ok(SpatialAdaptiveExactEvaluation {
5502 obs,
5503 adaptive_states,
5504 adaptive_penalty_value: penalty_value,
5505 adaptive_penaltygradient: penaltygradient,
5506 adaptive_penaltyhessian: penaltyhessian,
5507 fixed_quadraticvalue,
5508 fixed_quadraticgradient,
5509 fixed_quadratic_hessian: self.fixed_quadratic_hessian.clone(),
5510 })
5511 }
5512
5513 fn exact_evaluation(
5514 &self,
5515 beta: &Array1<f64>,
5516 ) -> Result<Arc<SpatialAdaptiveExactEvaluation>, String> {
5517 {
5518 let cache = self
5519 .exact_eval_cache
5520 .lock()
5521 .map_err(|_| "spatial adaptive exact-evaluation cache lock poisoned".to_string())?;
5522 if let Some(cached) = cache.as_ref()
5523 && cached.beta.len() == beta.len()
5524 && cached
5525 .beta
5526 .iter()
5527 .zip(beta.iter())
5528 .all(|(&left, &right)| left == right)
5529 {
5530 return Ok(Arc::clone(&cached.eval));
5531 }
5532 }
5533
5534 let eval = Arc::new(self.exact_evaluation_uncached(beta)?);
5535 let mut cache = self
5536 .exact_eval_cache
5537 .lock()
5538 .map_err(|_| "spatial adaptive exact-evaluation cache lock poisoned".to_string())?;
5539 *cache = Some(CachedSpatialAdaptiveExactEvaluation {
5540 beta: beta.clone(),
5541 eval: Arc::clone(&eval),
5542 });
5543 Ok(eval)
5544 }
5545
5546 fn exacthessian_directional_derivative_from_evaluation(
5547 &self,
5548 beta: &Array1<f64>,
5549 eval: &SpatialAdaptiveExactEvaluation,
5550 direction: &Array1<f64>,
5551 ) -> Result<Array2<f64>, String> {
5552 assert_eq!(
5553 beta.len(),
5554 direction.len(),
5555 "beta/direction length mismatch",
5556 );
5557 let d_eta = gam_linalg::faer_ndarray::fast_av(self.design.as_ref(), direction);
5558 let mut total = xt_diag_x_dense(
5559 self.design.view(),
5560 (&eval.obs.neghessian_eta_derivative * &d_eta).view(),
5561 )?;
5562 for (cache_idx, cache) in self.runtime_caches.iter().enumerate() {
5563 let params = self.adaptive_params.get(cache_idx).ok_or_else(|| {
5564 format!(
5565 "missing adaptive parameter block for cache {}",
5566 cache.termname
5567 )
5568 })?;
5569 let state = eval
5570 .adaptive_states
5571 .get(cache_idx)
5572 .ok_or_else(|| format!("missing adaptive state for cache {}", cache.termname))?;
5573 let direction_local = direction.slice(s![cache.coeff_global_range.clone()]);
5574 let d0_u = cache.d0.dot(&direction_local);
5575 let d1_u = cache.d1.dot(&direction_local);
5576 let d2_u = cache.d2.dot(&direction_local);
5577 let h0 =
5578 scalar_operatorhessian(&cache.d0, &state.magnitude.directionalhessian_diag(&d0_u))
5579 .mapv(|v| params.lambda[0] * v);
5580 let hg = grouped_operatorhessian(
5581 &cache.d1,
5582 cache.dimension,
5583 &state.gradient.directionalhessian_blocks(
5584 &collocationgradient_blocks(&d1_u, cache.dimension)
5585 .map_err(|e| e.to_string())?,
5586 ),
5587 )
5588 .map_err(|e| e.to_string())?
5589 .mapv(|v| params.lambda[1] * v);
5590 let hc = grouped_operatorhessian(
5591 &cache.d2,
5592 cache.dimension * cache.dimension,
5593 &state.curvature.directionalhessian_blocks(
5594 &collocationhessian_blocks(&d2_u, cache.dimension)
5595 .map_err(|e| e.to_string())?,
5596 ),
5597 )
5598 .map_err(|e| e.to_string())?
5599 .mapv(|v| params.lambda[2] * v);
5600 let range = cache.coeff_global_range.clone();
5601 let mut local = total.slice_mut(s![range.clone(), range]);
5602 local += &h0;
5603 local += &hg;
5604 local += &hc;
5605 }
5606 Ok(total)
5607 }
5608
5609 fn exacthessian_second_directional_derivative_from_evaluation(
5630 &self,
5631 eval: &SpatialAdaptiveExactEvaluation,
5632 direction_u: &Array1<f64>,
5633 direction_v: &Array1<f64>,
5634 ) -> Result<Option<Array2<f64>>, String> {
5635 let p = self.design.ncols();
5636 if eval.obs.neghessian_eta_derivative.iter().any(|&w| w != 0.0) {
5638 return Ok(None);
5639 }
5640 let mut total = Array2::<f64>::zeros((p, p));
5641 for (cache_idx, cache) in self.runtime_caches.iter().enumerate() {
5642 let params = self.adaptive_params.get(cache_idx).ok_or_else(|| {
5643 format!(
5644 "missing adaptive parameter block for cache {}",
5645 cache.termname
5646 )
5647 })?;
5648 let state = eval
5649 .adaptive_states
5650 .get(cache_idx)
5651 .ok_or_else(|| format!("missing adaptive state for cache {}", cache.termname))?;
5652 let u_local = direction_u.slice(s![cache.coeff_global_range.clone()]);
5653 let v_local = direction_v.slice(s![cache.coeff_global_range.clone()]);
5654
5655 let q0_u = cache.d0.dot(&u_local);
5657 let q0_v = cache.d0.dot(&v_local);
5658 let h0 = scalar_operatorhessian(
5659 &cache.d0,
5660 &state.magnitude.second_directionalhessian_diag(&q0_u, &q0_v),
5661 )
5662 .mapv(|x| params.lambda[0] * x);
5663
5664 let a1 = collocationgradient_blocks(&cache.d1.dot(&u_local), cache.dimension)
5666 .map_err(|e| e.to_string())?;
5667 let b1 = collocationgradient_blocks(&cache.d1.dot(&v_local), cache.dimension)
5668 .map_err(|e| e.to_string())?;
5669 let hg = grouped_operatorhessian(
5670 &cache.d1,
5671 cache.dimension,
5672 &state.gradient.second_directionalhessian_blocks(&a1, &b1),
5673 )
5674 .map_err(|e| e.to_string())?
5675 .mapv(|x| params.lambda[1] * x);
5676
5677 let a2 = collocationhessian_blocks(&cache.d2.dot(&u_local), cache.dimension)
5679 .map_err(|e| e.to_string())?;
5680 let b2 = collocationhessian_blocks(&cache.d2.dot(&v_local), cache.dimension)
5681 .map_err(|e| e.to_string())?;
5682 let hc = grouped_operatorhessian(
5683 &cache.d2,
5684 cache.dimension * cache.dimension,
5685 &state.curvature.second_directionalhessian_blocks(&a2, &b2),
5686 )
5687 .map_err(|e| e.to_string())?
5688 .mapv(|x| params.lambda[2] * x);
5689
5690 let range = cache.coeff_global_range.clone();
5691 let mut local = total.slice_mut(s![range.clone(), range]);
5692 local += &h0;
5693 local += &hg;
5694 local += &hc;
5695 }
5696 Ok(Some(total))
5697 }
5698}
5699
5700impl CustomFamily for SpatialAdaptiveExactFamily {
5701 fn joint_jeffreys_term_required(&self) -> bool {
5705 true
5706 }
5707
5708 fn joint_jeffreys_information_with_specs(
5745 &self,
5746 block_states: &[ParameterBlockState],
5747 specs: &[ParameterBlockSpec],
5748 ) -> Result<Option<Array2<f64>>, String> {
5749 let spec = expect_single_blockspec(specs, "spatial adaptive exact family")?;
5750 let beta = &expect_single_block_state(block_states, "spatial adaptive exact family")?.beta;
5751 if spec.design.ncols() != beta.len() {
5752 return Err(SmoothError::dimension_mismatch(format!(
5753 "spatial adaptive Jeffreys information: spec design has {} columns, beta has {}",
5754 spec.design.ncols(),
5755 beta.len()
5756 ))
5757 .into());
5758 }
5759 let eval = self.exact_evaluation(beta)?;
5760 Ok(Some(xt_diag_x_dense(
5761 self.design.view(),
5762 eval.obs.neghessian_eta.view(),
5763 )?))
5764 }
5765
5766 fn joint_jeffreys_information_directional_derivative_with_specs(
5767 &self,
5768 block_states: &[ParameterBlockState],
5769 specs: &[ParameterBlockSpec],
5770 d_beta_flat: &Array1<f64>,
5771 ) -> Result<Option<Array2<f64>>, String> {
5772 let spec = expect_single_blockspec(specs, "spatial adaptive exact family")?;
5778 let beta = &expect_single_block_state(block_states, "spatial adaptive exact family")?.beta;
5779 if spec.design.ncols() != d_beta_flat.len() {
5780 return Err(SmoothError::dimension_mismatch(format!(
5781 "spatial adaptive Jeffreys directional derivative: spec design has {} columns, direction has {}",
5782 spec.design.ncols(),
5783 d_beta_flat.len()
5784 ))
5785 .into());
5786 }
5787 let eval = self.exact_evaluation(beta)?;
5788 let d_eta = gam_linalg::faer_ndarray::fast_av(self.design.as_ref(), d_beta_flat);
5789 Ok(Some(xt_diag_x_dense(
5790 self.design.view(),
5791 (&eval.obs.neghessian_eta_derivative * &d_eta).view(),
5792 )?))
5793 }
5794
5795 fn joint_jeffreys_information_second_directional_derivative_with_specs(
5796 &self,
5797 block_states: &[ParameterBlockState],
5798 specs: &[ParameterBlockSpec],
5799 d_beta_u_flat: &Array1<f64>,
5800 d_betav_flat: &Array1<f64>,
5801 ) -> Result<Option<Array2<f64>>, String> {
5802 let spec = expect_single_blockspec(specs, "spatial adaptive exact family")?;
5809 let beta = &expect_single_block_state(block_states, "spatial adaptive exact family")?.beta;
5810 if spec.design.ncols() != beta.len()
5811 || d_beta_u_flat.len() != beta.len()
5812 || d_betav_flat.len() != beta.len()
5813 {
5814 return Err(SmoothError::dimension_mismatch(format!(
5815 "spatial adaptive Jeffreys second-direction length mismatch: spec cols={}, dirs=({}, {}), expected {}",
5816 spec.design.ncols(),
5817 d_beta_u_flat.len(),
5818 d_betav_flat.len(),
5819 beta.len()
5820 ))
5821 .into());
5822 }
5823 let eval = self.exact_evaluation(beta)?;
5824 if eval.obs.neghessian_eta_derivative.iter().any(|&w| w != 0.0) {
5825 return Ok(None);
5826 }
5827 Ok(Some(Array2::<f64>::zeros((beta.len(), beta.len()))))
5828 }
5829
5830 fn joint_jeffreys_information_matches_observed_hessian(&self) -> bool {
5831 false
5836 }
5837
5838 fn joint_jeffreys_information_depends_on_psi(&self) -> bool {
5839 false
5848 }
5849
5850 fn evaluate(&self, block_states: &[ParameterBlockState]) -> Result<FamilyEvaluation, String> {
5851 let beta = &expect_single_block_state(block_states, "spatial adaptive exact family")?.beta;
5852 let eval = self.exact_evaluation(beta)?;
5853 let mut gradient = fast_atv(&self.design, &eval.obs.score);
5854 gradient -= &eval.total_penaltygradient();
5855 let mut hessian = xt_diag_x_dense(self.design.view(), eval.obs.neghessian_eta.view())?;
5856 hessian += &eval.total_penaltyhessian();
5857 Ok(FamilyEvaluation {
5858 log_likelihood: eval.obs.log_likelihood - eval.total_penalty_value(),
5859 blockworking_sets: vec![BlockWorkingSet::ExactNewton {
5860 gradient,
5861 hessian: SymmetricMatrix::Dense(hessian),
5862 }],
5863 })
5864 }
5865
5866 fn log_likelihood_only(&self, block_states: &[ParameterBlockState]) -> Result<f64, String> {
5867 let state = expect_single_block_state(block_states, "spatial adaptive exact family")?;
5868 let beta = &state.beta;
5869 let obs = evaluate_standard_familyobservations(
5870 self.family.clone(),
5871 self.latent_cloglog_state.as_ref(),
5872 self.mixture_link_state.as_ref(),
5873 self.sas_link_state.as_ref(),
5874 &self.y,
5875 &self.weights,
5876 &state.eta,
5877 )
5878 .map_err(|e| e.to_string())?;
5879 let adaptive_penalty = self.adaptive_penalty_value_only(beta)?;
5880 let (fixed_quadratic, _) = self.fixed_quadratic_terms(beta)?;
5881 Ok(obs.log_likelihood - adaptive_penalty - fixed_quadratic)
5882 }
5883
5884 fn exact_newton_outerobjective(&self) -> ExactNewtonOuterObjective {
5885 ExactNewtonOuterObjective::StrictPseudoLaplace
5886 }
5887
5888 fn exact_newton_joint_hessian(
5889 &self,
5890 block_states: &[ParameterBlockState],
5891 ) -> Result<Option<Array2<f64>>, String> {
5892 let beta = &expect_single_block_state(block_states, "spatial adaptive exact family")?.beta;
5893 let eval = self.exact_evaluation(beta)?;
5894 Ok(Some(eval.totalobjectivehessian(&self.design)?))
5895 }
5896
5897 fn exact_newton_hessian_directional_derivative(
5898 &self,
5899 block_states: &[ParameterBlockState],
5900 block_idx: usize,
5901 d_beta: &Array1<f64>,
5902 ) -> Result<Option<Array2<f64>>, String> {
5903 expect_block_idx_zero(block_idx, "spatial adaptive exact family", "")?;
5904 self.exact_newton_joint_hessian_directional_derivative(block_states, d_beta)
5905 }
5906
5907 fn exact_newton_joint_hessian_directional_derivative(
5908 &self,
5909 block_states: &[ParameterBlockState],
5910 d_beta_flat: &Array1<f64>,
5911 ) -> Result<Option<Array2<f64>>, String> {
5912 let beta = &expect_single_block_state(block_states, "spatial adaptive exact family")?.beta;
5913 if d_beta_flat.len() != beta.len() {
5914 return Err(SmoothError::dimension_mismatch(format!(
5915 "spatial adaptive exact family direction length mismatch: got {}, expected {}",
5916 d_beta_flat.len(),
5917 beta.len()
5918 ))
5919 .into());
5920 }
5921 let eval = self.exact_evaluation(beta)?;
5922 Ok(Some(
5923 self.exacthessian_directional_derivative_from_evaluation(beta, &eval, d_beta_flat)?,
5924 ))
5925 }
5926
5927 fn exact_newton_joint_hessiansecond_directional_derivative(
5928 &self,
5929 block_states: &[ParameterBlockState],
5930 d_beta_u_flat: &Array1<f64>,
5931 d_betav_flat: &Array1<f64>,
5932 ) -> Result<Option<Array2<f64>>, String> {
5933 let beta = &expect_single_block_state(block_states, "spatial adaptive exact family")?.beta;
5934 if d_beta_u_flat.len() != beta.len() || d_betav_flat.len() != beta.len() {
5935 return Err(SmoothError::dimension_mismatch(format!(
5936 "spatial adaptive exact family second-direction length mismatch: got ({}, {}), expected {}",
5937 d_beta_u_flat.len(),
5938 d_betav_flat.len(),
5939 beta.len()
5940 ))
5941 .into());
5942 }
5943 let eval = self.exact_evaluation(beta)?;
5944 self.exacthessian_second_directional_derivative_from_evaluation(
5945 &eval,
5946 d_beta_u_flat,
5947 d_betav_flat,
5948 )
5949 }
5950
5951 fn block_linear_constraints(
5952 &self,
5953 block_states: &[ParameterBlockState],
5954 block_idx: usize,
5955 block_spec: &ParameterBlockSpec,
5956 ) -> Result<Option<ConstraintSet>, String> {
5957 assert!(!block_states.is_empty(), "block_states must be non-empty");
5958 assert!(
5959 !block_spec.name.is_empty(),
5960 "block spec name must be non-empty",
5961 );
5962 expect_block_idx_zero(block_idx, "spatial adaptive exact family", "")?;
5963 Ok(self.linear_constraints.clone().map(ConstraintSet::Dense))
5964 }
5965
5966 fn exact_newton_joint_psi_terms(
5967 &self,
5968 block_states: &[ParameterBlockState],
5969 specs: &[ParameterBlockSpec],
5970 hyper_layout: &crate::custom_family::CustomFamilyHyperLayout,
5971 psi_index: usize,
5972 ) -> Result<Option<ExactNewtonJointPsiTerms>, String> {
5973 if hyper_layout.family_axis_count() != 0 {
5974 return Err(
5975 "spatial adaptive exact family does not declare family-owned hyper axes"
5976 .to_string(),
5977 );
5978 }
5979 let derivative_blocks = hyper_layout.design_derivative_blocks();
5980 if block_states.len() != 1 || specs.len() != 1 || derivative_blocks.len() != 1 {
5981 return Err(SmoothError::dimension_mismatch(format!(
5982 "spatial adaptive exact family expects one block/state/spec/psi payload, got states={} specs={} deriv_blocks={}",
5983 block_states.len(),
5984 specs.len(),
5985 derivative_blocks.len()
5986 ))
5987 .into());
5988 }
5989 derivative_blocks[0]
5990 .get(psi_index)
5991 .ok_or_else(|| format!("adaptive psi index {} out of bounds", psi_index))?;
5992 let hyper = self
5993 .hyperspecs
5994 .get(psi_index)
5995 .ok_or_else(|| format!("adaptive psi index {} out of bounds", psi_index))?;
5996 let beta = &block_states[0].beta;
5997 let eval = self.exact_evaluation(beta)?;
5998 let (direct, beta_mixed, betahessian_explicit) =
5999 self.adaptive_hyper_parts(&eval, *hyper)?;
6000
6001 Ok(Some(ExactNewtonJointPsiTerms {
6022 objective_psi: direct,
6023 score_psi: beta_mixed,
6024 hessian_psi: betahessian_explicit,
6025 hessian_psi_operator: None,
6026 }))
6027 }
6028
6029 fn exact_newton_joint_psisecond_order_terms(
6030 &self,
6031 block_states: &[ParameterBlockState],
6032 specs: &[ParameterBlockSpec],
6033 hyper_layout: &crate::custom_family::CustomFamilyHyperLayout,
6034 psi_i: usize,
6035 psi_j: usize,
6036 ) -> Result<Option<gam_problem::ExactNewtonJointPsiSecondOrderTerms>, String> {
6037 if hyper_layout.family_axis_count() != 0 {
6038 return Err(
6039 "spatial adaptive exact family does not declare family-owned hyper axes"
6040 .to_string(),
6041 );
6042 }
6043 let derivative_blocks = hyper_layout.design_derivative_blocks();
6044 if block_states.len() != 1 || specs.len() != 1 || derivative_blocks.len() != 1 {
6045 return Err(SmoothError::dimension_mismatch(format!(
6046 "spatial adaptive exact family expects one block/state/spec/psi payload, got states={} specs={} deriv_blocks={}",
6047 block_states.len(),
6048 specs.len(),
6049 derivative_blocks.len()
6050 ))
6051 .into());
6052 }
6053 derivative_blocks[0]
6054 .get(psi_i)
6055 .ok_or_else(|| format!("adaptive psi index {} out of bounds", psi_i))?;
6056 derivative_blocks[0]
6057 .get(psi_j)
6058 .ok_or_else(|| format!("adaptive psi index {} out of bounds", psi_j))?;
6059 let hyper_i = self
6060 .hyperspecs
6061 .get(psi_i)
6062 .ok_or_else(|| format!("adaptive psi index {} out of bounds", psi_i))?;
6063 let hyper_j = self
6064 .hyperspecs
6065 .get(psi_j)
6066 .ok_or_else(|| format!("adaptive psi index {} out of bounds", psi_j))?;
6067 let beta = &block_states[0].beta;
6068 let eval = self.exact_evaluation(beta)?;
6069 let (objective_psi_psi, score_psi_psi, hessian_psi_psi) =
6070 self.adaptive_explicit_second_order_parts(&eval, *hyper_i, *hyper_j)?;
6071
6072 Ok(Some(gam_problem::ExactNewtonJointPsiSecondOrderTerms {
6073 objective_psi_psi,
6074 score_psi_psi,
6075 hessian_psi_psi,
6076 hessian_psi_psi_operator: None,
6077 }))
6078 }
6079
6080 fn exact_newton_joint_psihessian_directional_derivative(
6081 &self,
6082 block_states: &[ParameterBlockState],
6083 specs: &[ParameterBlockSpec],
6084 hyper_layout: &crate::custom_family::CustomFamilyHyperLayout,
6085 psi_index: usize,
6086 direction: &Array1<f64>,
6087 ) -> Result<Option<Array2<f64>>, String> {
6088 if hyper_layout.family_axis_count() != 0 {
6089 return Err(
6090 "spatial adaptive exact family does not declare family-owned hyper axes"
6091 .to_string(),
6092 );
6093 }
6094 let derivative_blocks = hyper_layout.design_derivative_blocks();
6095 if block_states.len() != 1 || specs.len() != 1 || derivative_blocks.len() != 1 {
6096 return Err(SmoothError::dimension_mismatch(format!(
6097 "spatial adaptive exact family expects one block/state/spec/psi payload, got states={} specs={} deriv_blocks={}",
6098 block_states.len(),
6099 specs.len(),
6100 derivative_blocks.len()
6101 ))
6102 .into());
6103 }
6104 let beta = &block_states[0].beta;
6105 if direction.len() != beta.len() {
6106 return Err(SmoothError::dimension_mismatch(format!(
6107 "spatial adaptive exact family direction length mismatch: got {}, expected {}",
6108 direction.len(),
6109 beta.len()
6110 ))
6111 .into());
6112 }
6113 derivative_blocks[0]
6114 .get(psi_index)
6115 .ok_or_else(|| format!("adaptive psi index {} out of bounds", psi_index))?;
6116 let hyper = self
6117 .hyperspecs
6118 .get(psi_index)
6119 .ok_or_else(|| format!("adaptive psi index {} out of bounds", psi_index))?;
6120 let eval = self.exact_evaluation(beta)?;
6121 let drift = match hyper.kind {
6122 SpatialAdaptiveHyperKind::LogLambdaMagnitude
6123 | SpatialAdaptiveHyperKind::LogLambdaGradient
6124 | SpatialAdaptiveHyperKind::LogLambdaCurvature => self.adaptive_block_drift_eval(
6125 &eval,
6126 hyper.cache_index,
6127 AdaptiveComponent::from_index(hyper.kind.component_index())?,
6128 HyperDriftKind::Rho,
6129 direction,
6130 )?,
6131 SpatialAdaptiveHyperKind::LogEpsilonMagnitude
6132 | SpatialAdaptiveHyperKind::LogEpsilonGradient
6133 | SpatialAdaptiveHyperKind::LogEpsilonCurvature => self
6134 .adaptive_shared_log_epsilon_drift(
6135 &eval,
6136 hyper.kind.component_index(),
6137 direction,
6138 )?,
6139 };
6140 Ok(Some(drift))
6141 }
6142}
6143
6144fn expect_single_block_state<'a>(
6145 block_states: &'a [ParameterBlockState],
6146 family_name: &str,
6147) -> Result<&'a ParameterBlockState, String> {
6148 crate::block_layout::block_count::validate_block_count::<SmoothError>(
6149 family_name,
6150 1,
6151 block_states.len(),
6152 )?;
6153 Ok(&block_states[0])
6154}
6155
6156fn expect_single_blockspec<'a>(
6157 specs: &'a [ParameterBlockSpec],
6158 family_name: &str,
6159) -> Result<&'a ParameterBlockSpec, String> {
6160 crate::block_layout::block_count::validate_block_count::<SmoothError>(
6161 family_name,
6162 1,
6163 specs.len(),
6164 )?;
6165 Ok(&specs[0])
6166}
6167
6168fn expect_block_idx_zero(block_idx: usize, family_name: &str, context: &str) -> Result<(), String> {
6169 if block_idx != 0 {
6170 return Err(SmoothError::invalid_index(format!(
6171 "{family_name} expects block_idx 0{context}, got {block_idx}"
6172 ))
6173 .into());
6174 }
6175 Ok::<(), _>(())
6176}
6177
6178impl BoundedLinearFamily {
6179 fn bounded_term_derivative_data(
6180 &self,
6181 latent_beta: &Array1<f64>,
6182 ) -> Result<
6183 (
6184 Array1<f64>,
6185 Array1<f64>,
6186 Array1<f64>,
6187 Array1<f64>,
6188 Array1<f64>,
6189 ),
6190 String,
6191 > {
6192 let p = latent_beta.len();
6193 if p != self.design.ncols() || latent_beta.iter().any(|value| !value.is_finite()) {
6194 return Err(format!(
6195 "bounded coefficient geometry requires {} finite latent coefficients, got {}",
6196 self.design.ncols(),
6197 p
6198 ));
6199 }
6200 let mut beta_user = latent_beta.clone();
6201 let mut jac_diag = Array1::<f64>::ones(p);
6202 let mut second_diag = Array1::<f64>::zeros(p);
6203 let mut third_diag = Array1::<f64>::zeros(p);
6204 let mut priorthird = Array1::<f64>::zeros(p);
6205 for term in &self.bounded_terms {
6206 let width = term.max - term.min;
6207 if term.col_idx >= p
6208 || !term.min.is_finite()
6209 || !term.max.is_finite()
6210 || !(width.is_finite() && width > 0.0)
6211 {
6212 return Err(format!(
6213 "bounded coefficient geometry has invalid column/bounds: col={}, p={p}, bounds=({}, {})",
6214 term.col_idx, term.min, term.max
6215 ));
6216 }
6217 let (beta, _, db_dtheta, d2b_dtheta2, d3b_dtheta3) =
6218 bounded_latent_derivatives(latent_beta[term.col_idx], term.min, term.max);
6219 if [beta, db_dtheta, d2b_dtheta2, d3b_dtheta3]
6220 .iter()
6221 .any(|value| !value.is_finite())
6222 {
6223 return Err(format!(
6224 "bounded coefficient transform is not representable at column {} and theta={}",
6225 term.col_idx, latent_beta[term.col_idx]
6226 ));
6227 }
6228 beta_user[term.col_idx] = beta;
6229 jac_diag[term.col_idx] = db_dtheta;
6230 second_diag[term.col_idx] = d2b_dtheta2;
6231 third_diag[term.col_idx] = d3b_dtheta3;
6232 let (_, _, _, prior_neghess_derivative) =
6233 bounded_prior_terms(latent_beta[term.col_idx], &term.prior)?;
6234 priorthird[term.col_idx] = prior_neghess_derivative;
6235 }
6236 Ok((beta_user, jac_diag, second_diag, third_diag, priorthird))
6237 }
6238
6239 fn user_beta_and_jacobian(
6240 &self,
6241 latent_beta: &Array1<f64>,
6242 ) -> Result<(Array1<f64>, Array1<f64>), String> {
6243 let (beta_user, jac_diag, _, _, _) = self.bounded_term_derivative_data(latent_beta)?;
6244 Ok((beta_user, jac_diag))
6245 }
6246
6247 fn nonlinear_offset_from_latent(
6248 &self,
6249 latent_beta: &Array1<f64>,
6250 ) -> Result<Array1<f64>, String> {
6251 self.bounded_term_derivative_data(latent_beta)?;
6252 let mut offset = self.offset.clone();
6253 for term in &self.bounded_terms {
6254 let (beta, _, _) =
6255 bounded_latent_to_user(latent_beta[term.col_idx], term.min, term.max);
6256 offset.scaled_add(beta, &self.design.column(term.col_idx));
6257 }
6258 if offset.iter().any(|value| !value.is_finite()) {
6259 return Err("bounded nonlinear offset is not representable".to_string());
6260 }
6261 Ok(offset)
6262 }
6263
6264 fn effective_design_for_latent(&self, jac_diag: &Array1<f64>) -> Array2<f64> {
6265 let mut x_eff = self.design.clone();
6266 for term in &self.bounded_terms {
6267 x_eff
6268 .column_mut(term.col_idx)
6269 .mapv_inplace(|v| v * jac_diag[term.col_idx]);
6270 }
6271 x_eff
6272 }
6273
6274 fn exacthessian_andgradient(
6275 &self,
6276 latent_beta: &Array1<f64>,
6277 ) -> Result<
6278 (
6279 StandardFamilyObservationState,
6280 Array2<f64>,
6281 Array1<f64>,
6282 f64,
6283 Array1<f64>,
6284 Array1<f64>,
6285 Array1<f64>,
6286 ),
6287 String,
6288 > {
6289 let (_, jac_diag, second_diag, third_diag, priorthird) =
6290 self.bounded_term_derivative_data(latent_beta)?;
6291 let x_eff = self.effective_design_for_latent(&jac_diag);
6292 let eta =
6293 self.designzeroed.dot(latent_beta) + self.nonlinear_offset_from_latent(latent_beta)?;
6294 let obs = evaluate_resolved_standard_family_observations(
6295 &self.likelihood,
6296 self.latent_cloglog_state.as_ref(),
6297 self.mixture_link_state.as_ref(),
6298 self.sas_link_state.as_ref(),
6299 &self.y,
6300 &self.weights,
6301 &eta,
6302 )
6303 .map_err(|e| e.to_string())?;
6304
6305 let mut priorgrad = Array1::<f64>::zeros(latent_beta.len());
6306 let mut prior_neghess = Array2::<f64>::zeros((latent_beta.len(), latent_beta.len()));
6307 let mut prior_loglik = 0.0;
6308 for term in &self.bounded_terms {
6309 let (logp, grad, neghess, _) =
6310 bounded_prior_terms(latent_beta[term.col_idx], &term.prior)?;
6311 prior_loglik += logp;
6312 priorgrad[term.col_idx] += grad;
6313 prior_neghess[[term.col_idx, term.col_idx]] += neghess;
6314 }
6315
6316 let mut hessian = xt_diag_x_dense(x_eff.view(), obs.neghessian_eta.view())?;
6317 let mut gradient = fast_atv(&x_eff, &obs.score);
6318 for term in &self.bounded_terms {
6319 let score_beta = self.design.column(term.col_idx).dot(&obs.score);
6320 hessian[[term.col_idx, term.col_idx]] -= score_beta * second_diag[term.col_idx];
6321 }
6322 hessian += &prior_neghess;
6323 gradient += &priorgrad;
6324
6325 Ok((
6326 obs,
6327 hessian,
6328 gradient,
6329 prior_loglik,
6330 second_diag,
6331 third_diag,
6332 priorthird,
6333 ))
6334 }
6335
6336 fn evaluation_from_latent(
6337 &self,
6338 latent_beta: &Array1<f64>,
6339 ) -> Result<
6340 (
6341 StandardFamilyObservationState,
6342 Array2<f64>,
6343 Array1<f64>,
6344 f64,
6345 ),
6346 String,
6347 > {
6348 let (obs, hessian, gradient, prior_loglik, _, _, _) =
6349 self.exacthessian_andgradient(latent_beta)?;
6350 Ok((obs, hessian, gradient, prior_loglik))
6351 }
6352}
6353
6354impl CustomFamily for BoundedLinearFamily {
6355 fn joint_jeffreys_term_required(&self) -> bool {
6359 true
6360 }
6361
6362 fn evaluate(&self, block_states: &[ParameterBlockState]) -> Result<FamilyEvaluation, String> {
6363 let latent_beta = &expect_single_block_state(block_states, "bounded linear family")?.beta;
6364 let (obs, hessian, gradient, prior_loglik) = self.evaluation_from_latent(latent_beta)?;
6365 Ok(FamilyEvaluation {
6366 log_likelihood: obs.log_likelihood + prior_loglik,
6367 blockworking_sets: vec![BlockWorkingSet::ExactNewton {
6368 gradient,
6369 hessian: SymmetricMatrix::Dense(hessian),
6370 }],
6371 })
6372 }
6373
6374 fn exact_newton_joint_hessian(
6375 &self,
6376 block_states: &[ParameterBlockState],
6377 ) -> Result<Option<Array2<f64>>, String> {
6378 let latent_beta = &expect_single_block_state(block_states, "bounded linear family")?.beta;
6379 let (_, hessian, _, _) = self.evaluation_from_latent(latent_beta)?;
6380 Ok(Some(hessian))
6381 }
6382
6383 fn exact_newton_hessian_directional_derivative(
6384 &self,
6385 block_states: &[ParameterBlockState],
6386 block_idx: usize,
6387 d_beta: &Array1<f64>,
6388 ) -> Result<Option<Array2<f64>>, String> {
6389 expect_block_idx_zero(block_idx, "bounded linear family", "")?;
6390 self.exact_newton_joint_hessian_directional_derivative(block_states, d_beta)
6391 }
6392
6393 fn exact_newton_joint_hessian_directional_derivative(
6394 &self,
6395 block_states: &[ParameterBlockState],
6396 d_beta_flat: &Array1<f64>,
6397 ) -> Result<Option<Array2<f64>>, String> {
6398 let latent_beta = &expect_single_block_state(block_states, "bounded linear family")?.beta;
6399 if d_beta_flat.len() != latent_beta.len() {
6400 return Err(SmoothError::dimension_mismatch(format!(
6401 "bounded linear family directional derivative length mismatch: got {}, expected {}",
6402 d_beta_flat.len(),
6403 latent_beta.len()
6404 ))
6405 .into());
6406 }
6407
6408 let (obs, _, _, _, second_diag, third_diag, priorthird) =
6409 self.exacthessian_andgradient(latent_beta)?;
6410
6411 let (_, jac_diag, _, _, _) = self.bounded_term_derivative_data(latent_beta)?;
6412 let x_eff = self.effective_design_for_latent(&jac_diag);
6413 let deta = x_eff.dot(d_beta_flat);
6414 let d_neghess_eta = &obs.neghessian_eta_derivative * &deta;
6415
6416 let mut dx_eff = Array2::<f64>::zeros(x_eff.raw_dim());
6417 for term in &self.bounded_terms {
6418 let scale = second_diag[term.col_idx] * d_beta_flat[term.col_idx];
6419 if scale != 0.0 {
6420 let mut col = dx_eff.column_mut(term.col_idx);
6421 col.assign(&self.design.column(term.col_idx));
6422 col.mapv_inplace(|v| v * scale);
6423 }
6424 }
6425
6426 let mut dhessian = xt_diag_x_dense(x_eff.view(), d_neghess_eta.view())?;
6427 let mut wxdx = Array2::<f64>::zeros((x_eff.ncols(), x_eff.ncols()));
6428 for i in 0..x_eff.nrows() {
6429 let wi = obs.neghessian_eta[i];
6430 if wi == 0.0 {
6431 continue;
6432 }
6433 for a in 0..x_eff.ncols() {
6434 let xa = x_eff[[i, a]];
6435 for b in 0..x_eff.ncols() {
6436 wxdx[[a, b]] += wi * (dx_eff[[i, a]] * x_eff[[i, b]] + xa * dx_eff[[i, b]]);
6437 }
6438 }
6439 }
6440 dhessian += &wxdx;
6441
6442 let d_score = -&obs.neghessian_eta * &deta;
6443 for term in &self.bounded_terms {
6444 let score_beta = self.design.column(term.col_idx).dot(&obs.score);
6445 let d_score_beta = self.design.column(term.col_idx).dot(&d_score);
6446 dhessian[[term.col_idx, term.col_idx]] -= d_score_beta * second_diag[term.col_idx]
6447 + score_beta * third_diag[term.col_idx] * d_beta_flat[term.col_idx];
6448 dhessian[[term.col_idx, term.col_idx]] +=
6449 priorthird[term.col_idx] * d_beta_flat[term.col_idx];
6450 }
6451
6452 Ok(Some(dhessian))
6453 }
6454
6455 fn block_geometry(
6456 &self,
6457 block_states: &[ParameterBlockState],
6458 spec: &ParameterBlockSpec,
6459 ) -> Result<(DesignMatrix, Array1<f64>), String> {
6460 if block_states.is_empty() {
6461 return Ok((
6462 DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(
6463 self.designzeroed.clone(),
6464 )),
6465 self.offset.clone(),
6466 ));
6467 }
6468 let offset = self.nonlinear_offset_from_latent(
6469 &expect_single_block_state(block_states, "bounded linear family")?.beta,
6470 )?;
6471 let x = if spec.design.ncols() == self.designzeroed.ncols() {
6472 self.designzeroed.clone()
6473 } else {
6474 return Err(SmoothError::dimension_mismatch(
6475 "bounded linear family design column mismatch",
6476 )
6477 .into());
6478 };
6479 Ok((
6480 DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(x)),
6481 offset,
6482 ))
6483 }
6484
6485 fn block_geometry_is_dynamic(&self) -> bool {
6486 true
6487 }
6488
6489 fn block_geometry_directional_derivative(
6490 &self,
6491 block_states: &[ParameterBlockState],
6492 block_idx: usize,
6493 spec: &ParameterBlockSpec,
6494 d_beta: &Array1<f64>,
6495 ) -> Result<Option<BlockGeometryDirectionalDerivative>, String> {
6496 expect_block_idx_zero(
6497 block_idx,
6498 "bounded linear family",
6499 " for geometry derivative",
6500 )?;
6501 expect_single_block_state(block_states, "bounded linear family")?;
6502 if d_beta.len() != spec.design.ncols() {
6503 return Err(SmoothError::dimension_mismatch(format!(
6504 "bounded linear family geometry derivative direction mismatch: got {}, expected {}",
6505 d_beta.len(),
6506 spec.design.ncols()
6507 ))
6508 .into());
6509 }
6510 let (_, jac_diag, _, _, _) = self.bounded_term_derivative_data(&block_states[0].beta)?;
6511 let mut d_offset = Array1::<f64>::zeros(self.offset.len());
6512 let has_drift = self
6513 .bounded_terms
6514 .iter()
6515 .any(|term| jac_diag[term.col_idx] != 0.0 && d_beta[term.col_idx] != 0.0);
6516 if !has_drift {
6517 return Ok(Some(BlockGeometryDirectionalDerivative {
6518 d_design: None,
6519 d_offset,
6520 }));
6521 }
6522 for term in &self.bounded_terms {
6523 let col = term.col_idx;
6524 let drift = jac_diag[col] * d_beta[col];
6525 if drift != 0.0 {
6526 d_offset.scaled_add(drift, &self.design.column(col));
6527 }
6528 }
6529 Ok(Some(BlockGeometryDirectionalDerivative {
6530 d_design: None,
6531 d_offset,
6532 }))
6533 }
6534}
6535
6536#[inline]
6537fn dense_diag_gram_chunkrows(p: usize) -> usize {
6538 const MIN_ROWS: usize = 512;
6539 const MAX_ROWS: usize = 2048;
6540 const TARGET_BYTES: usize = 2 * 1024 * 1024;
6541 let bytes_per_row = p.max(1) * std::mem::size_of::<f64>();
6542 (TARGET_BYTES / bytes_per_row).clamp(MIN_ROWS, MAX_ROWS)
6543}
6544
6545fn xt_diag_x_dense(x: ArrayView2<'_, f64>, w: ArrayView1<'_, f64>) -> Result<Array2<f64>, String> {
6546 if x.nrows() != w.len() {
6547 return Err(SmoothError::dimension_mismatch("xt_diag_x_dense row mismatch").into());
6548 }
6549 let (n, p) = x.dim();
6550 if n == 0 || p == 0 {
6551 return Ok(Array2::<f64>::zeros((p, p)));
6552 }
6553
6554 const STREAMING_BYTES_THRESHOLD: usize = 8 * 1024 * 1024;
6555 let dense_work_bytes = n
6556 .checked_mul(p)
6557 .and_then(|cells| cells.checked_mul(std::mem::size_of::<f64>()))
6558 .unwrap_or(usize::MAX);
6559 if dense_work_bytes <= STREAMING_BYTES_THRESHOLD {
6560 let mut weighted = x.to_owned();
6561 ndarray::Zip::from(weighted.rows_mut())
6562 .and(w)
6563 .par_for_each(|mut row, wi| row *= *wi);
6564 return Ok(fast_atb(&x, &weighted));
6565 }
6566
6567 let chunkrows = dense_diag_gram_chunkrows(p).min(n);
6568 let mut weighted_chunk = Array2::<f64>::zeros((chunkrows, p));
6569 let mut out = Array2::<f64>::zeros((p, p));
6570 for row_start in (0..n).step_by(chunkrows) {
6571 let rows = (n - row_start).min(chunkrows);
6572 let x_chunk = x.slice(s![row_start..row_start + rows, ..]);
6573 {
6574 let mut chunk = weighted_chunk.slice_mut(s![0..rows, ..]);
6575 for local_row in 0..rows {
6576 let scale = w[row_start + local_row];
6577 if scale == 0.0 {
6578 chunk.row_mut(local_row).fill(0.0);
6579 continue;
6580 }
6581 for col in 0..p {
6582 chunk[[local_row, col]] = x_chunk[[local_row, col]] * scale;
6583 }
6584 }
6585 }
6586 out += &fast_atb(&x_chunk, &weighted_chunk.slice(s![0..rows, ..]));
6587 }
6588 Ok(out)
6589}
6590
6591fn trace_of_dense_product(a: &Array2<f64>, b: &Array2<f64>) -> Result<f64, String> {
6592 if a.nrows() != a.ncols() || b.nrows() != b.ncols() || a.nrows() != b.nrows() {
6593 return Err(
6594 SmoothError::dimension_mismatch("trace_of_dense_product dimension mismatch").into(),
6595 );
6596 }
6597 if a.iter().chain(b.iter()).any(|value| !value.is_finite()) {
6598 return Err("trace_of_dense_product requires finite matrices".to_string());
6599 }
6600 let mut trace = gam_linalg::utils::KahanSum::default();
6601 for i in 0..a.nrows() {
6602 for j in 0..a.ncols() {
6603 let term = a[[i, j]] * b[[j, i]];
6604 if !term.is_finite() {
6605 return Err(format!(
6606 "trace_of_dense_product term ({i}, {j}) is not representable"
6607 ));
6608 }
6609 trace.add(term);
6610 }
6611 }
6612 let trace = trace.sum();
6613 if !trace.is_finite() {
6614 return Err("trace_of_dense_product sum is not representable".to_string());
6615 }
6616 Ok(trace)
6617}
6618
6619fn certify_bounded_edf_interval(
6620 value: f64,
6621 lower: f64,
6622 upper: f64,
6623 dimension: usize,
6624 label: &str,
6625) -> Result<f64, EstimationError> {
6626 if !(value.is_finite() && lower.is_finite() && upper.is_finite() && lower <= upper) {
6627 crate::bail_invalid_estim!(
6628 "{label} has invalid EDF interval/value: value={value}, interval=[{lower}, {upper}]"
6629 );
6630 }
6631 let scale = 1.0_f64.max(value.abs()).max(lower.abs()).max(upper.abs());
6632 let allowed = 256.0 * f64::EPSILON * (dimension.max(1) as f64).powi(2) * scale;
6636 if value < lower {
6637 if lower - value <= allowed {
6638 return Ok(lower);
6639 }
6640 } else if value > upper {
6641 if value - upper <= allowed {
6642 return Ok(upper);
6643 }
6644 } else {
6645 return Ok(value);
6646 }
6647 crate::bail_invalid_estim!(
6648 "{label}={value} lies outside [{lower}, {upper}] by more than the dense-trace backward-error allowance {allowed}"
6649 )
6650}
6651
6652fn exact_bounded_edf(
6653 penalties: &[PenaltySpec],
6654 lambdas: &Array1<f64>,
6655 latent_cov: &Array2<f64>,
6656) -> Result<(Vec<f64>, Vec<f64>, f64), EstimationError> {
6657 if penalties.len() != lambdas.len() {
6658 crate::bail_invalid_estim!(
6659 "bounded EDF penalty/lambda mismatch: {} penalties vs {} lambdas",
6660 penalties.len(),
6661 lambdas.len()
6662 );
6663 }
6664 if latent_cov.nrows() != latent_cov.ncols() {
6665 crate::bail_invalid_estim!("bounded EDF covariance must be square");
6666 }
6667
6668 let p = latent_cov.nrows();
6669 let mut s_lambda = Array2::<f64>::zeros((p, p));
6670 let mut edf_by_block = Vec::with_capacity(penalties.len());
6671 let mut penalty_block_trace = Vec::with_capacity(penalties.len());
6673 let mut trace_sum = gam_linalg::utils::KahanSum::default();
6674
6675 for (k, ps) in penalties.iter().enumerate() {
6676 let lambda_k = lambdas[k];
6677 if !(lambda_k.is_finite() && lambda_k >= 0.0) {
6678 crate::bail_invalid_estim!(
6679 "bounded EDF smoothing strength at block {k} must be finite and non-negative, got {lambda_k}"
6680 );
6681 }
6682 match ps {
6683 PenaltySpec::Block {
6684 local, col_range, ..
6685 } => {
6686 s_lambda
6687 .slice_mut(ndarray::s![col_range.clone(), col_range.clone()])
6688 .scaled_add(lambda_k, local);
6689 let penalty_rank =
6691 local
6692 .nrows()
6693 .saturating_sub(estimate_penalty_nullity(local).map_err(|e| {
6694 EstimationError::InvalidInput(format!("bounded EDF rank failed: {e}"))
6695 })?);
6696 let cov_block = latent_cov.slice(ndarray::s![col_range.clone(), col_range.clone()]);
6698 let trace_k = lambda_k
6699 * trace_of_dense_product(&cov_block.to_owned(), local)
6700 .map_err(EstimationError::InvalidInput)?;
6701 trace_sum.add(trace_k);
6702 penalty_block_trace.push(trace_k);
6703 let p_k = penalty_rank as f64;
6704 edf_by_block.push(certify_bounded_edf_interval(
6705 p_k - trace_k,
6706 0.0,
6707 p_k,
6708 p,
6709 &format!("bounded EDF block {k}"),
6710 )?);
6711 }
6712 PenaltySpec::Dense(m) | PenaltySpec::DenseWithMean { matrix: m, .. } => {
6713 s_lambda.scaled_add(lambda_k, m);
6714 let penalty_rank = p.saturating_sub(estimate_penalty_nullity(m).map_err(|e| {
6715 EstimationError::InvalidInput(format!("bounded EDF rank failed: {e}"))
6716 })?);
6717 let trace_k = lambda_k
6718 * trace_of_dense_product(latent_cov, m)
6719 .map_err(EstimationError::InvalidInput)?;
6720 trace_sum.add(trace_k);
6721 penalty_block_trace.push(trace_k);
6722 let p_k = penalty_rank as f64;
6723 edf_by_block.push(certify_bounded_edf_interval(
6724 p_k - trace_k,
6725 0.0,
6726 p_k,
6727 p,
6728 &format!("bounded EDF block {k}"),
6729 )?);
6730 }
6731 }
6732 }
6733
6734 let nullity_total = estimate_penalty_nullity(&s_lambda)
6735 .map_err(|e| EstimationError::InvalidInput(format!("bounded EDF nullity failed: {e}")))?
6736 as f64;
6737 let trace_sum = trace_sum.sum();
6738 let edf_total = certify_bounded_edf_interval(
6739 p as f64 - trace_sum,
6740 nullity_total,
6741 p as f64,
6742 p,
6743 "bounded total EDF",
6744 )?;
6745 Ok((edf_by_block, penalty_block_trace, edf_total))
6746}
6747
6748fn certified_bounded_posterior_covariance(
6754 precision: &Array2<f64>,
6755 label: &'static str,
6756) -> Result<Array2<f64>, EstimationError> {
6757 gam_linalg::utils::certified_spd_inverse(precision, label)
6758 .map(gam_linalg::utils::CertifiedSpdInverse::into_inverse)
6759 .map_err(|error| {
6760 EstimationError::InvalidInput(format!(
6761 "bounded posterior covariance requires an exact SPD precision: {error}"
6762 ))
6763 })
6764}
6765
6766fn transform_bounded_latent_precision_to_user_internal(
6767 latent_precision: &Array2<f64>,
6768 jac_diag: &Array1<f64>,
6769) -> Result<Array2<f64>, EstimationError> {
6770 let p = latent_precision.nrows();
6771 if latent_precision.ncols() != p || jac_diag.len() != p {
6772 crate::bail_invalid_estim!(
6773 "bounded precision transform dimension mismatch: precision is {}x{}, jacobian has {} entries",
6774 latent_precision.nrows(),
6775 latent_precision.ncols(),
6776 jac_diag.len()
6777 );
6778 }
6779 let mut out = latent_precision.clone();
6780 for i in 0..p {
6781 let scale = jac_diag[i];
6782 if !scale.is_finite() || scale <= 0.0 {
6783 crate::bail_invalid_estim!(
6784 "bounded precision transform requires a positive finite coefficient jacobian; column {i} has {scale}"
6785 );
6786 }
6787 if scale != 1.0 {
6788 out.row_mut(i).mapv_inplace(|v| v / scale);
6789 out.column_mut(i).mapv_inplace(|v| v / scale);
6790 }
6791 }
6792 Ok(out)
6793}
6794
6795fn fit_bounded_term_collection_with_design(
6796 y: ArrayView1<'_, f64>,
6797 weights: ArrayView1<'_, f64>,
6798 offset: ArrayView1<'_, f64>,
6799 spec: &TermCollectionSpec,
6800 design: &TermCollectionDesign,
6801 heuristic_lambdas: Option<&[f64]>,
6802 family: LikelihoodSpec,
6803 options: &FitOptions,
6804) -> Result<FittedTermCollection, EstimationError> {
6805 let conditioning_cols: Vec<usize> = spec
6806 .linear_terms
6807 .iter()
6808 .enumerate()
6809 .filter_map(|(j, linear)| {
6810 (!linear.double_penalty).then_some(design.intercept_range.end + j)
6811 })
6812 .collect();
6813 let conditioning = LinearFitConditioning::from_columns(design, &conditioning_cols);
6814 let dense_design = design.design.to_dense_cow();
6815 let fit_design = conditioning.apply_to_design(&dense_design);
6816 let fit_penalties = conditioning
6817 .transform_blockwise_penalties_to_internal(&design.penalties, design.design.ncols());
6818 if design.linear_constraints.is_some() {
6819 crate::bail_invalid_estim!(
6820 "bounded() terms are not yet compatible with explicit linear constraints"
6821 );
6822 }
6823 let mut bounded_terms = Vec::<BoundedLinearTermMeta>::new();
6824 for (j, term) in spec.linear_terms.iter().enumerate() {
6825 if term.double_penalty
6826 && matches!(
6827 term.coefficient_geometry,
6828 LinearCoefficientGeometry::Bounded { .. }
6829 )
6830 {
6831 crate::bail_invalid_estim!(
6832 "bounded linear term '{}' cannot also use double_penalty",
6833 term.name
6834 );
6835 }
6836 if let LinearCoefficientGeometry::Bounded { min, max, prior } =
6837 term.coefficient_geometry.clone()
6838 {
6839 let col_idx = design.intercept_range.end + j;
6840 let (min_internal, max_internal) = conditioning.internal_bounds_for(col_idx, min, max);
6841 bounded_terms.push(BoundedLinearTermMeta {
6842 col_idx,
6843 min: min_internal,
6844 max: max_internal,
6845 prior,
6846 });
6847 }
6848 }
6849 if bounded_terms.is_empty() {
6850 crate::bail_invalid_estim!("internal bounded fit path called with no bounded terms");
6851 }
6852
6853 let mut designzeroed = fit_design.clone();
6854 let mut initial_beta = Array1::<f64>::zeros(fit_design.ncols());
6855 for term in &bounded_terms {
6856 designzeroed.column_mut(term.col_idx).fill(0.0);
6857 initial_beta[term.col_idx] = 0.0;
6858 }
6859
6860 let initial_log_lambdas = heuristic_lambdas
6861 .map(|vals| Array1::from_vec(vals.to_vec()))
6862 .unwrap_or_else(|| Array1::zeros(fit_penalties.len()));
6863 if initial_log_lambdas.len() != fit_penalties.len() {
6864 crate::bail_invalid_estim!(
6865 "heuristic lambda length mismatch for bounded model: got {}, expected {}",
6866 initial_log_lambdas.len(),
6867 fit_penalties.len()
6868 );
6869 }
6870
6871 let glm_likelihood = gam_spec::GlmLikelihoodSpec::canonical(family);
6872 let resolved_likelihood_scale = glm_likelihood
6873 .resolved_scale()
6874 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
6875 let is_beta_logistic = glm_likelihood.spec.is_binomial_beta_logistic();
6876 let family_adapter = BoundedLinearFamily {
6877 likelihood: glm_likelihood.clone(),
6878 latent_cloglog_state: options.latent_cloglog,
6879 mixture_link_state: options
6880 .mixture_link
6881 .clone()
6882 .as_ref()
6883 .map(state_fromspec)
6884 .transpose()
6885 .map_err(EstimationError::InvalidInput)?,
6886 sas_link_state: options
6887 .sas_link
6888 .map(|spec| {
6889 if is_beta_logistic {
6890 state_from_beta_logisticspec(spec)
6891 } else {
6892 state_from_sasspec(spec)
6893 }
6894 })
6895 .transpose()
6896 .map_err(EstimationError::InvalidInput)?,
6897 y: y.to_owned(),
6898 weights: weights.to_owned(),
6899 design: fit_design.clone(),
6900 designzeroed: designzeroed.clone(),
6901 offset: offset.to_owned(),
6902 bounded_terms: bounded_terms.clone(),
6903 };
6904 let blockspec = ParameterBlockSpec {
6905 name: "eta".to_string(),
6906 design: DesignMatrix::Dense(gam_linalg::matrix::DenseDesignMatrix::from(designzeroed)),
6907 offset: offset.to_owned(),
6908 penalties: fit_penalties
6909 .iter()
6910 .map(|ps| match ps {
6911 PenaltySpec::Block {
6912 local, col_range, ..
6913 } => PenaltyMatrix::Blockwise {
6914 local: local.clone(),
6915 col_range: col_range.clone(),
6916 total_dim: design.design.ncols(),
6917 },
6918 PenaltySpec::Dense(m) | PenaltySpec::DenseWithMean { matrix: m, .. } => {
6919 PenaltyMatrix::Dense(m.clone())
6920 }
6921 })
6922 .collect(),
6923 nullspace_dims: design.nullspace_dims.clone(),
6924 initial_log_lambdas,
6925 initial_beta: Some(initial_beta),
6926 gauge_priority: 100,
6927 jacobian_callback: Some(Arc::new(BoundedEffectiveJacobian {
6933 design: fit_design.clone(),
6934 bounded_terms: bounded_terms.clone(),
6935 })),
6936 stacked_design: None,
6937 stacked_offset: None,
6938 };
6939 let fit = fit_custom_family(
6940 &family_adapter,
6941 &[blockspec],
6942 &BlockwiseFitOptions {
6943 inner_max_cycles: options.max_iter,
6944 inner_tol: options.tol,
6945 outer_max_iter: options.max_iter,
6946 outer_tol: options.tol,
6947 compute_covariance: false,
6957 ..BlockwiseFitOptions::default()
6958 },
6959 )
6960 .map_err(EstimationError::CustomFamily)?;
6961
6962 let latent_beta = fit.block_states[0].beta.clone();
6963 let (beta_user_internal, jac_diag) = family_adapter
6964 .user_beta_and_jacobian(&latent_beta)
6965 .map_err(EstimationError::InvalidInput)?;
6966 let beta_user = conditioning.backtransform_beta(&beta_user_internal);
6967
6968 let (eta_state, h_data, _, _) = family_adapter
6969 .evaluation_from_latent(&latent_beta)
6970 .map_err(EstimationError::InvalidInput)?;
6971 let p_fit = fit_design.ncols();
6972 let mut s_lambda_internal = Array2::<f64>::zeros((p_fit, p_fit));
6973 for (k, penalty) in fit_penalties.iter().enumerate() {
6974 match penalty {
6975 PenaltySpec::Block {
6976 local, col_range, ..
6977 } => {
6978 s_lambda_internal
6979 .slice_mut(ndarray::s![col_range.clone(), col_range.clone()])
6980 .scaled_add(fit.lambdas[k], local);
6981 }
6982 PenaltySpec::Dense(m) | PenaltySpec::DenseWithMean { matrix: m, .. } => {
6983 s_lambda_internal.scaled_add(fit.lambdas[k], m);
6984 }
6985 }
6986 }
6987 let mut latent_precision = h_data.clone();
6988 latent_precision += &s_lambda_internal;
6989 let user_precision_internal =
6990 transform_bounded_latent_precision_to_user_internal(&latent_precision, &jac_diag)?;
6991 let penalized_hessian =
6992 conditioning.transform_penalized_hessian_to_original(&user_precision_internal);
6993
6994 let beta_covariance_unscaled = if options.compute_inference {
7021 Some(certified_bounded_posterior_covariance(
7022 &penalized_hessian,
7023 "bounded user-scale posterior precision",
7024 )?)
7025 } else {
7026 None
7027 };
7028 let latent_cov = if options.compute_inference {
7034 Some(certified_bounded_posterior_covariance(
7035 &latent_precision,
7036 "bounded latent posterior precision",
7037 )?)
7038 } else {
7039 None
7040 };
7041 let s_lambda_original = weighted_blockwise_penalty_sum(
7042 &design.penalties,
7043 fit.lambdas.as_slice().unwrap(),
7044 design.design.ncols(),
7045 );
7046 let penalty_term = beta_user.dot(&s_lambda_original.dot(&beta_user));
7047 let deviance = -2.0 * eta_state.log_likelihood;
7048 let (edf_by_block, penalty_block_trace, edf_total) = if let Some(cov) = latent_cov.as_ref() {
7049 exact_bounded_edf(&fit_penalties, &fit.lambdas, cov)?
7050 } else {
7051 (
7052 vec![0.0; fit_penalties.len()],
7053 vec![0.0; fit_penalties.len()],
7054 0.0,
7055 )
7056 };
7057
7058 let profiled_gaussian_standard_deviation = if matches!(
7070 resolved_likelihood_scale,
7071 gam_spec::ResolvedLikelihoodScale::ProfiledGaussian
7072 ) {
7073 let residual_dof = if options.compute_inference {
7074 y.len() as f64 - edf_total
7075 } else {
7076 y.len() as f64
7077 };
7078 if !(residual_dof.is_finite() && residual_dof > 0.0) {
7079 return Err(EstimationError::InvalidInput(format!(
7080 "bounded Gaussian residual degrees of freedom must be finite and positive, got n={} minus edf={edf_total} = {residual_dof}",
7081 y.len()
7082 )));
7083 }
7084 if !(deviance.is_finite() && deviance >= 0.0) {
7085 return Err(EstimationError::InvalidInput(format!(
7086 "bounded Gaussian deviance must be finite and non-negative, got {deviance}"
7087 )));
7088 }
7089 let variance = deviance / residual_dof;
7090 if !variance.is_finite() {
7091 return Err(EstimationError::InvalidInput(format!(
7092 "bounded Gaussian residual variance is not representable: {deviance}/{residual_dof}"
7093 )));
7094 }
7095 Some(variance.sqrt())
7096 } else {
7097 None
7098 };
7099 let dispersion = gam_solve::estimate::dispersion_from_likelihood(
7100 &glm_likelihood,
7101 profiled_gaussian_standard_deviation,
7102 )?;
7103 let standard_deviation = dispersion.phi().sqrt();
7104 let cov_scale = glm_likelihood
7105 .coefficient_covariance_scale(dispersion.phi())
7106 .map_err(|error| EstimationError::InvalidInput(error.to_string()))?;
7107 let beta_covariance = beta_covariance_unscaled.map(|mut cov| {
7113 if cov_scale != 1.0 {
7114 cov.mapv_inplace(|v| v * cov_scale);
7115 }
7116 cov
7117 });
7118 if let Some(covariance) = beta_covariance.as_ref()
7119 && covariance.iter().any(|value| !value.is_finite())
7120 {
7121 return Err(EstimationError::InvalidInput(
7122 "bounded coefficient covariance scaling produced a non-finite value".to_string(),
7123 ));
7124 }
7125 let beta_standard_errors = beta_covariance
7126 .as_ref()
7127 .map(gam_problem::se_from_covariance)
7128 .transpose()
7129 .map_err(|err| {
7130 EstimationError::InvalidInput(format!(
7131 "bounded coefficient covariance cannot produce standard errors: {err}"
7132 ))
7133 })?;
7134 let working_response = exact_standard_working_response(&eta_state)?;
7135
7136 let geometry = Some(gam_solve::estimate::FitGeometry {
7137 coefficient_gauge: gam_problem::gauge::Gauge::identity(&[beta_user.len()]),
7138 penalized_hessian: penalized_hessian.clone().into(),
7139 constrained_posterior: None,
7140 working: Some(gam_solve::estimate::WorkingGeometry {
7141 weights: eta_state.fisherweight.clone(),
7142 response: working_response,
7143 }),
7144 });
7145 let max_abs_eta = eta_state
7146 .eta
7147 .iter()
7148 .fold(0.0_f64, |acc, &v| acc.max(v.abs()));
7149 Ok(FittedTermCollection {
7150 fit: {
7151 let log_lambdas =
7152 checked_fit_log_lambdas(&fit.lambdas, "final fitted term collection")?;
7153 let inf = FitInference {
7154 edf_by_block,
7155 penalty_block_trace,
7156 edf_total,
7157 smoothing_correction: None,
7158 smoothing_correction_method: None,
7159 smoothing_correction_first_order: None,
7160 smoothing_correction_method_first_order: None,
7161 penalized_hessian: penalized_hessian.clone().into(),
7164 reparam_qs: None,
7165 dispersion,
7166 beta_covariance: beta_covariance
7167 .clone()
7168 .map(gam_problem::dispersion_cov::PhiScaledCovariance::from),
7169 beta_standard_errors,
7170 beta_covariance_corrected: None,
7171 beta_standard_errors_corrected: None,
7172 beta_covariance_frequentist: None,
7173 coefficient_influence: None,
7174 weighted_gram: None,
7175 bias_correction_beta: None,
7176 bias_correction_jacobian: None,
7177 };
7178 let covariance_conditional = beta_covariance;
7179 let pirls_status_val = gam_solve::pirls::PirlsStatus::Converged;
7182 UnifiedFitResult::try_from_parts(UnifiedFitResultParts {
7183 blocks: vec![gam_solve::estimate::FittedBlock {
7184 beta: beta_user.clone(),
7185 role: gam_problem::BlockRole::Mean,
7186 edf: edf_total,
7187 lambdas: fit.lambdas.clone(),
7188 }],
7189 log_lambdas,
7190 lambdas: fit.lambdas,
7191 likelihood_scale: glm_likelihood.scale,
7192 likelihood_family: Some(glm_likelihood.spec),
7193 log_likelihood_normalization: gam_spec::LogLikelihoodNormalization::UserProvided,
7194 log_likelihood: eta_state.log_likelihood,
7195 deviance,
7196 reml_score: fit.penalized_objective,
7197 stable_penalty_term: penalty_term,
7198 penalized_objective: fit.penalized_objective,
7199 used_device: false,
7200 outer_iterations: fit.outer_iterations,
7201 outer_converged: true,
7203 outer_gradient_norm: fit.outer_gradient_norm,
7204 standard_deviation,
7205 covariance_conditional,
7206 covariance_corrected: None,
7207 inference: Some(inf),
7208 fitted_link: gam_solve::estimate::FittedLinkState::Standard(None),
7209 geometry,
7210 block_states: Vec::new(),
7211 pirls_status: pirls_status_val,
7212 max_abs_eta,
7213 constraint_kkt: None,
7214 artifacts: gam_solve::estimate::FitArtifacts {
7215 pirls: None,
7216 ..Default::default()
7217 },
7218 inner_cycles: 0,
7219 })?
7220 },
7221 design: design.clone(),
7222 adaptive_diagnostics: None,
7223 })
7224}
7225
7226fn enforce_term_constraint_feasibility(
7227 design: &TermCollectionDesign,
7228 fit: &UnifiedFitResult,
7229) -> Result<(), EstimationError> {
7230 const CONSTRAINT_FEASIBILITY_RAW_TOL: f64 = 1e-7;
7244 let tol = CONSTRAINT_FEASIBILITY_RAW_TOL;
7245 let smooth_start = design
7246 .design
7247 .ncols()
7248 .saturating_sub(design.smooth.total_smooth_cols());
7249 let mut violations: Vec<String> = Vec::new();
7250 for term in &design.smooth.terms {
7251 let gr = (smooth_start + term.coeff_range.start)..(smooth_start + term.coeff_range.end);
7252 let beta_local = fit.beta.slice(s![gr.clone()]).to_owned();
7253 if let Some(lb) = term.lower_bounds_local.as_ref() {
7254 let mut worst = 0.0_f64;
7255 let mut worst_idx = 0usize;
7256 for i in 0..lb.len().min(beta_local.len()) {
7257 if lb[i].is_finite() {
7258 let viol = (lb[i] - beta_local[i]).max(0.0);
7259 if viol > worst {
7260 worst = viol;
7261 worst_idx = i;
7262 }
7263 }
7264 }
7265 if worst > tol {
7266 violations.push(format!(
7267 "term='{}' kind=lower-bound maxviolation={:.3e} coeff_index={}",
7268 term.name, worst, worst_idx
7269 ));
7270 }
7271 }
7272 if let Some(lin) = term.linear_constraints_local.as_ref() {
7273 let mut worst = 0.0_f64;
7274 let mut worstrow = 0usize;
7275 for i in 0..lin.a.nrows() {
7276 let norm = lin.a.row(i).dot(&lin.a.row(i)).sqrt();
7277 let inv = if norm > 0.0 { 1.0 / norm } else { 0.0 };
7278 let s = (lin.a.row(i).dot(&beta_local) - lin.b[i]) * inv;
7279 let viol = (-s).max(0.0);
7280 if viol > worst {
7281 worst = viol;
7282 worstrow = i;
7283 }
7284 }
7285 if worst > tol {
7286 violations.push(format!(
7287 "term='{}' kind=linear-inequality maxviolation={:.3e} row={}",
7288 term.name, worst, worstrow
7289 ));
7290 }
7291 }
7292 }
7293
7294 if !violations.is_empty() {
7295 let mut msg = format!(
7296 "constraint violation after fit ({} violating term constraints): {}",
7297 violations.len(),
7298 violations.join(" | ")
7299 );
7300 if let Some(kkt) = fit.constraint_kkt.as_ref() {
7301 msg.push_str(&format!(
7302 "; KKT[primal={:.3e}, dual={:.3e}, comp={:.3e}, stat={:.3e}]",
7303 kkt.primal_feasibility, kkt.dual_feasibility, kkt.complementarity, kkt.stationarity
7304 ));
7305 }
7306 return Err(EstimationError::ParameterConstraintViolation(msg));
7307 }
7308 Ok(())
7309}
7310
7311fn stratified_spatial_subsample(
7312 data: ArrayView2<'_, f64>,
7313 spec: &TermCollectionSpec,
7314 target_size: usize,
7315) -> Vec<usize> {
7316 use rand::SeedableRng;
7317 use rand::rngs::StdRng;
7318 use rand::seq::SliceRandom;
7319
7320 let n = data.nrows();
7321 if n <= target_size {
7322 return (0..n).collect();
7323 }
7324
7325 let spatial_cols: Option<Vec<usize>> =
7326 spec.smooth_terms.iter().find_map(|term| match &term.basis {
7327 SmoothBasisSpec::ThinPlate { feature_cols, .. }
7328 | SmoothBasisSpec::Matern { feature_cols, .. }
7329 | SmoothBasisSpec::Duchon { feature_cols, .. } => {
7330 if !feature_cols.is_empty() {
7331 Some(feature_cols.clone())
7332 } else {
7333 None
7334 }
7335 }
7336 _ => None,
7337 });
7338
7339 let cols = match spatial_cols {
7340 Some(c) if !c.is_empty() => c,
7341 _ => {
7342 let mut rng = StdRng::seed_from_u64(spatial_subsample_seed(data, &[], target_size));
7343 let mut indices: Vec<usize> = (0..n).collect();
7344 indices.shuffle(&mut rng);
7345 indices.truncate(target_size);
7346 indices.sort_unstable();
7347 return indices;
7348 }
7349 };
7350 let mut rng = StdRng::seed_from_u64(spatial_subsample_seed(data, &cols, target_size));
7351
7352 let d = cols.len();
7353 let mut mins = vec![f64::INFINITY; d];
7354 let mut maxs = vec![f64::NEG_INFINITY; d];
7355 for i in 0..n {
7356 for (ax, &col) in cols.iter().enumerate() {
7357 let v = data[[i, col]];
7358 if v < mins[ax] {
7359 mins[ax] = v;
7360 }
7361 if v > maxs[ax] {
7362 maxs[ax] = v;
7363 }
7364 }
7365 }
7366
7367 const TARGET_POINTS_PER_CELL: usize = 5;
7371 let total_cells_target = (target_size / TARGET_POINTS_PER_CELL).max(1);
7372 let cells_per_axis = ((total_cells_target as f64).powf(1.0 / d as f64)).ceil() as usize;
7373 let cells_per_axis = cells_per_axis.max(1);
7374
7375 let mut cell_members: std::collections::HashMap<Vec<usize>, Vec<usize>> =
7376 std::collections::HashMap::new();
7377 for i in 0..n {
7378 let mut cell_key = Vec::with_capacity(d);
7379 for (ax, &col) in cols.iter().enumerate() {
7380 let range = maxs[ax] - mins[ax];
7381 let cell = if range <= 0.0 {
7382 0
7383 } else {
7384 let frac = (data[[i, col]] - mins[ax]) / range;
7385 (frac * cells_per_axis as f64).floor() as usize
7386 };
7387 cell_key.push(cell.min(cells_per_axis - 1));
7388 }
7389 cell_members.entry(cell_key).or_default().push(i);
7390 }
7391
7392 let mut selected: Vec<usize> = Vec::with_capacity(target_size);
7393 let mut remaining_budget = target_size;
7394 let mut remaining_population = n;
7395
7396 let mut cells: Vec<(Vec<usize>, Vec<usize>)> = cell_members.into_iter().collect();
7397 cells.sort_by(|a, b| a.0.cmp(&b.0));
7398
7399 for (_, members) in &mut cells {
7400 if remaining_budget == 0 {
7401 break;
7402 }
7403 let alloc = ((members.len() as f64 / remaining_population as f64) * remaining_budget as f64)
7404 .round() as usize;
7405 let alloc = alloc.max(1).min(members.len()).min(remaining_budget);
7406 members.shuffle(&mut rng);
7407 selected.extend_from_slice(&members[..alloc]);
7408 remaining_budget = remaining_budget.saturating_sub(alloc);
7409 remaining_population = remaining_population.saturating_sub(members.len());
7410 }
7411
7412 if selected.len() > target_size {
7413 selected.shuffle(&mut rng);
7414 selected.truncate(target_size);
7415 }
7416
7417 selected.sort_unstable();
7418 selected
7419}
7420
7421fn spatial_subsample_seed(
7422 data: ArrayView2<'_, f64>,
7423 spatial_cols: &[usize],
7424 target_size: usize,
7425) -> u64 {
7426 let mut state = 0x5350_4154_4941_4C53_u64;
7427 spatial_seed_mix(&mut state, data.nrows() as u64);
7428 spatial_seed_mix(&mut state, data.ncols() as u64);
7429 spatial_seed_mix(&mut state, target_size as u64);
7430 spatial_seed_mix(&mut state, spatial_cols.len() as u64);
7431 for &col in spatial_cols {
7432 spatial_seed_mix(&mut state, col as u64);
7433 }
7434
7435 if data.nrows() > 0 {
7436 let mid = data.nrows() / 2;
7437 let last = data.nrows() - 1;
7438 for &row in &[0usize, mid, last] {
7439 for &col in spatial_cols {
7440 let value = data[[row, col]];
7441 spatial_seed_mix(&mut state, value.to_bits());
7442 }
7443 }
7444 }
7445 state
7446}
7447
7448#[inline]
7449fn spatial_seed_mix(state: &mut u64, value: u64) {
7450 let mut s = value.wrapping_add(*state);
7453 let z = gam_linalg::utils::splitmix64(&mut s);
7454 *state ^= z;
7455 *state = (*state).rotate_left(27).wrapping_mul(0x3C79_AC49_2BA7_B653);
7456}
7457
7458fn sampled_rows(data: ArrayView2<'_, f64>, indices: &[usize]) -> Array2<f64> {
7459 let mut sampled = Array2::<f64>::zeros((indices.len(), data.ncols()));
7460 for (new_row, &orig_row) in indices.iter().enumerate() {
7461 sampled.row_mut(new_row).assign(&data.row(orig_row));
7462 }
7463 sampled
7464}
7465
7466fn spatial_term_user_centers(term: &SmoothTermSpec) -> Option<ArrayView2<'_, f64>> {
7467 match spatial_term_center_strategy(term) {
7468 Some(CenterStrategy::UserProvided(centers)) => Some(centers.view()),
7469 _ => None,
7470 }
7471}
7472
7473fn finite_centered_axis_contrasts(values: &[f64], expected_dim: usize) -> Option<Vec<f64>> {
7474 if values.len() != expected_dim || expected_dim <= 1 {
7475 return None;
7476 }
7477 if values.iter().any(|value| !value.is_finite()) {
7478 return None;
7479 }
7480 Some(center_aniso_log_scales(values))
7481}
7482
7483fn blended_pilot_axis_contrasts(
7484 pilot_data: ArrayView2<'_, f64>,
7485 term: &SmoothTermSpec,
7486 centers: ArrayView2<'_, f64>,
7487) -> Result<Option<Vec<f64>>, BasisError> {
7488 let d = centers.ncols();
7489 if d <= 1 {
7490 return Ok(None);
7491 }
7492 let center_eta = initial_aniso_contrasts(centers);
7493 let standardized_data = standardized_spatial_term_data(pilot_data, term)?;
7494 let data_eta = finite_centered_axis_contrasts(
7495 &initial_aniso_contrasts(standardized_data.view()),
7496 d,
7497 );
7498 let Some(center_eta) = finite_centered_axis_contrasts(¢er_eta, d) else {
7499 return Ok(None);
7500 };
7501 let blended = match data_eta {
7502 Some(data_eta) => center_eta
7503 .iter()
7504 .zip(data_eta.iter())
7505 .map(|(&from_centers, &from_data)| 0.5 * (from_centers + from_data))
7506 .collect::<Vec<_>>(),
7507 None => center_eta,
7508 };
7509 Ok(finite_centered_axis_contrasts(&blended, d))
7510}
7511
7512fn apply_pilot_spatial_psi_reseed(
7513 pilot_data: ArrayView2<'_, f64>,
7514 spec: &TermCollectionSpec,
7515 spatial_terms: &[usize],
7516 kappa_options: &SpatialLengthScaleOptimizationOptions,
7517) -> Result<TermCollectionSpec, EstimationError> {
7518 let dims_per_term = spatial_dims_per_term(spec, spatial_terms);
7519 let use_aniso = has_aniso_terms(spec, spatial_terms);
7520 let log_kappa0 = if use_aniso {
7521 SpatialLogKappaCoords::from_length_scales_aniso(spec, spatial_terms, kappa_options)
7522 } else {
7523 SpatialLogKappaCoords::from_length_scales(spec, spatial_terms, kappa_options)
7524 };
7525 let log_kappa0 = log_kappa0
7526 .reseed_from_data(pilot_data, spec, spatial_terms, kappa_options)
7527 .map_err(EstimationError::BasisError)?;
7528 let log_kappa_lower = if use_aniso {
7529 SpatialLogKappaCoords::lower_bounds_aniso_from_data(
7530 pilot_data,
7531 spec,
7532 spatial_terms,
7533 &dims_per_term,
7534 kappa_options,
7535 )
7536 } else {
7537 SpatialLogKappaCoords::lower_bounds_from_data(
7538 pilot_data,
7539 spec,
7540 spatial_terms,
7541 kappa_options,
7542 )
7543 }
7544 .map_err(EstimationError::BasisError)?;
7545 let log_kappa_upper = if use_aniso {
7546 SpatialLogKappaCoords::upper_bounds_aniso_from_data(
7547 pilot_data,
7548 spec,
7549 spatial_terms,
7550 &dims_per_term,
7551 kappa_options,
7552 )
7553 } else {
7554 SpatialLogKappaCoords::upper_bounds_from_data(
7555 pilot_data,
7556 spec,
7557 spatial_terms,
7558 kappa_options,
7559 )
7560 }
7561 .map_err(EstimationError::BasisError)?;
7562 log_kappa0
7563 .clamp_to_bounds(&log_kappa_lower, &log_kappa_upper)
7564 .apply_tospec(spec, spatial_terms)
7565}
7566
7567pub(crate) fn apply_spatial_anisotropy_pilot_initializer(
7568 data: ArrayView2<'_, f64>,
7569 spec: &mut TermCollectionSpec,
7570 spatial_terms: &[usize],
7571 target_size: usize,
7572 kappa_options: &SpatialLengthScaleOptimizationOptions,
7573) -> Result<usize, EstimationError> {
7574 if target_size == 0 || data.nrows() <= target_size.saturating_mul(2) || spatial_terms.is_empty()
7575 {
7576 return Ok(0);
7577 }
7578 if !has_aniso_terms(spec, spatial_terms) {
7579 return Ok(0);
7580 }
7581 let indices = stratified_spatial_subsample(data, spec, target_size);
7582 let pilot_data = sampled_rows(data, &indices);
7583 let mut working = spec.clone();
7584 let mut updated_terms = 0usize;
7585 const GEOMETRY_UPDATES: usize = 2;
7586
7587 for pass in 0..GEOMETRY_UPDATES {
7588 let planned_terms = plan_joint_spatial_centers_for_term_blocks(
7589 pilot_data.view(),
7590 &[working.smooth_terms.clone()],
7591 )
7592 .and_then(|mut blocks| {
7593 blocks.pop().ok_or_else(|| {
7594 BasisError::InvalidInput(
7595 "pilot geometry initializer produced no smooth-term block".to_string(),
7596 )
7597 })
7598 })
7599 .map_err(EstimationError::BasisError)?;
7600
7601 for &term_idx in spatial_terms {
7602 let Some(current_eta) = get_spatial_aniso_log_scales(&working, term_idx) else {
7603 continue;
7604 };
7605 let Some(d) = get_spatial_feature_dim(&working, term_idx) else {
7606 continue;
7607 };
7608 if d <= 1 || current_eta.len() != d {
7609 continue;
7610 }
7611 let Some(planned_term) = planned_terms.get(term_idx) else {
7612 continue;
7613 };
7614 let Some(centers) = spatial_term_user_centers(planned_term) else {
7615 continue;
7616 };
7617 let Some(eta) = blended_pilot_axis_contrasts(
7618 pilot_data.view(),
7619 planned_term,
7620 centers,
7621 )
7622 .map_err(EstimationError::BasisError)?
7623 else {
7624 continue;
7625 };
7626 set_spatial_aniso_log_scales(&mut working, term_idx, eta)?;
7627 updated_terms += usize::from(pass == 0);
7628 }
7629
7630 working = apply_pilot_spatial_psi_reseed(
7631 pilot_data.view(),
7632 &working,
7633 spatial_terms,
7634 kappa_options,
7635 )?;
7636 }
7637
7638 if updated_terms > 0 {
7639 log::info!(
7640 "[spatial-kappa] initialized anisotropy from {}-row pilot geometry for {} spatial term(s); proceeding to full-data optimization",
7641 indices.len(),
7642 updated_terms
7643 );
7644 *spec = working;
7645 }
7646 Ok(updated_terms)
7647}
7648
7649pub(crate) fn spatial_length_scale_term_indices(spec: &TermCollectionSpec) -> Vec<usize> {
7650 spec.smooth_terms
7651 .iter()
7652 .enumerate()
7653 .filter_map(|(idx, _)| spatial_term_supports_hyper_optimization(spec, idx).then_some(idx))
7654 .collect()
7655}
7656
7657fn fit_score(fit: &UnifiedFitResult) -> f64 {
7669 if fit.reml_score.is_finite() {
7670 return fit.reml_score;
7671 }
7672 let score = 0.5 * fit.deviance + 0.5 * fit.stable_penalty_term;
7673 if score.is_finite() {
7674 score
7675 } else {
7676 f64::INFINITY
7677 }
7678}
7679
7680fn is_recoverable_trial_point_error(err: &EstimationError) -> bool {
7702 matches!(err, EstimationError::BasisError(_))
7703 || err.is_inner_solve_retreat()
7704 || is_recoverable_fit_inference_finiteness_error(err)
7705}
7706
7707fn is_recoverable_fit_inference_finiteness_error(err: &EstimationError) -> bool {
7708 let EstimationError::InvalidInput(message) = err else {
7709 return false;
7710 };
7711
7712 message.contains("must be finite")
7713 && [
7714 "fit_result.beta_covariance_frequentist",
7715 "fit_result.coefficient_influence",
7716 "fit_result.weighted_gram",
7717 ]
7718 .iter()
7719 .any(|field| message.contains(field))
7720}
7721
7722#[cfg(test)]
7723mod spatial_trial_recovery_tests {
7724 use super::*;
7725
7726 #[test]
7727 fn nonfinite_frequentist_covariance_is_recoverable_trial_point() {
7728 let err = EstimationError::InvalidInput(
7729 "fit_result.beta_covariance_frequentist[0] must be finite, got NaN".to_string(),
7730 );
7731
7732 assert!(
7733 is_recoverable_trial_point_error(&err),
7734 "singular trial-point curvature should make spatial κ retreat, not abort"
7735 );
7736 }
7737
7738 #[test]
7739 fn arbitrary_invalid_input_remains_fatal_trial_point_error() {
7740 let err = EstimationError::InvalidInput("outer rho bounds are invalid".to_string());
7741
7742 assert!(
7743 !is_recoverable_trial_point_error(&err),
7744 "the spatial κ recovery gate must not mask unrelated invalid inputs"
7745 );
7746 }
7747}
7748
7749fn require_successful_spatial_optimization_result<T>(
7750 initial_score: f64,
7751 result: Result<Option<(T, f64)>, EstimationError>,
7752) -> Result<T, EstimationError> {
7753 match result {
7754 Ok(Some((value, exact_score))) => {
7755 const SCORE_DRIFT_ABS_TOL: f64 = 1e-6;
7764 const SCORE_DRIFT_REL_TOL: f64 = 1e-8;
7765 let tol = SCORE_DRIFT_ABS_TOL.max(initial_score.abs() * SCORE_DRIFT_REL_TOL);
7766 if exact_score <= initial_score + tol {
7767 Ok(value)
7768 } else {
7769 Err(EstimationError::RemlOptimizationFailed(format!(
7770 "spatial kappa optimization made REML score worse ({initial_score:.6e} -> {exact_score:.6e})"
7771 )))
7772 }
7773 }
7774 Ok(None) => Err(EstimationError::RemlOptimizationFailed(
7775 "spatial kappa optimization is unavailable for one or more eligible spatial terms"
7776 .to_string(),
7777 )),
7778 Err(err) => Err(EstimationError::RemlOptimizationFailed(format!(
7779 "spatial kappa optimization failed: {err}"
7780 ))),
7781 }
7782}
7783
7784fn external_opts_for_design(
7785 family: &LikelihoodSpec,
7786 design: &TermCollectionDesign,
7787 options: &FitOptions,
7788) -> ExternalOptimOptions {
7789 ExternalOptimOptions {
7790 family: family.clone(),
7791 latent_cloglog: options.latent_cloglog,
7792 mixture_link: options.mixture_link.clone(),
7793 optimize_mixture: options.optimize_mixture,
7794 sas_link: options.sas_link,
7795 optimize_sas: options.optimize_sas,
7796 compute_inference: options.compute_inference,
7797 skip_rho_posterior_inference: options.skip_rho_posterior_inference,
7798 max_iter: options.max_iter,
7799 tol: options.tol,
7800 nullspace_dims: design.nullspace_dims.clone(),
7801 linear_constraints: design.linear_constraints.clone(),
7802 firth_bias_reduction: Some(options.firth_bias_reduction),
7803 penalty_shrinkage_floor: options.penalty_shrinkage_floor,
7804 rho_prior: options.rho_prior.clone(),
7805 kronecker_penalty_system: design.kronecker_penalty_system(),
7808 kronecker_factored: design
7809 .smooth
7810 .terms
7811 .iter()
7812 .find_map(|t| t.kronecker_factored.clone()),
7813 persist_warm_start_disk: options.persist_warm_start_disk,
7814 }
7815}
7816
7817fn evaluate_joint_reml_outer_eval_at_theta(
7825 evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>,
7826 design: &TermCollectionDesign,
7827 theta: &Array1<f64>,
7828 rho_dim: usize,
7829 hyper_dirs: Vec<gam_solve::estimate::reml::DirectionalHyperParam>,
7830 warm_start_beta: Option<ArrayView1<'_, f64>>,
7831 order: gam_solve::rho_optimizer::OuterEvalOrder,
7832 design_revision: Option<u64>,
7833) -> Result<(f64, Array1<f64>, gam_problem::HessianValue), EstimationError> {
7834 evaluator.evaluate_with_order(
7835 &design.design,
7836 &design.penalties,
7837 &design.nullspace_dims,
7838 design.linear_constraints.clone(),
7839 theta,
7840 rho_dim,
7841 hyper_dirs,
7842 warm_start_beta,
7843 "evaluate_joint_reml_outer_eval_at_theta",
7844 order,
7845 design_revision,
7846 )
7847}
7848
7849fn evaluate_joint_reml_efs_at_theta(
7850 evaluator: &mut gam_solve::estimate::ExternalJointHyperEvaluator<'_>,
7851 design: &TermCollectionDesign,
7852 theta: &Array1<f64>,
7853 rho_dim: usize,
7854 hyper_dirs: Vec<gam_solve::estimate::reml::DirectionalHyperParam>,
7855 warm_start_beta: Option<ArrayView1<'_, f64>>,
7856 design_revision: Option<u64>,
7857) -> Result<gam_problem::EfsEval, EstimationError> {
7858 evaluator.evaluate_efs(
7859 &design.design,
7860 &design.penalties,
7861 &design.nullspace_dims,
7862 design.linear_constraints.clone(),
7863 theta,
7864 rho_dim,
7865 hyper_dirs,
7866 warm_start_beta,
7867 "evaluate_joint_reml_efs_at_theta",
7868 design_revision,
7869 )
7870}
7871
7872fn exact_joint_spatial_outer_hessian_available(
7873 family: &LikelihoodSpec,
7874 design: &TermCollectionDesign,
7875) -> bool {
7876 let family_supported = match &family.response {
7899 ResponseFamily::Gaussian
7900 | ResponseFamily::Binomial
7901 | ResponseFamily::Poisson
7902 | ResponseFamily::Tweedie { .. }
7903 | ResponseFamily::NegativeBinomial { .. }
7904 | ResponseFamily::Beta { .. }
7905 | ResponseFamily::Gamma
7906 | ResponseFamily::RoystonParmar => true,
7907 };
7908 family_supported && design.design.ncols() > 0
7911}
7912
7913fn try_build_spatial_term_log_kappa_derivativeinfo(
7914 data: ArrayView2<'_, f64>,
7915 resolvedspec: &TermCollectionSpec,
7916 design: &TermCollectionDesign,
7917 term_idx: usize,
7918) -> Result<Option<SpatialPsiDerivative>, EstimationError> {
7919 let Some((
7920 global_range,
7921 total_p,
7922 x_psi_local,
7923 s_psi_local_check,
7924 x_psi_psi_local,
7925 s_psi_psi_local,
7926 s_psi_components_local,
7927 s_psi_psi_components_local,
7928 implicit_operator,
7929 )) = try_build_spatial_term_log_kappa_derivative(data, resolvedspec, design, term_idx)?
7930 else {
7931 return Ok(None);
7932 };
7933 let Some(penalty_range) = design
7934 .smooth_term_penalty_range(term_idx)
7935 .map_err(EstimationError::InvalidInput)?
7936 else {
7937 return Ok(None);
7938 };
7939 let penalty_start = penalty_range.start;
7940 if s_psi_components_local.is_empty() || s_psi_psi_components_local.is_empty() {
7941 return Ok(None);
7942 }
7943 if s_psi_components_local.len() != s_psi_psi_components_local.len() {
7944 return Ok(None);
7945 }
7946 let penalty_indices = (0..s_psi_components_local.len())
7947 .map(|j| penalty_start + j)
7948 .collect::<Vec<_>>();
7949 let penalty_index = penalty_indices[0];
7950 if s_psi_local_check.nrows() == 0 || s_psi_psi_local.nrows() == 0 {
7951 return Ok(None);
7952 }
7953 Ok(Some(SpatialPsiDerivative {
7954 penalty_index,
7955 penalty_indices,
7956 global_range,
7957 total_p,
7958 x_psi_local,
7959 s_psi_components_local,
7960 x_psi_psi_local,
7961 s_psi_psi_components_local,
7962 aniso_group_id: None,
7963 aniso_cross_designs: None,
7964 aniso_cross_penalty_provider: None,
7965 implicit_operator,
7966 implicit_axis: 0,
7967 }))
7968}
7969
7970pub(crate) fn try_build_spatial_log_kappa_derivativeinfo_list(
7971 data: ArrayView2<'_, f64>,
7972 resolvedspec: &TermCollectionSpec,
7973 design: &TermCollectionDesign,
7974 spatial_terms: &[usize],
7975) -> Result<Option<Vec<SpatialPsiDerivative>>, EstimationError> {
7976 let mut out = Vec::new();
7977 let mut aniso_gid = 0usize;
7978 for &term_idx in spatial_terms {
7979 if spatial_term_uses_per_axis_psi(resolvedspec, term_idx) {
7980 if let Some(entries) = try_build_spatial_term_log_kappa_aniso_derivativeinfos(
7981 data,
7982 resolvedspec,
7983 design,
7984 term_idx,
7985 aniso_gid,
7986 )? {
7987 aniso_gid += 1;
7988 out.extend(entries);
7989 continue;
7990 } else {
7991 return Ok(None);
7992 }
7993 }
7994 let Some(info) =
7995 try_build_spatial_term_log_kappa_derivativeinfo(data, resolvedspec, design, term_idx)?
7996 else {
7997 return Ok(None);
7998 };
7999 out.push(info);
8000 }
8001 Ok(Some(out))
8002}
8003
8004fn try_build_spatial_term_log_kappa_aniso_derivativeinfos(
8006 data: ArrayView2<'_, f64>,
8007 resolvedspec: &TermCollectionSpec,
8008 design: &TermCollectionDesign,
8009 term_idx: usize,
8010 aniso_group_id: usize,
8011) -> Result<Option<Vec<SpatialPsiDerivative>>, EstimationError> {
8012 let Some(smooth_term) = design.smooth.terms.get(term_idx) else {
8013 return Ok(None);
8014 };
8015 let Some(termspec) = resolvedspec.smooth_terms.get(term_idx) else {
8016 return Ok(None);
8017 };
8018 let mut aniso_result = match &termspec.basis {
8019 SmoothBasisSpec::Sphere { .. } => return Ok(None),
8020 SmoothBasisSpec::Matern {
8021 feature_cols,
8022 spec,
8023 input_scale,
8024 } => {
8025 let mut x = select_columns(data, feature_cols).map_err(EstimationError::from)?;
8026 let mut spec_operator = spec.clone();
8027 if let Some(scale) = input_scale {
8028 scale.standardize(&mut x);
8029 let length_scale = spec.length_scale.resolved().ok_or_else(|| {
8030 EstimationError::InvalidInput(
8031 "anisotropic Matérn Auto length_scale reached derivative construction \
8032 unresolved"
8033 .to_string(),
8034 )
8035 })?;
8036 spec_operator
8037 .length_scale
8038 .set_resolved(scale.to_standardized_units(length_scale));
8039 }
8040 spec_operator.double_penalty = false;
8049 build_matern_basis_log_kappa_aniso_derivatives(x.view(), &spec_operator)
8050 .map_err(EstimationError::from)?
8051 }
8052 SmoothBasisSpec::MeasureJet {
8058 feature_cols,
8059 spec,
8060 input_scale,
8061 } => {
8062 let mut x = select_columns(data, feature_cols).map_err(EstimationError::from)?;
8063 if let Some(scale) = input_scale {
8064 scale.standardize(&mut x);
8065 }
8066 build_measure_jet_basis_psi_derivatives(x.view(), spec)
8067 .map_err(EstimationError::from)?
8068 }
8069 _ => return Ok(None),
8070 };
8071 let d = if let Some(ref op) = aniso_result.implicit_operator {
8074 op.n_axes()
8075 } else if !aniso_result.design_first.is_empty() {
8076 aniso_result.design_first.len()
8077 } else {
8078 0
8079 };
8080 if d == 0 {
8081 return Ok(None);
8082 }
8083 let Some(penalty_range) = design
8084 .smooth_term_penalty_range(term_idx)
8085 .map_err(EstimationError::InvalidInput)?
8086 else {
8087 return Ok(None);
8088 };
8089 let penalty_start = penalty_range.start;
8090 let p_total = design.design.ncols();
8091 let smooth_start = p_total.saturating_sub(design.smooth.total_smooth_cols());
8092 let global_range = (smooth_start + smooth_term.coeff_range.start)
8093 ..(smooth_start + smooth_term.coeff_range.end);
8094 let num_penalties = aniso_result.penalties_first[0].len();
8095 let penalty_indices: Vec<usize> = (0..num_penalties).map(|j| penalty_start + j).collect();
8096 let penalties_cross_provider = aniso_result.penalties_cross_provider.clone();
8097
8098 let use_implicit_design = aniso_result.design_first.is_empty();
8102 let implicit_op_arc = aniso_result
8103 .implicit_operator
8104 .as_ref()
8105 .map(|op| std::sync::Arc::new(op.clone()));
8106
8107 let mut entries = Vec::with_capacity(d);
8108 for a in 0..d {
8109 let (x_psi_local, x_psi_psi_local) = if use_implicit_design {
8110 (Array2::<f64>::zeros((0, 0)), Array2::<f64>::zeros((0, 0)))
8116 } else {
8117 let x_first = std::mem::take(&mut aniso_result.design_first[a]);
8122 let x_second = std::mem::take(&mut aniso_result.design_second_diag[a]);
8123 if x_first.ncols() != smooth_term.coeff_range.len() {
8124 return Ok(None);
8125 }
8126 (x_first, x_second)
8127 };
8128 let s_psi_components = std::mem::take(&mut aniso_result.penalties_first[a]);
8129 let s_psi_psi_components = std::mem::take(&mut aniso_result.penalties_second_diag[a]);
8130 let cross_designs = if implicit_op_arc.is_some() {
8136 let mut cd = Vec::with_capacity(d - 1);
8137 for b in 0..d {
8138 if b == a {
8139 continue;
8140 }
8141 cd.push((b, Array2::<f64>::zeros((0, 0))));
8142 }
8143 cd
8144 } else if !aniso_result.design_second_cross.is_empty() {
8145 let mut cd = Vec::new();
8146 for (cross_idx, &(pa, pb)) in aniso_result.design_second_cross_pairs.iter().enumerate()
8147 {
8148 if pa == a {
8149 cd.push((pb, aniso_result.design_second_cross[cross_idx].clone()));
8150 } else if pb == a {
8151 cd.push((pa, aniso_result.design_second_cross[cross_idx].clone()));
8152 }
8153 }
8154 cd
8155 } else {
8156 Vec::new()
8157 };
8158 let cross_penalty_provider = if d > 1 {
8159 let penalties_cross_provider = penalties_cross_provider.clone();
8160 Some(std::sync::Arc::new(
8161 move |b_axis: usize| -> Result<Vec<Array2<f64>>, EstimationError> {
8162 if b_axis == a {
8163 return Ok(Vec::new());
8164 }
8165 let (axis_lo, axis_hi) = if a < b_axis { (a, b_axis) } else { (b_axis, a) };
8166 if let Some(provider) = penalties_cross_provider.as_ref() {
8167 provider
8168 .evaluate(axis_lo, axis_hi)
8169 .map_err(EstimationError::from)
8170 } else {
8171 Ok(Vec::new())
8175 }
8176 },
8177 )
8178 as std::sync::Arc<
8179 dyn Fn(usize) -> Result<Vec<Array2<f64>>, EstimationError>
8180 + Send
8181 + Sync
8182 + 'static,
8183 >)
8184 } else {
8185 None
8186 };
8187
8188 entries.push(SpatialPsiDerivative {
8189 penalty_index: penalty_indices[0],
8190 penalty_indices: penalty_indices.clone(),
8191 global_range: global_range.clone(),
8192 total_p: p_total,
8193 x_psi_local,
8194 s_psi_components_local: s_psi_components,
8195 x_psi_psi_local,
8196 s_psi_psi_components_local: s_psi_psi_components,
8197 aniso_group_id: Some(aniso_group_id),
8198 aniso_cross_designs: if cross_designs.is_empty() {
8199 None
8200 } else {
8201 Some(cross_designs)
8202 },
8203 aniso_cross_penalty_provider: cross_penalty_provider,
8204 implicit_operator: implicit_op_arc.clone(),
8205 implicit_axis: a,
8206 });
8207 }
8208 Ok(Some(entries))
8209}
8210
8211#[cfg(test)]
8212mod glm_eta_observation_fd_tests {
8213 use super::*;
8219 use ndarray::array;
8220
8221 fn one_obs_weight(
8222 spec: &LikelihoodSpec,
8223 y: f64,
8224 weight: f64,
8225 eta: f64,
8226 ) -> StandardFamilyObservationState {
8227 let yv = Array1::from_vec(vec![y]);
8228 let wv = Array1::from_vec(vec![weight]);
8229 let ev = Array1::from_vec(vec![eta]);
8230 evaluate_standard_familyobservations(spec.clone(), None, None, None, &yv, &wv, &ev)
8231 .expect("standard family observation state assembles")
8232 }
8233
8234 fn one_obs(spec: &LikelihoodSpec, y: f64, eta: f64) -> StandardFamilyObservationState {
8235 one_obs_weight(spec, y, 1.0, eta)
8236 }
8237
8238 fn one_obs_resolved(
8239 likelihood: &gam_spec::GlmLikelihoodSpec,
8240 y: f64,
8241 weight: f64,
8242 eta: f64,
8243 ) -> StandardFamilyObservationState {
8244 evaluate_resolved_standard_family_observations(
8245 likelihood,
8246 None,
8247 None,
8248 None,
8249 &array![y],
8250 &array![weight],
8251 &array![eta],
8252 )
8253 .expect("resolved standard family observation state assembles")
8254 }
8255
8256 #[test]
8257 fn bounded_gamma_and_tweedie_use_the_resolved_likelihood_scale() {
8258 let gamma_unit = gam_spec::GlmLikelihoodSpec {
8259 spec: LikelihoodSpec::gamma_log(),
8260 scale: gam_spec::LikelihoodScaleMetadata::FixedGammaShape { shape: 1.0 },
8261 };
8262 let gamma_scaled = gam_spec::GlmLikelihoodSpec {
8263 spec: LikelihoodSpec::gamma_log(),
8264 scale: gam_spec::LikelihoodScaleMetadata::FixedGammaShape { shape: 8.0 },
8265 };
8266 let unit = one_obs_resolved(&gamma_unit, 2.3, 0.7, 0.2);
8267 let scaled = one_obs_resolved(&gamma_scaled, 2.3, 0.7, 0.2);
8268 for (label, actual, base) in [
8269 ("Gamma score", scaled.score[0], unit.score[0]),
8270 (
8271 "Gamma Fisher weight",
8272 scaled.fisherweight[0],
8273 unit.fisherweight[0],
8274 ),
8275 (
8276 "Gamma observed Hessian",
8277 scaled.neghessian_eta[0],
8278 unit.neghessian_eta[0],
8279 ),
8280 (
8281 "Gamma Hessian derivative",
8282 scaled.neghessian_eta_derivative[0],
8283 unit.neghessian_eta_derivative[0],
8284 ),
8285 (
8286 "Gamma log likelihood",
8287 scaled.log_likelihood,
8288 unit.log_likelihood,
8289 ),
8290 ] {
8291 let expected = 8.0 * base;
8292 assert!(
8293 (actual - expected).abs() <= 32.0 * f64::EPSILON * expected.abs().max(1.0),
8294 "{label} scale mismatch: actual={actual}, expected={expected}"
8295 );
8296 }
8297
8298 let tweedie_unit = gam_spec::GlmLikelihoodSpec {
8299 spec: LikelihoodSpec::tweedie_log(1.5),
8300 scale: gam_spec::LikelihoodScaleMetadata::FixedDispersion { phi: 1.0 },
8301 };
8302 let tweedie_scaled = gam_spec::GlmLikelihoodSpec {
8303 spec: LikelihoodSpec::tweedie_log(1.5),
8304 scale: gam_spec::LikelihoodScaleMetadata::FixedDispersion { phi: 0.25 },
8305 };
8306 let unit = one_obs_resolved(&tweedie_unit, 1.7, 0.8, -0.1);
8307 let scaled = one_obs_resolved(&tweedie_scaled, 1.7, 0.8, -0.1);
8308 for (actual, base) in [
8309 (scaled.score[0], unit.score[0]),
8310 (scaled.fisherweight[0], unit.fisherweight[0]),
8311 (scaled.neghessian_eta[0], unit.neghessian_eta[0]),
8312 (
8313 scaled.neghessian_eta_derivative[0],
8314 unit.neghessian_eta_derivative[0],
8315 ),
8316 (scaled.log_likelihood, unit.log_likelihood),
8317 ] {
8318 let expected = 4.0 * base;
8319 assert!((actual - expected).abs() <= 32.0 * f64::EPSILON * expected.abs().max(1.0));
8320 }
8321 }
8322
8323 #[test]
8324 fn bounded_zero_rows_are_dormant_and_weight_preflight_is_atomic() {
8325 let likelihood = gam_spec::GlmLikelihoodSpec::canonical(LikelihoodSpec::poisson_log());
8326 let dormant = evaluate_resolved_standard_family_observations(
8327 &likelihood,
8328 None,
8329 None,
8330 None,
8331 &array![f64::NAN, 2.0],
8332 &array![0.0, 1.0],
8333 &array![f64::NAN, 0.2],
8334 )
8335 .expect("zero-weight response and predictor are dormant");
8336 assert_eq!(dormant.score[0], 0.0);
8337 assert_eq!(dormant.fisherweight[0], 0.0);
8338
8339 let error = evaluate_resolved_standard_family_observations(
8340 &likelihood,
8341 None,
8342 None,
8343 None,
8344 &array![f64::NAN, 2.0],
8345 &array![1.0, f64::NAN],
8346 &array![f64::NAN, 0.2],
8347 )
8348 .expect_err("later invalid weight must refuse before row evaluation");
8349 assert!(
8350 error.to_string().contains("row 2 has invalid prior weight"),
8351 "unexpected atomic preflight error: {error}"
8352 );
8353 }
8354
8355 fn check_fd(label: &str, spec: &LikelihoodSpec, y: f64, eta: f64) {
8356 let h = 1e-5;
8357 let s0 = one_obs(spec, y, eta);
8358 let sp = one_obs(spec, y, eta + h);
8359 let sm = one_obs(spec, y, eta - h);
8360
8361 let score_fd = (sp.log_likelihood - sm.log_likelihood) / (2.0 * h);
8363 let score = s0.score[0];
8364 assert!(
8365 (score - score_fd).abs() <= 1e-4 * (1.0 + score.abs()),
8366 "{label}: score {score} vs FD {score_fd}"
8367 );
8368
8369 let neghess_fd = -(sp.score[0] - sm.score[0]) / (2.0 * h);
8371 let neghess = s0.neghessian_eta[0];
8372 assert!(
8373 (neghess - neghess_fd).abs() <= 1e-3 * (1.0 + neghess.abs()),
8374 "{label}: neghessian_eta {neghess} vs FD {neghess_fd}"
8375 );
8376
8377 let nhd_fd = (sp.neghessian_eta[0] - sm.neghessian_eta[0]) / (2.0 * h);
8379 let nhd = s0.neghessian_eta_derivative[0];
8380 assert!(
8381 (nhd - nhd_fd).abs() <= 1e-2 * (1.0 + nhd.abs()),
8382 "{label}: neghessian_eta_derivative {nhd} vs FD {nhd_fd}"
8383 );
8384 }
8385
8386 #[test]
8387 fn poisson_gamma_nb_tweedie_arms_match_finite_differences_1615_1616() {
8388 let log = InverseLink::Standard(StandardLink::Log);
8389 let poisson = LikelihoodSpec {
8390 response: ResponseFamily::Poisson,
8391 link: log.clone(),
8392 };
8393 check_fd("poisson y=3", &poisson, 3.0, 0.4);
8394 check_fd("poisson y=0", &poisson, 0.0, -0.2);
8395
8396 let gamma = LikelihoodSpec {
8397 response: ResponseFamily::Gamma,
8398 link: log.clone(),
8399 };
8400 check_fd("gamma y=2.5", &gamma, 2.5, 0.3);
8401 check_fd("gamma y=0.7", &gamma, 0.7, -0.1);
8402
8403 let nb = LikelihoodSpec {
8404 response: ResponseFamily::NegativeBinomial {
8405 theta: 1.5,
8406 theta_fixed: true,
8407 },
8408 link: log.clone(),
8409 };
8410 check_fd("negbin y=4", &nb, 4.0, 0.5);
8411 check_fd("negbin y=0", &nb, 0.0, -0.3);
8412
8413 let tweedie = LikelihoodSpec {
8414 response: ResponseFamily::Tweedie { p: 1.5 },
8415 link: log.clone(),
8416 };
8417 check_fd("tweedie y=2", &tweedie, 2.0, 0.25);
8418 check_fd("tweedie y=0.5", &tweedie, 0.5, -0.15);
8419 }
8420
8421 #[test]
8422 fn binomial_natural_coordinate_towers_match_finite_differences() {
8423 for (label, family, eta) in [
8424 ("logit", LikelihoodSpec::binomial_logit(), 0.7),
8425 ("probit", LikelihoodSpec::binomial_probit(), -1.1),
8426 ("cloglog", LikelihoodSpec::binomial_cloglog(), 0.4),
8427 (
8428 "loglog",
8429 LikelihoodSpec::try_new(
8430 ResponseFamily::Binomial,
8431 InverseLink::Standard(StandardLink::LogLog),
8432 )
8433 .unwrap(),
8434 -0.35,
8435 ),
8436 (
8437 "cauchit",
8438 LikelihoodSpec::try_new(
8439 ResponseFamily::Binomial,
8440 InverseLink::Standard(StandardLink::Cauchit),
8441 )
8442 .unwrap(),
8443 1.25,
8444 ),
8445 ] {
8446 check_fd(label, &family, 0.37, eta);
8447 }
8448 }
8449
8450 #[test]
8451 fn logit_observation_geometry_carries_the_prior_weight_everywhere() {
8452 let eta = 1.75;
8453 let y = 0.3;
8454 let weight = 7.25;
8455 let state = one_obs_weight(&LikelihoodSpec::binomial_logit(), y, weight, eta);
8456 let jet = logit_inverse_link_jet5(eta);
8457 for (got, expected) in [
8458 (state.fisherweight[0], weight * jet.d1),
8459 (state.neghessian_eta[0], weight * jet.d1),
8460 (state.neghessian_eta_derivative[0], weight * jet.d2),
8461 (state.score[0], weight * (y - jet.mu)),
8462 ] {
8463 assert!((got - expected).abs() <= 4.0 * f64::EPSILON * (1.0 + expected.abs()));
8464 }
8465 }
8466
8467 #[test]
8468 fn tiny_positive_and_zero_weights_are_not_projected() {
8469 let tiny = 1e-200;
8470 let logit = one_obs_weight(&LikelihoodSpec::binomial_logit(), 0.4, tiny, 0.0);
8471 assert!((logit.fisherweight[0] / tiny - 0.25).abs() <= 2.0 * f64::EPSILON);
8472 assert!(logit.fisherweight[0] < 1e-190);
8473
8474 let zero = one_obs_weight(&LikelihoodSpec::gaussian_identity(), 3.0, 0.0, -2.0);
8475 assert_eq!(zero.score[0], 0.0);
8476 assert_eq!(zero.fisherweight[0], 0.0);
8477 assert_eq!(zero.neghessian_eta[0], 0.0);
8478 assert_eq!(zero.neghessian_eta_derivative[0], 0.0);
8479 assert_eq!(zero.log_likelihood, 0.0);
8480 assert_eq!(exact_standard_working_response(&zero).unwrap()[0], -2.0);
8481 }
8482
8483 #[test]
8484 fn log_link_tails_balance_tiny_weights_before_certification() {
8485 let poisson = one_obs_weight(&LikelihoodSpec::poisson_log(), 0.0, 1e-300, 700.0);
8486 assert!(poisson.fisherweight[0].is_finite() && poisson.fisherweight[0] > 1.0);
8487 assert!(poisson.score[0].is_finite());
8488 assert!(poisson.log_likelihood.is_finite());
8489
8490 let gamma = one_obs_weight(&LikelihoodSpec::gamma_log(), 1.0, 1e-300, -700.0);
8491 assert!(gamma.neghessian_eta[0].is_finite() && gamma.neghessian_eta[0] > 1.0);
8492 assert!(gamma.score[0].is_finite());
8493 assert!(gamma.log_likelihood.is_finite());
8494 }
8495
8496 #[test]
8497 fn invalid_weights_and_nonfinite_inputs_are_refused_in_row_order() {
8498 let family = LikelihoodSpec::gaussian_identity();
8499 let y = array![1.0, 2.0];
8500 let eta = array![0.0, 0.0];
8501 for weights in [array![-1.0, 1.0], array![f64::NAN, 1.0]] {
8502 let err = evaluate_standard_familyobservations(
8503 family.clone(),
8504 None,
8505 None,
8506 None,
8507 &y,
8508 &weights,
8509 &eta,
8510 )
8511 .expect_err("invalid prior weight must be refused");
8512 assert!(err.to_string().contains("row 0"), "{err}");
8513 }
8514
8515 let err = evaluate_standard_familyobservations(
8516 family,
8517 None,
8518 None,
8519 None,
8520 &array![f64::NAN],
8521 &array![0.0],
8522 &array![0.0],
8523 )
8524 .expect_err("a non-finite response may not hide behind zero weight");
8525 assert!(err.to_string().contains("row 0"), "{err}");
8526 }
8527
8528 #[test]
8529 fn unrepresentable_cloglog_curvature_is_refused_without_a_floor() {
8530 let err = evaluate_standard_familyobservations(
8531 LikelihoodSpec::binomial_cloglog(),
8532 None,
8533 None,
8534 None,
8535 &array![1.0],
8536 &array![1.0],
8537 &array![18.0],
8538 )
8539 .expect_err("mathematically sub-f64 Fisher information must be refused");
8540 assert!(err.to_string().contains("Fisher weight"), "{err}");
8541 }
8542
8543 #[test]
8544 fn bounded_covariance_requires_a_certified_strict_spd_precision() {
8545 let covariance = certified_bounded_posterior_covariance(
8546 &array![[4.0, 1.0], [1.0, 3.0]],
8547 "bounded covariance regression",
8548 )
8549 .expect("strict SPD precision");
8550 assert!((covariance[[0, 0]] - 3.0 / 11.0).abs() < 1e-14);
8551 assert!((covariance[[0, 1]] + 1.0 / 11.0).abs() < 1e-14);
8552 assert!((covariance[[1, 1]] - 4.0 / 11.0).abs() < 1e-14);
8553
8554 for invalid in [
8555 array![[1.0, 1.0], [1.0, 1.0]],
8556 array![[1.0, 2.0], [2.0, 1.0]],
8557 ] {
8558 assert!(
8559 certified_bounded_posterior_covariance(
8560 &invalid,
8561 "invalid bounded covariance regression"
8562 )
8563 .is_err(),
8564 "singular/indefinite precision must not become a pseudo-covariance"
8565 );
8566 }
8567 }
8568}