1fn try_build_spatial_term_log_kappa_derivative(
2 data: ArrayView2<'_, f64>,
3 resolvedspec: &TermCollectionSpec,
4 design: &TermCollectionDesign,
5 term_idx: usize,
6) -> Result<
7 Option<(
8 Range<usize>,
9 usize,
10 Array2<f64>,
11 Array2<f64>,
12 Array2<f64>,
13 Array2<f64>,
14 Vec<Array2<f64>>,
15 Vec<Array2<f64>>,
16 Option<std::sync::Arc<gam_terms::basis::ImplicitDesignPsiDerivative>>,
17 )>,
18 EstimationError,
19> {
20 let Some(smooth_term) = design.smooth.terms.get(term_idx) else {
21 return Ok(None);
22 };
23 let Some(termspec) = resolvedspec.smooth_terms.get(term_idx) else {
24 return Ok(None);
25 };
26
27 let derivative_bundle = match &termspec.basis {
28 SmoothBasisSpec::ThinPlate {
29 feature_cols,
30 spec,
31 input_scale,
32 } => {
33 let mut x = select_columns(data, feature_cols).map_err(EstimationError::from)?;
34 let mut spec_local = spec.clone();
35 if let Some(scale) = input_scale {
36 scale.standardize(&mut x);
37 spec_local.length_scale = scale
38 .to_standardized_units(gam_terms::OriginalUnits::new(spec.length_scale))
39 .standardized_value();
40 }
41 build_thin_plate_basis_log_kappa_derivatives(x.view(), &spec_local)
42 .map_err(EstimationError::from)?
43 }
44 SmoothBasisSpec::Sphere { .. } => return Ok(None),
45 SmoothBasisSpec::ConstantCurvature { feature_cols, spec } => {
54 let x = select_columns(data, feature_cols).map_err(EstimationError::from)?;
55 build_constant_curvature_basis_kappa_derivatives(x.view(), spec)
56 .map_err(EstimationError::from)?
57 }
58 SmoothBasisSpec::MeasureJet { .. } => return Ok(None),
64 SmoothBasisSpec::Matern {
65 feature_cols,
66 spec,
67 input_scale,
68 } => {
69 let mut x = select_columns(data, feature_cols).map_err(EstimationError::from)?;
70 let mut spec_local = spec.clone();
71 if let Some(scale) = input_scale {
72 scale.standardize(&mut x);
73 let length_scale = spec.length_scale.resolved().ok_or_else(|| {
74 EstimationError::InvalidInput(
75 "Matérn Auto length_scale reached derivative construction unresolved"
76 .to_string(),
77 )
78 })?;
79 spec_local.length_scale.set_resolved(
80 scale
81 .to_standardized_units(gam_terms::OriginalUnits::new(length_scale))
82 .standardized_value(),
83 );
84 }
85 spec_local.double_penalty = false;
100 build_matern_basis_log_kappa_derivatives(x.view(), &spec_local)
101 .map_err(EstimationError::from)?
102 }
103 SmoothBasisSpec::Duchon {
104 feature_cols,
105 spec,
106 input_scale,
107 } => {
108 let mut x = select_columns(data, feature_cols).map_err(EstimationError::from)?;
109 let mut spec_local = spec.clone();
110 if let Some(scale) = input_scale {
111 scale.standardize(&mut x);
112 spec_local.length_scale = spec.length_scale.map(|length| {
113 scale
114 .to_standardized_units(gam_terms::OriginalUnits::new(length))
115 .standardized_value()
116 });
117 }
118 let BasisMetadata::Duchon {
119 centers,
120 identifiability_transform,
121 operator_collocation_points,
122 radial_reparam,
123 ..
124 } = &smooth_term.metadata
125 else {
126 return Ok(None);
127 };
128 if spec_local.radial_reparam.is_none() {
131 spec_local.radial_reparam = radial_reparam.clone();
132 }
133 gam_terms::basis::build_duchon_basis_log_kappa_derivativeswith_collocationwithworkspace(
134 x.view(),
135 &spec_local,
136 centers.view(),
137 identifiability_transform.as_ref(),
138 operator_collocation_points
139 .as_ref()
140 .map(|points| points.view()),
141 &mut BasisWorkspace::default(),
142 )
143 .map_err(EstimationError::from)?
144 }
145 SmoothBasisSpec::BSpline1D { .. }
146 | SmoothBasisSpec::TensorBSpline { .. }
147 | SmoothBasisSpec::ByVariable { .. }
148 | SmoothBasisSpec::FactorSumToZero { .. }
149 | SmoothBasisSpec::BySmooth { .. }
150 | SmoothBasisSpec::FactorSmooth { .. }
151 | SmoothBasisSpec::Pca { .. } => {
152 return Ok(None);
153 }
154 };
155 let mut implicit_operator = derivative_bundle.implicit_operator;
156 let BasisPsiDerivativeResult {
157 design_derivative: mut local_x_psi,
158 penalties_derivative: mut local_s_psi,
159 implicit_operator: local_implicit_first_unused,
160 } = derivative_bundle.first;
161 let BasisPsiSecondDerivativeResult {
162 designsecond_derivative: mut local_x_psi_psi,
163 penaltiessecond_derivative: mut local_s_psi_psi,
164 implicit_operator: local_implicit_second_unused,
165 } = derivative_bundle.second;
166 assert!(local_implicit_first_unused.is_none());
167 assert!(local_implicit_second_unused.is_none());
168
169 if let Some(rotation) = smooth_term.joint_null_rotation.as_ref() {
170 let q = &rotation.rotation;
171 if let Some(op) = implicit_operator.take() {
172 implicit_operator = Some(op.append_full_transform(q).map_err(EstimationError::from)?);
173 } else {
174 if local_x_psi.ncols() != q.nrows() || local_x_psi_psi.ncols() != q.nrows() {
175 return Ok(None);
176 }
177 local_x_psi = fast_ab(&local_x_psi, q);
178 local_x_psi_psi = fast_ab(&local_x_psi_psi, q);
179 }
180 let rotate_penalty = |s_local: Array2<f64>| -> Option<Array2<f64>> {
181 if s_local.nrows() != q.nrows() || s_local.ncols() != q.nrows() {
182 return None;
183 }
184 let qt_s = gam_linalg::faer_ndarray::fast_atb(q, &s_local);
185 Some(gam_linalg::faer_ndarray::fast_ab(&qt_s, q))
186 };
187 let Some(rotated_s_psi) = local_s_psi
188 .into_iter()
189 .map(|s| rotate_penalty(s))
190 .collect::<Option<Vec<_>>>()
191 else {
192 return Ok(None);
193 };
194 local_s_psi = rotated_s_psi;
195 let Some(rotated_s_psi_psi) = local_s_psi_psi
196 .into_iter()
197 .map(|s| rotate_penalty(s))
198 .collect::<Option<Vec<_>>>()
199 else {
200 return Ok(None);
201 };
202 local_s_psi_psi = rotated_s_psi_psi;
203 }
204 let implicit_operator = implicit_operator.map(std::sync::Arc::new);
205
206 if let Some(ref op) = implicit_operator {
207 if op.p_out() != smooth_term.coeff_range.len() {
208 return Ok(None);
209 }
210 } else {
211 if local_x_psi.ncols() != smooth_term.coeff_range.len() {
212 return Ok(None);
213 }
214 if local_x_psi_psi.ncols() != smooth_term.coeff_range.len() {
215 return Ok(None);
216 }
217 }
218 if local_s_psi.is_empty() || local_s_psi.len() != local_s_psi_psi.len() {
219 return Ok(None);
220 }
221 if local_s_psi.iter().any(|s| {
222 s.nrows() != smooth_term.coeff_range.len() || s.ncols() != smooth_term.coeff_range.len()
223 }) {
224 return Ok(None);
225 }
226 if local_s_psi_psi.iter().any(|s| {
227 s.nrows() != smooth_term.coeff_range.len() || s.ncols() != smooth_term.coeff_range.len()
228 }) {
229 return Ok(None);
230 }
231
232 let p_total = design.design.ncols();
233 let smooth_start = p_total.saturating_sub(design.smooth.total_smooth_cols());
234 let global_range = (smooth_start + smooth_term.coeff_range.start)
235 ..(smooth_start + smooth_term.coeff_range.end);
236
237 Ok(Some((
238 global_range,
239 p_total,
240 local_x_psi,
241 local_s_psi.iter().fold(
242 Array2::<f64>::zeros((smooth_term.coeff_range.len(), smooth_term.coeff_range.len())),
243 |acc, m| acc + m,
244 ),
245 local_x_psi_psi,
246 local_s_psi_psi.iter().fold(
247 Array2::<f64>::zeros((smooth_term.coeff_range.len(), smooth_term.coeff_range.len())),
248 |acc, m| acc + m,
249 ),
250 local_s_psi,
251 local_s_psi_psi,
252 implicit_operator,
253 )))
254}
255
256fn try_build_spatial_log_kappa_hyper_dirs(
257 data: ArrayView2<'_, f64>,
258 resolvedspec: &TermCollectionSpec,
259 design: &TermCollectionDesign,
260 spatial_terms: &[usize],
261) -> Result<Option<Vec<DirectionalHyperParam>>, EstimationError> {
262 let Some(info_list) =
269 try_build_spatial_log_kappa_derivativeinfo_list(data, resolvedspec, design, spatial_terms)?
270 else {
271 return Ok(None);
272 };
273 Ok(Some(spatial_log_kappa_hyper_dirs_frominfo_list(info_list)?))
274}
275
276pub(crate) fn try_build_latent_coord_hyper_dirs(
277 latent: std::sync::Arc<gam_terms::latent::LatentCoordValues>,
278 resolvedspec: &TermCollectionSpec,
279 design: &TermCollectionDesign,
280 latent_terms: &[gam_problem::types::SmoothTermIdx],
281 analytic_rho_count: usize,
282) -> Result<Option<Vec<DirectionalHyperParam>>, EstimationError> {
283 if latent_terms.is_empty() || latent.is_empty() {
284 return Ok(None);
285 }
286 if latent_terms.len() != 1 {
287 crate::bail_invalid_estim!(
288 "LatentCoord standard-fit hyper_dirs currently require exactly one latent smooth term"
289 .to_string(),
290 );
291 }
292 let term_idx = latent_terms[0];
293 let smooth_term = design.smooth.terms.get(term_idx.get()).ok_or_else(|| {
294 EstimationError::InvalidInput(format!(
295 "LatentCoord term index {term_idx} out of bounds for realized smooth design"
296 ))
297 })?;
298 let termspec = resolvedspec
299 .smooth_terms
300 .get(term_idx.get())
301 .ok_or_else(|| {
302 EstimationError::InvalidInput(format!(
303 "LatentCoord term index {term_idx} out of bounds for resolved smooth spec"
304 ))
305 })?;
306 let p_total = design.design.ncols();
307 let smooth_start = p_total.saturating_sub(design.smooth.total_smooth_cols());
308 let global_range = (smooth_start + smooth_term.coeff_range.start)
309 ..(smooth_start + smooth_term.coeff_range.end);
310
311 let operator = match (&termspec.basis, &smooth_term.metadata) {
316 (
317 SmoothBasisSpec::Matern { .. },
318 BasisMetadata::Matern {
319 centers,
320 length_scale,
321 nu,
322 include_intercept,
323 identifiability_transform,
324 input_scale,
325 ..
326 },
327 ) => gam_terms::basis::LatentCoordDesignDerivative::new_matern(
328 latent.clone(),
329 std::sync::Arc::new(centers.clone()),
330 *input_scale,
334 *length_scale,
335 *nu,
336 *include_intercept,
337 identifiability_transform.clone(),
338 )
339 .map_err(EstimationError::from)?,
340 (
341 SmoothBasisSpec::Duchon { .. },
342 BasisMetadata::Duchon {
343 centers,
344 length_scale,
345 power,
346 nullspace_order,
347 identifiability_transform,
348 input_scale,
349 ..
350 },
351 ) => gam_terms::basis::LatentCoordDesignDerivative::new_duchon(
352 latent.clone(),
353 std::sync::Arc::new(centers.clone()),
354 *input_scale,
356 *length_scale,
357 *power,
358 *nullspace_order,
359 identifiability_transform.clone(),
360 )
361 .map_err(EstimationError::from)?,
362 (
363 SmoothBasisSpec::Sphere { .. },
364 BasisMetadata::Sphere {
365 centers,
366 penalty_order,
367 method,
368 constraint_transform,
369 ..
370 },
371 ) if matches!(*method, gam_terms::basis::SphereMethod::Wahba) => {
372 gam_terms::basis::LatentCoordDesignDerivative::new_sphere(
373 latent.clone(),
374 std::sync::Arc::new(centers.clone()),
375 *penalty_order,
376 constraint_transform.clone(),
377 )
378 .map_err(EstimationError::from)?
379 }
380 (
381 SmoothBasisSpec::BSpline1D { spec, .. },
382 BasisMetadata::BSpline1D {
383 knots,
384 identifiability_transform,
385 periodic,
386 degree: meta_degree,
387 ..
388 },
389 ) => {
390 let effective_degree = meta_degree.unwrap_or(spec.degree);
394 if let Some((domain_start, period, num_basis)) = periodic {
395 gam_terms::basis::LatentCoordDesignDerivative::new_periodic_bspline(
396 latent.clone(),
397 (*domain_start, *domain_start + *period),
398 effective_degree,
399 *num_basis,
400 identifiability_transform.clone(),
401 )
402 .map_err(EstimationError::from)?
403 } else {
404 gam_terms::basis::LatentCoordDesignDerivative::new_tensor_bspline(
405 latent.clone(),
406 vec![knots.clone()],
407 vec![effective_degree],
408 identifiability_transform.clone(),
409 )
410 .map_err(EstimationError::from)?
411 }
412 }
413 (
414 SmoothBasisSpec::TensorBSpline { .. },
415 BasisMetadata::TensorBSpline {
416 knots,
417 degrees,
418 identifiability_transform,
419 ..
420 },
421 ) => gam_terms::basis::LatentCoordDesignDerivative::new_tensor_bspline(
422 latent.clone(),
423 knots.clone(),
424 degrees.clone(),
425 identifiability_transform.clone(),
426 )
427 .map_err(EstimationError::from)?,
428 (SmoothBasisSpec::Pca { .. }, BasisMetadata::Pca { basis_matrix, .. }) => {
429 gam_terms::basis::LatentCoordDesignDerivative::new_pca(
430 latent.clone(),
431 std::sync::Arc::new(basis_matrix.clone()),
432 )
433 .map_err(EstimationError::from)?
434 }
435 _ => return Ok(None),
436 };
437 if operator.p_out() != global_range.len() {
438 crate::bail_invalid_estim!(
439 "LatentCoord derivative width mismatch for term '{}': operator p={}, coeff range={}",
440 smooth_term.name,
441 operator.p_out(),
442 global_range.len()
443 );
444 }
445 let operator = std::sync::Arc::new(operator);
446 let mut hyper_dirs = Vec::with_capacity(operator.n_axes());
447 for flat_axis in 0..operator.n_axes() {
448 let dir = DirectionalHyperParam::new_compact(
449 gam_solve::estimate::reml::HyperDesignDerivative::from_latent_coord(
450 operator.clone(),
451 flat_axis,
452 global_range.clone(),
453 p_total,
454 ),
455 Vec::new(),
456 None,
457 None,
458 )?
459 .not_penalty_like();
460 hyper_dirs.push(dir);
461 }
462 let direct_dim = latent_coord_direct_hyper_count(latent.id_mode(), latent.latent_dim());
463 if analytic_rho_count + direct_dim > 0 {
464 let zero_x = gam_solve::estimate::reml::HyperDesignDerivative::from(Array2::<f64>::zeros(
465 (design.design.nrows(), p_total),
466 ));
467 for _ in 0..analytic_rho_count {
468 hyper_dirs.push(
469 DirectionalHyperParam::new_compact(zero_x.clone(), Vec::new(), None, None)?
470 .not_penalty_like(),
471 );
472 }
473 for _ in 0..direct_dim {
474 hyper_dirs.push(
475 DirectionalHyperParam::new_compact(zero_x.clone(), Vec::new(), None, None)?
476 .not_penalty_like(),
477 );
478 }
479 }
480 Ok(Some(hyper_dirs))
481}
482
483fn latent_coord_direct_hyper_count(
484 id_mode: &gam_terms::latent::LatentIdMode,
485 latent_dim: usize,
486) -> usize {
487 use gam_terms::latent::{AuxPriorStrength, LatentIdMode};
488 match id_mode {
489 LatentIdMode::AuxPrior { strength, .. } => match strength {
490 AuxPriorStrength::Auto => 1,
491 AuxPriorStrength::Fixed(_) => 0,
492 },
493 LatentIdMode::AuxPriorDimSelection { strength, .. } => {
494 latent_dim
495 + match strength {
496 AuxPriorStrength::Auto => 1,
497 AuxPriorStrength::Fixed(_) => 0,
498 }
499 }
500 LatentIdMode::DimSelection { .. } => latent_dim,
501 LatentIdMode::IsometryToReference { strength, .. } => match strength {
504 AuxPriorStrength::Auto => 1,
505 AuxPriorStrength::Fixed(_) => 0,
506 },
507 LatentIdMode::AuxOutcome { head, .. } => head.n_coeffs(latent_dim) + latent_dim,
510 LatentIdMode::None => 0,
511 }
512}
513
514fn latent_coord_initial_direct_hypers(
515 id_mode: &gam_terms::latent::LatentIdMode,
516 latent_dim: usize,
517) -> Result<Array1<f64>, EstimationError> {
518 use gam_terms::latent::{AuxPriorStrength, LatentIdMode};
519 let mut values = Vec::with_capacity(latent_coord_direct_hyper_count(id_mode, latent_dim));
520 match id_mode {
521 LatentIdMode::AuxPrior { strength, .. } => {
522 if matches!(strength, AuxPriorStrength::Auto) {
523 values.push(0.0);
524 }
525 }
526 LatentIdMode::AuxPriorDimSelection {
527 strength,
528 init_log_precision,
529 ..
530 } => {
531 if matches!(strength, AuxPriorStrength::Auto) {
532 values.push(0.0);
533 }
534 append_latent_ard_seed(&mut values, init_log_precision.as_ref(), latent_dim)?;
535 }
536 LatentIdMode::DimSelection { init_log_precision } => {
537 append_latent_ard_seed(&mut values, init_log_precision.as_ref(), latent_dim)?;
538 }
539 LatentIdMode::IsometryToReference { strength, .. } => {
540 if matches!(strength, AuxPriorStrength::Auto) {
541 values.push(0.0);
542 }
543 }
544 LatentIdMode::AuxOutcome {
545 head,
546 init_log_precision,
547 } => {
548 values.extend(std::iter::repeat_n(0.0, head.n_coeffs(latent_dim)));
552 append_latent_ard_seed(&mut values, init_log_precision.as_ref(), latent_dim)?;
553 }
554 LatentIdMode::None => {}
555 }
556 Ok(Array1::from_vec(values))
557}
558
559fn append_latent_ard_seed(
560 values: &mut Vec<f64>,
561 init: Option<&Array1<f64>>,
562 latent_dim: usize,
563) -> Result<(), EstimationError> {
564 if let Some(init) = init {
565 if init.len() != latent_dim {
566 crate::bail_invalid_estim!(
567 "latent dim_selection init_log_precision length mismatch: got {}, expected {}",
568 init.len(),
569 latent_dim
570 );
571 }
572 values.extend(init.iter().copied());
573 } else {
574 values.extend(std::iter::repeat_n(0.0, latent_dim));
575 }
576 Ok(())
577}
578
579struct LatentIdObjectiveContribution {
580 cost: f64,
581 gradient: Array1<f64>,
582}
583
584fn latent_id_objective_contribution(
585 theta: &Array1<f64>,
586 rho_dim: usize,
587 analytic_rho_count: usize,
588 latent: &gam_terms::latent::LatentCoordValues,
589) -> Result<LatentIdObjectiveContribution, EstimationError> {
590 use gam_terms::latent::{AuxPriorStrength, LatentIdMode, aux_prior_targets};
591 let n_obs = latent.n_obs();
592 let latent_dim = latent.latent_dim();
593 let flat_len = latent.len();
594 let mut gradient = Array1::<f64>::zeros(theta.len());
595 let t_start = rho_dim;
596 let direct_start = t_start + flat_len + analytic_rho_count;
597 if theta.len() < direct_start {
598 crate::bail_invalid_estim!(
599 "latent-coordinate theta too short for id objective: got {}, need at least {}",
600 theta.len(),
601 direct_start
602 );
603 }
604 let t = latent.as_matrix();
605 let mut cost = 0.0;
606 let mut cursor = direct_start;
607
608 match latent.id_mode() {
609 LatentIdMode::AuxPrior {
610 u,
611 family,
612 strength,
613 }
614 | LatentIdMode::AuxPriorDimSelection {
615 u,
616 family,
617 strength,
618 ..
619 } => {
620 let (log_mu, mu) = match strength {
621 AuxPriorStrength::Fixed(mu) => (
622 gam_problem::checked_log_strength(*mu).map_err(|error| {
623 EstimationError::InvalidInput(format!(
624 "fixed latent auxiliary-prior precision is outside the canonical physical-strength domain: {error}"
625 ))
626 })?,
627 *mu,
628 ),
629 AuxPriorStrength::Auto => {
630 let log_mu = *theta.get(cursor).ok_or_else(|| {
631 EstimationError::InvalidInput(format!(
632 "latent auxiliary-prior precision coordinate {cursor} is missing from theta length {}",
633 theta.len(),
634 ))
635 })?;
636 cursor += 1;
637 let mu = gam_problem::checked_exp_log_strength(log_mu).map_err(|error| {
638 EstimationError::InvalidInput(format!(
639 "latent auxiliary-prior log precision is outside the canonical log-strength domain: {error}"
640 ))
641 })?;
642 (log_mu, mu)
643 }
644 };
645 let targets = aux_prior_targets(t.view(), u.view(), *family)
646 .map_err(EstimationError::InvalidInput)?;
647 let residual = &t - &targets;
648 let q = residual.iter().map(|v| v * v).sum::<f64>();
649 let k = (n_obs * latent_dim) as f64;
656 cost += 0.5 * mu * q - 0.5 * k * log_mu;
657
658 let projected_residual = aux_prior_targets(residual.view(), u.view(), *family)
659 .map_err(EstimationError::InvalidInput)?;
660 let grad_base = residual - projected_residual;
661 for n in 0..n_obs {
662 for axis in 0..latent_dim {
663 gradient[t_start + n * latent_dim + axis] += mu * grad_base[[n, axis]];
664 }
665 }
666 if matches!(strength, AuxPriorStrength::Auto) {
667 gradient[direct_start] += 0.5 * mu * q - 0.5 * k;
668 }
669 }
670 LatentIdMode::IsometryToReference {
671 reference,
672 strength,
673 } => {
674 if reference.dim() != (n_obs, latent_dim) {
681 crate::bail_invalid_estim!(
682 "IsometryToReference reference shape {:?} must equal (n_obs, latent_dim) = ({}, {})",
683 reference.dim(),
684 n_obs,
685 latent_dim
686 );
687 }
688 let mu_slot = cursor;
689 let (log_mu, mu) = match strength {
690 AuxPriorStrength::Fixed(mu) => (
691 gam_problem::checked_log_strength(*mu).map_err(|error| {
692 EstimationError::InvalidInput(format!(
693 "fixed latent isometry precision is outside the canonical physical-strength domain: {error}"
694 ))
695 })?,
696 *mu,
697 ),
698 AuxPriorStrength::Auto => {
699 let log_mu = *theta.get(cursor).ok_or_else(|| {
700 EstimationError::InvalidInput(format!(
701 "latent isometry precision coordinate {cursor} is missing from theta length {}",
702 theta.len(),
703 ))
704 })?;
705 cursor += 1;
706 let mu = gam_problem::checked_exp_log_strength(log_mu).map_err(|error| {
707 EstimationError::InvalidInput(format!(
708 "latent isometry log precision is outside the canonical log-strength domain: {error}"
709 ))
710 })?;
711 (log_mu, mu)
712 }
713 };
714 let residual = &t - reference;
715 let q = residual.iter().map(|v| v * v).sum::<f64>();
716 let k = (n_obs * latent_dim) as f64;
720 cost += 0.5 * mu * q - 0.5 * k * log_mu;
721 for n in 0..n_obs {
722 for axis in 0..latent_dim {
723 gradient[t_start + n * latent_dim + axis] += mu * residual[[n, axis]];
724 }
725 }
726 if matches!(strength, AuxPriorStrength::Auto) {
727 gradient[mu_slot] += 0.5 * mu * q - 0.5 * k;
728 }
729 }
730 LatentIdMode::AuxOutcome { head, .. } => {
731 let n_coeffs = head.n_coeffs(latent_dim);
739 if cursor + n_coeffs > theta.len() {
740 crate::bail_invalid_estim!(
741 "latent auxiliary-outcome coefficient block overruns theta: start={cursor}, width={n_coeffs}, theta_len={}",
742 theta.len(),
743 );
744 }
745 let coeffs = theta
746 .slice(ndarray::s![cursor..cursor + n_coeffs])
747 .to_owned();
748 let (head_nll, grad_coeffs, grad_t) = head
749 .neg_loglik_and_grad(t.view(), coeffs.view())
750 .map_err(EstimationError::InvalidInput)?;
751 cost += head_nll;
752 for (offset, &g) in grad_coeffs.iter().enumerate() {
753 gradient[cursor + offset] += g;
754 }
755 for n in 0..n_obs {
756 for axis in 0..latent_dim {
757 gradient[t_start + n * latent_dim + axis] += grad_t[[n, axis]];
758 }
759 }
760 cursor += n_coeffs;
761 }
762 LatentIdMode::DimSelection { .. } | LatentIdMode::None => {}
763 }
764
765 match latent.id_mode() {
766 LatentIdMode::AuxPriorDimSelection { .. }
767 | LatentIdMode::DimSelection { .. }
768 | LatentIdMode::AuxOutcome { .. } => {
769 if cursor + latent_dim > theta.len() {
770 crate::bail_invalid_estim!(
771 "latent dimension-selection precision block overruns theta: start={cursor}, width={latent_dim}, theta_len={}",
772 theta.len(),
773 );
774 }
775 let alphas = gam_problem::checked_exp_log_strengths(
776 theta.slice(s![cursor..cursor + latent_dim]).iter().copied(),
777 )
778 .map_err(|error| {
779 EstimationError::InvalidInput(format!(
780 "latent dimension-selection log precision is outside the canonical log-strength domain: {error}"
781 ))
782 })?;
783 for axis in 0..latent_dim {
784 let log_alpha = theta[cursor + axis];
785 let alpha = alphas[axis];
786 let mut q_axis = 0.0;
787 for n in 0..n_obs {
788 let flat_idx = n * latent_dim + axis;
789 let value = latent.as_flat()[flat_idx];
790 q_axis += value * value;
791 gradient[t_start + flat_idx] += alpha * value;
792 }
793 cost += 0.5 * alpha * q_axis - 0.5 * n_obs as f64 * log_alpha;
794 gradient[cursor + axis] += 0.5 * alpha * q_axis - 0.5 * n_obs as f64;
795 }
796 cursor += latent_dim;
797 }
798 LatentIdMode::AuxPrior { .. }
799 | LatentIdMode::IsometryToReference { .. }
800 | LatentIdMode::None => {}
801 }
802
803 if cursor != theta.len() {
804 crate::bail_invalid_estim!(
805 "latent-coordinate direct hyperparameter length mismatch: consumed {}, theta len {}",
806 cursor,
807 theta.len()
808 );
809 }
810 Ok(LatentIdObjectiveContribution { cost, gradient })
811}
812
813fn add_latent_id_objective_to_eval(
814 theta: &Array1<f64>,
815 rho_dim: usize,
816 analytic_rho_count: usize,
817 latent: &gam_terms::latent::LatentCoordValues,
818 eval: &mut (f64, Array1<f64>, gam_problem::HessianValue),
819) -> Result<(), EstimationError> {
820 let contribution =
821 latent_id_objective_contribution(theta, rho_dim, analytic_rho_count, latent)?;
822 eval.0 += contribution.cost;
823 if eval.1.len() != contribution.gradient.len() {
824 crate::bail_invalid_estim!(
825 "latent-coordinate REML gradient length mismatch: base={}, id={}",
826 eval.1.len(),
827 contribution.gradient.len()
828 );
829 }
830 eval.1 += &contribution.gradient;
831 if eval.2.is_analytic() {
832 eval.2 = gam_problem::HessianValue::Unavailable;
833 }
834 Ok(())
835}
836
837fn analytic_penalty_objective_contribution(
838 theta: &Array1<f64>,
839 rho_dim: usize,
840 latent: &gam_terms::latent::LatentCoordValues,
841 registry: &gam_terms::AnalyticPenaltyRegistry,
842) -> Result<LatentIdObjectiveContribution, EstimationError> {
843 let flat_len = latent.len();
844 let t_start = rho_dim;
845 let t_end = t_start + flat_len;
846 let rho_start = t_end;
847 let rho_end = rho_start + registry.total_rho_count();
848 if theta.len() < rho_end {
849 crate::bail_invalid_estim!(
850 "latent-coordinate theta too short for analytic penalties: got {}, need at least {}",
851 theta.len(),
852 rho_end
853 );
854 }
855 let target_t = theta.slice(s![t_start..t_end]);
856 let rho = theta.slice(s![rho_start..rho_end]);
857 registry
858 .validate_rho(rho)
859 .map_err(EstimationError::InvalidInput)?;
860 let mut cost = 0.0_f64;
861 let mut gradient = Array1::<f64>::zeros(theta.len());
862 for (penalty, (rho_slice, tier, name)) in registry.penalties.iter().zip(registry.rho_layout()) {
863 let rho_local = rho.slice(s![rho_slice.clone()]);
864 match tier {
865 gam_terms::PenaltyTier::Psi => {
866 cost += penalty.value(target_t.view(), rho_local);
867 let grad = penalty.grad_target(target_t.view(), rho_local);
868 if grad.len() != flat_len {
869 crate::bail_invalid_estim!(
870 "analytic penalty {name:?} gradient length mismatch: got {}, expected {}",
871 grad.len(),
872 flat_len
873 );
874 }
875 for i in 0..flat_len {
876 gradient[t_start + i] += grad[i];
877 }
878 let grad_rho_local = penalty.grad_rho(target_t.view(), rho_local);
879 if grad_rho_local.len() != rho_slice.len() {
880 crate::bail_invalid_estim!(
881 "analytic penalty {name:?} rho-gradient length mismatch: got {}, expected {}",
882 grad_rho_local.len(),
883 rho_slice.len()
884 );
885 }
886 for local_idx in 0..grad_rho_local.len() {
887 gradient[rho_start + rho_slice.start + local_idx] += grad_rho_local[local_idx];
888 }
889 }
890 gam_terms::PenaltyTier::Beta => {}
891 gam_terms::PenaltyTier::Rho => {}
892 }
893 }
894 Ok(LatentIdObjectiveContribution { cost, gradient })
895}
896
897fn add_analytic_penalty_hessian_to_eval(
898 theta: &Array1<f64>,
899 rho_dim: usize,
900 latent: &gam_terms::latent::LatentCoordValues,
901 registry: &gam_terms::AnalyticPenaltyRegistry,
902 eval: &mut (f64, Array1<f64>, gam_problem::HessianValue),
903) -> Result<(), EstimationError> {
904 let flat_len = latent.len();
905 let t_start = rho_dim;
906 let t_end = t_start + flat_len;
907 let rho_start = t_end;
908 let rho_end = rho_start + registry.total_rho_count();
909 if theta.len() < rho_end {
910 crate::bail_invalid_estim!(
911 "latent-coordinate theta too short for analytic penalty Hessian: got {}, need at least {}",
912 theta.len(),
913 rho_end
914 );
915 }
916 let gam_problem::HessianValue::Dense(hessian) = &mut eval.2 else {
917 if eval.2.is_analytic() {
918 eval.2 = gam_problem::HessianValue::Unavailable;
919 }
920 return Ok(());
921 };
922 if hessian.dim() != (theta.len(), theta.len()) {
923 crate::bail_invalid_estim!(
924 "analytic penalty Hessian target shape mismatch: got {}x{}, expected {}x{}",
925 hessian.nrows(),
926 hessian.ncols(),
927 theta.len(),
928 theta.len()
929 );
930 }
931 let target_t = theta.slice(s![t_start..t_end]);
932 let rho = theta.slice(s![rho_start..rho_end]);
933 registry
934 .validate_rho(rho)
935 .map_err(EstimationError::InvalidInput)?;
936 for (penalty, (rho_slice, tier, _name)) in registry.penalties.iter().zip(registry.rho_layout())
937 {
938 let rho_local = rho.slice(s![rho_slice]);
939 if !matches!(tier, gam_terms::PenaltyTier::Psi) {
940 continue;
941 }
942 if let Some(diag) = penalty.hessian_diag(target_t.view(), rho_local) {
943 if diag.len() != flat_len {
944 crate::bail_invalid_estim!(
945 "analytic penalty Hessian diagonal length mismatch: got {}, expected {}",
946 diag.len(),
947 flat_len
948 );
949 }
950 for i in 0..flat_len {
951 hessian[[t_start + i, t_start + i]] += diag[i];
952 }
953 continue;
954 }
955 let mut probe = Array1::<f64>::zeros(flat_len);
956 for col in 0..flat_len {
957 probe[col] = 1.0;
958 let hv = penalty.hvp(target_t.view(), rho_local, probe.view());
959 if hv.len() != flat_len {
960 crate::bail_invalid_estim!(
961 "analytic penalty Hessian-vector length mismatch: got {}, expected {}",
962 hv.len(),
963 flat_len
964 );
965 }
966 for row in 0..flat_len {
967 hessian[[t_start + row, t_start + col]] += hv[row];
968 }
969 probe[col] = 0.0;
970 }
971 }
972 Ok(())
973}
974
975fn add_analytic_penalty_objective_to_eval(
976 theta: &Array1<f64>,
977 rho_dim: usize,
978 latent: &gam_terms::latent::LatentCoordValues,
979 registry: &gam_terms::AnalyticPenaltyRegistry,
980 eval: &mut (f64, Array1<f64>, gam_problem::HessianValue),
981) -> Result<(), EstimationError> {
982 let contribution = analytic_penalty_objective_contribution(theta, rho_dim, latent, registry)?;
983 eval.0 += contribution.cost;
984 if eval.1.len() != contribution.gradient.len() {
985 crate::bail_invalid_estim!(
986 "latent-coordinate REML gradient length mismatch: base={}, analytic_penalty={}",
987 eval.1.len(),
988 contribution.gradient.len()
989 );
990 }
991 eval.1 += &contribution.gradient;
992 add_analytic_penalty_hessian_to_eval(theta, rho_dim, latent, registry, eval)?;
993 Ok(())
994}
995
996fn spatial_log_kappa_hyper_dirs_frominfo_list(
997 info_list: Vec<SpatialPsiDerivative>,
998) -> Result<Vec<DirectionalHyperParam>, EstimationError> {
999 use gam_solve::estimate::reml::ImplicitDerivLevel;
1000 use std::collections::HashMap;
1001
1002 let log_kappa_dim = info_list.len();
1003 let group_ids: Vec<Option<usize>> = info_list.iter().map(|e| e.aniso_group_id).collect();
1009 let mut group_indices_map: HashMap<usize, Vec<usize>> = HashMap::new();
1010 for (idx, gid) in group_ids.iter().enumerate() {
1011 if let Some(g) = gid {
1012 group_indices_map.entry(*g).or_default().push(idx);
1013 }
1014 }
1015
1016 let mut hyper_dirs = Vec::with_capacity(log_kappa_dim);
1017 for (i, info) in info_list.into_iter().enumerate() {
1018 let SpatialPsiDerivative {
1019 penalty_index: _,
1020 penalty_indices,
1021 global_range,
1022 total_p,
1023 x_psi_local,
1024 s_psi_components_local,
1025 x_psi_psi_local,
1026 s_psi_psi_components_local,
1027 aniso_group_id,
1028 aniso_cross_designs,
1029 aniso_cross_penalty_provider,
1030 implicit_operator,
1031 implicit_axis,
1032 } = info;
1033
1034 let mut xsecond = vec![None; log_kappa_dim];
1035 xsecond[i] = Some(if let Some(ref op) = implicit_operator {
1037 gam_solve::estimate::reml::HyperDesignDerivative::from_implicit(
1038 op.clone(),
1039 ImplicitDerivLevel::SecondDiag(implicit_axis),
1040 global_range.clone(),
1041 total_p,
1042 )
1043 } else {
1044 gam_solve::estimate::reml::HyperDesignDerivative::from_embedded(
1045 x_psi_psi_local,
1046 global_range.clone(),
1047 total_p,
1048 )
1049 });
1050 if let Some(cross_designs) = aniso_cross_designs {
1052 if let Some(gid) = aniso_group_id {
1056 let base = group_indices_map
1057 .get(&gid)
1058 .and_then(|v| v.first().copied())
1059 .unwrap_or(i);
1060 for (b_axis, cross_mat) in cross_designs.into_iter() {
1061 let j = base + b_axis;
1062 if j < log_kappa_dim {
1063 xsecond[j] = Some(if let Some(ref op) = implicit_operator {
1064 gam_solve::estimate::reml::HyperDesignDerivative::from_implicit(
1065 op.clone(),
1066 ImplicitDerivLevel::SecondCross(implicit_axis, b_axis),
1067 global_range.clone(),
1068 total_p,
1069 )
1070 } else {
1071 gam_solve::estimate::reml::HyperDesignDerivative::from_embedded(
1072 cross_mat,
1073 global_range.clone(),
1074 total_p,
1075 )
1076 });
1077 }
1078 }
1079 }
1080 }
1081 let s_components = penalty_indices
1082 .iter()
1083 .copied()
1084 .zip(s_psi_components_local.into_iter().map(|local| {
1085 gam_solve::estimate::reml::HyperPenaltyDerivative::from_embedded(
1086 local,
1087 global_range.clone(),
1088 total_p,
1089 )
1090 }))
1091 .collect::<Vec<_>>();
1092 let s2_components = penalty_indices
1093 .iter()
1094 .copied()
1095 .zip(s_psi_psi_components_local.into_iter().map(|local| {
1096 gam_solve::estimate::reml::HyperPenaltyDerivative::from_embedded(
1097 local,
1098 global_range.clone(),
1099 total_p,
1100 )
1101 }))
1102 .collect::<Vec<_>>();
1103 let mut ssecond_components = vec![None; log_kappa_dim];
1104 ssecond_components[i] = Some(s2_components);
1105 let mut penaltysecond_partner_indices: Option<Vec<usize>> = None;
1106 let penaltysecond_component_provider =
1107 if let (Some(provider), Some(gid)) = (aniso_cross_penalty_provider, aniso_group_id) {
1108 let group_indices = group_indices_map.get(&gid).cloned().unwrap_or_default();
1109 let axis_in_group =
1110 group_indices
1111 .iter()
1112 .position(|&idx| idx == i)
1113 .ok_or_else(|| {
1114 EstimationError::InvalidInput(format!(
1115 "missing spatial hyper axis {} in anisotropy group {}",
1116 i, gid
1117 ))
1118 })?;
1119 penaltysecond_partner_indices = Some(
1120 group_indices
1121 .iter()
1122 .copied()
1123 .filter(|&idx| idx != i)
1124 .collect(),
1125 );
1126 let penalty_indices_inner = penalty_indices.clone();
1127 let global_range_inner = global_range.clone();
1128 let total_p_inner = total_p;
1129 let group_indices_inner = group_indices;
1130 Some(std::sync::Arc::new(
1131 move |j: usize| -> Result<
1132 Option<Vec<gam_solve::estimate::reml::PenaltyDerivativeComponent>>,
1133 EstimationError,
1134 > {
1135 let Some(other_axis_in_group) =
1136 group_indices_inner.iter().position(|&idx| idx == j)
1137 else {
1138 return Ok(None);
1139 };
1140 if other_axis_in_group == axis_in_group {
1141 return Ok(None);
1142 }
1143 let cross_pens = provider(other_axis_in_group)?;
1144 if cross_pens.is_empty() {
1145 return Ok(None);
1146 }
1147 Ok(Some(
1148 penalty_indices_inner
1149 .iter()
1150 .copied()
1151 .zip(cross_pens.into_iter().map(|local| {
1152 gam_solve::estimate::reml::HyperPenaltyDerivative::from_embedded(
1153 local,
1154 global_range_inner.clone(),
1155 total_p_inner,
1156 )
1157 }))
1158 .map(|(penalty_index, matrix)| {
1159 gam_solve::estimate::reml::PenaltyDerivativeComponent {
1160 penalty_index,
1161 matrix,
1162 }
1163 })
1164 .collect(),
1165 ))
1166 },
1167 )
1168 as std::sync::Arc<
1169 dyn Fn(
1170 usize,
1171 ) -> Result<
1172 Option<Vec<gam_solve::estimate::reml::PenaltyDerivativeComponent>>,
1173 EstimationError,
1174 > + Send
1175 + Sync
1176 + 'static,
1177 >)
1178 } else {
1179 None
1180 };
1181 let x_first_hyper = if let Some(ref op) = implicit_operator {
1184 gam_solve::estimate::reml::HyperDesignDerivative::from_implicit(
1185 op.clone(),
1186 ImplicitDerivLevel::First(implicit_axis),
1187 global_range.clone(),
1188 total_p,
1189 )
1190 } else {
1191 gam_solve::estimate::reml::HyperDesignDerivative::from_embedded(
1192 x_psi_local,
1193 global_range.clone(),
1194 total_p,
1195 )
1196 };
1197 let mut dir = DirectionalHyperParam::new_compact(
1198 x_first_hyper,
1199 s_components,
1200 Some(xsecond),
1201 Some(ssecond_components),
1202 )?
1203 .not_penalty_like();
1204 if let Some(provider) = penaltysecond_component_provider {
1205 dir = dir.with_penaltysecond_component_provider(provider);
1206 }
1207 if let Some(partner_indices) = penaltysecond_partner_indices {
1208 dir = dir.with_penaltysecond_partner_indices(partner_indices);
1209 }
1210 hyper_dirs.push(dir);
1211 }
1212 Ok(hyper_dirs)
1213}
1214
1215pub(crate) fn spatial_dims_per_term(
1221 resolvedspec: &TermCollectionSpec,
1222 spatial_terms: &[usize],
1223) -> Vec<usize> {
1224 spatial_terms
1225 .iter()
1226 .map(|&term_idx| {
1227 if let Some(mj) = measure_jet_term_spec(resolvedspec, term_idx) {
1228 measure_jet_psi_dim(mj)
1231 } else if spatial_term_uses_per_axis_psi(resolvedspec, term_idx) {
1232 get_spatial_feature_dim(resolvedspec, term_idx).unwrap_or(1)
1233 } else {
1234 1
1235 }
1236 })
1237 .collect()
1238}
1239
1240fn has_aniso_terms(resolvedspec: &TermCollectionSpec, spatial_terms: &[usize]) -> bool {
1244 spatial_terms
1245 .iter()
1246 .any(|&term_idx| spatial_term_uses_per_axis_psi(resolvedspec, term_idx))
1247}
1248
1249macro_rules! impl_exact_joint_theta_memo {
1255 () => {
1256 fn memoized_cost(&self, theta: &Array1<f64>) -> Option<f64> {
1257 if self
1258 .current_theta
1259 .as_ref()
1260 .is_some_and(|cached| theta_values_match(cached, theta))
1261 {
1262 self.last_eval
1263 .as_ref()
1264 .map(|cached| cached.0)
1265 .or(self.last_cost)
1266 } else {
1267 None
1268 }
1269 }
1270
1271 fn memoized_eval(
1272 &self,
1273 theta: &Array1<f64>,
1274 ) -> Option<(f64, Array1<f64>, gam_problem::HessianValue)> {
1275 if self
1276 .current_theta
1277 .as_ref()
1278 .is_some_and(|cached| theta_values_match(cached, theta))
1279 {
1280 self.last_eval.clone()
1281 } else {
1282 None
1283 }
1284 }
1285
1286 fn store_eval(&mut self, eval: (f64, Array1<f64>, gam_problem::HessianValue)) {
1287 self.last_cost = Some(eval.0);
1288 self.last_eval = Some(eval);
1289 }
1290 };
1291}
1292
1293struct SingleBlockExactJointDesignCache<'d> {
1294 realizer: FrozenTermCollectionIncrementalRealizer<'d>,
1295 current_theta: Option<Array1<f64>>,
1296 last_eval_theta: Option<Array1<f64>>,
1303 last_cost: Option<f64>,
1304 last_eval: Option<(f64, Array1<f64>, gam_problem::HessianValue)>,
1305 cached_hyper_dirs: Option<(u64, Vec<DirectionalHyperParam>)>,
1317 spatial_terms: Vec<usize>,
1318 rho_dim: usize,
1319 dims_per_term: Vec<usize>,
1320}
1321
1322impl<'d> SingleBlockExactJointDesignCache<'d> {
1323 fn new_with_policy(
1324 data: ArrayView2<'d, f64>,
1325 spec: TermCollectionSpec,
1326 design: TermCollectionDesign,
1327 spatial_terms: Vec<usize>,
1328 rho_dim: usize,
1329 dims_per_term: Vec<usize>,
1330 policy: &gam_runtime::resource::ResourcePolicy,
1331 ) -> Result<Self, String> {
1332 Ok(Self {
1333 realizer: FrozenTermCollectionIncrementalRealizer::new_with_policy(
1334 data, spec, design, policy,
1335 )?,
1336 current_theta: None,
1337 last_eval_theta: None,
1338 last_cost: None,
1339 last_eval: None,
1340 cached_hyper_dirs: None,
1341 spatial_terms,
1342 rho_dim,
1343 dims_per_term,
1344 })
1345 }
1346
1347 fn design_revision(&self) -> u64 {
1348 self.realizer.design_revision()
1349 }
1350
1351 fn hyper_dirs_for_current_design(
1361 &mut self,
1362 data: ArrayView2<'_, f64>,
1363 kind: SpatialHyperKind,
1364 ) -> Result<Vec<DirectionalHyperParam>, EstimationError> {
1365 let revision = self.realizer.design_revision();
1366 if let Some((cached_rev, dirs)) = self.cached_hyper_dirs.as_ref()
1367 && *cached_rev == revision
1368 {
1369 return Ok(dirs.clone());
1370 }
1371 let dirs = try_build_spatial_log_kappa_hyper_dirs(
1372 data,
1373 self.realizer.spec(),
1374 self.realizer.design(),
1375 &self.spatial_terms,
1376 )?
1377 .ok_or_else(|| {
1378 EstimationError::InvalidInput(format!(
1379 "failed to build {} hyper_dirs at current {}",
1380 kind.adjective(),
1381 kind.coord_name(),
1382 ))
1383 })?;
1384 self.cached_hyper_dirs = Some((revision, dirs.clone()));
1385 Ok(dirs)
1386 }
1387
1388 fn nfree_tensor_gradient_hyper_dirs(
1389 &mut self,
1390 theta: &Array1<f64>,
1391 ) -> Result<Vec<DirectionalHyperParam>, EstimationError> {
1392 let psi = &theta.as_slice().ok_or_else(|| {
1393 EstimationError::InvalidInput(
1394 "nfree_tensor_gradient_hyper_dirs: theta is not contiguous".to_string(),
1395 )
1396 })?[self.rho_dim..];
1397 let (global_range, p_total, s_psi_components) = self
1398 .realizer
1399 .canonical_penalty_derivatives_at_psi(&self.spatial_terms, psi)
1400 .map_err(EstimationError::InvalidInput)?;
1401 let zero_x = gam_solve::estimate::reml::HyperDesignDerivative::zero(
1402 self.realizer.design().design.nrows(),
1403 p_total,
1404 );
1405 let components = s_psi_components
1406 .into_iter()
1407 .enumerate()
1408 .map(|(penalty_index, local)| {
1409 (
1410 penalty_index,
1411 gam_solve::estimate::reml::HyperPenaltyDerivative::from_embedded(
1412 local,
1413 global_range.clone(),
1414 p_total,
1415 ),
1416 )
1417 })
1418 .collect::<Vec<_>>();
1419 Ok(DirectionalHyperParam::new_compact(zero_x, components, None, None)?.not_penalty_like())
1420 .map(|dir| vec![dir])
1421 }
1422
1423 fn ensure_theta(&mut self, theta: &Array1<f64>) -> Result<(), EstimationError> {
1429 if self
1430 .current_theta
1431 .as_ref()
1432 .is_some_and(|cached| theta_values_match(cached, theta))
1433 {
1434 return Ok(());
1435 }
1436 let t_ensure = std::time::Instant::now();
1437 let log_kappa = SpatialLogKappaCoords::from_theta_tail_with_dims(
1438 theta,
1439 self.rho_dim,
1440 self.dims_per_term.clone(),
1441 );
1442 self.realizer
1443 .apply_log_kappa(&log_kappa, &self.spatial_terms)?;
1444 log::info!(
1445 "[STAGE] ensure_theta (apply_log_kappa, {} terms): {:.3}s",
1446 self.spatial_terms.len(),
1447 t_ensure.elapsed().as_secs_f64(),
1448 );
1449 self.current_theta = Some(theta.clone());
1450 self.last_eval_theta = None;
1451 self.last_cost = None;
1452 self.last_eval = None;
1453 Ok(())
1454 }
1455
1456 fn memoized_cost(&self, theta: &Array1<f64>) -> Option<f64> {
1463 if self
1464 .last_eval_theta
1465 .as_ref()
1466 .is_some_and(|cached| theta_values_match(cached, theta))
1467 {
1468 self.last_eval
1469 .as_ref()
1470 .map(|cached| cached.0)
1471 .or(self.last_cost)
1472 } else {
1473 None
1474 }
1475 }
1476
1477 fn memoized_eval(
1478 &self,
1479 theta: &Array1<f64>,
1480 ) -> Option<(f64, Array1<f64>, gam_problem::HessianValue)> {
1481 if self
1482 .last_eval_theta
1483 .as_ref()
1484 .is_some_and(|cached| theta_values_match(cached, theta))
1485 {
1486 self.last_eval.clone()
1487 } else {
1488 None
1489 }
1490 }
1491
1492 fn forget_eval_memo(&mut self) {
1501 self.last_eval_theta = None;
1502 self.last_cost = None;
1503 self.last_eval = None;
1504 }
1505
1506 fn store_eval_at(
1510 &mut self,
1511 theta: &Array1<f64>,
1512 eval: (f64, Array1<f64>, gam_problem::HessianValue),
1513 ) {
1514 self.last_eval_theta = Some(theta.clone());
1515 self.last_cost = Some(eval.0);
1516 self.last_eval = Some(eval);
1517 }
1518
1519 fn store_cost_at(&mut self, theta: &Array1<f64>, cost: f64) {
1522 self.last_eval_theta = Some(theta.clone());
1523 self.last_cost = Some(cost);
1524 self.last_eval = None;
1528 }
1529
1530 fn spec(&self) -> &TermCollectionSpec {
1531 self.realizer.spec()
1532 }
1533
1534 fn design(&self) -> &TermCollectionDesign {
1535 self.realizer.design()
1536 }
1537
1538 fn supports_nfree_penalty_rekey(&self) -> bool {
1544 self.realizer
1545 .supports_nfree_penalty_rekey(&self.spatial_terms)
1546 }
1547
1548 fn supports_nfree_gradient_only_routing(&self) -> bool {
1549 self.realizer
1550 .supports_nfree_gradient_only_routing(&self.spatial_terms)
1551 }
1552
1553 fn canonical_penalties_at(
1563 &mut self,
1564 theta: &Array1<f64>,
1565 ) -> Result<(Vec<gam_terms::construction::CanonicalPenalty>, Vec<usize>), String> {
1566 let psi = &theta
1567 .as_slice()
1568 .ok_or_else(|| "canonical_penalties_at: theta is not contiguous".to_string())?
1569 [self.rho_dim..];
1570 self.realizer
1571 .canonical_penalties_at_psi(&self.spatial_terms, psi)
1572 }
1573}
1574
1575struct SingleBlockLatentCoordDesignCache {
1576 data: Array2<f64>,
1577 spec: TermCollectionSpec,
1578 design: TermCollectionDesign,
1579 current_theta: Option<Array1<f64>>,
1580 current_latent: Option<std::sync::Arc<gam_terms::latent::LatentCoordValues>>,
1581 current_hyper_dirs: Option<Vec<gam_solve::estimate::reml::DirectionalHyperParam>>,
1582 current_design_cache_id: Option<u64>,
1583 latent_design_cache: gam_solve::latent_cache::LatentDesignCache,
1584 last_cost: Option<f64>,
1585 last_eval: Option<(f64, Array1<f64>, gam_problem::HessianValue)>,
1586 term_index: gam_problem::types::SmoothTermIdx,
1587 feature_cols: Vec<usize>,
1588 rho_dim: usize,
1589 n_obs: usize,
1590 latent_dim: usize,
1591 id_mode: gam_terms::latent::LatentIdMode,
1592 manifold: gam_terms::latent::LatentManifold,
1593 retraction_registry: gam_solve::latent_cache::LatentRetractionRegistry,
1594 latent_id: u64,
1595 analytic_penalties: Option<std::sync::Arc<gam_terms::AnalyticPenaltyRegistry>>,
1596 analytic_rho_count: usize,
1597 design_revision: u64,
1598}
1599
1600impl SingleBlockLatentCoordDesignCache {
1601 fn new(
1602 data: Array2<f64>,
1603 spec: TermCollectionSpec,
1604 design: TermCollectionDesign,
1605 latent: &StandardLatentCoordConfig,
1606 rho_dim: usize,
1607 ) -> Result<Self, String> {
1608 if latent.term_index.get() >= spec.smooth_terms.len() {
1609 return Err(SmoothError::dimension_mismatch(format!(
1610 "latent-coordinate term index {} out of bounds for {} smooth terms",
1611 latent.term_index,
1612 spec.smooth_terms.len()
1613 ))
1614 .into());
1615 }
1616 if latent.feature_cols.len() != latent.values.latent_dim() {
1617 return Err(SmoothError::dimension_mismatch(format!(
1618 "latent-coordinate feature width mismatch: feature_cols={}, latent_dim={}",
1619 latent.feature_cols.len(),
1620 latent.values.latent_dim()
1621 ))
1622 .into());
1623 }
1624 if latent.values.n_obs() != data.nrows() {
1625 return Err(SmoothError::dimension_mismatch(format!(
1626 "latent-coordinate row mismatch: latent n={}, data n={}",
1627 latent.values.n_obs(),
1628 data.nrows()
1629 ))
1630 .into());
1631 }
1632 let analytic_rho_count = latent
1633 .analytic_penalties
1634 .as_ref()
1635 .map_or(0, |registry| registry.total_rho_count());
1636 Ok(Self {
1637 data,
1638 spec,
1639 design,
1640 current_theta: None,
1641 current_latent: None,
1642 current_hyper_dirs: None,
1643 current_design_cache_id: None,
1644 latent_design_cache: gam_solve::latent_cache::LatentDesignCache::default(),
1645 last_cost: None,
1646 last_eval: None,
1647 term_index: latent.term_index,
1648 feature_cols: latent.feature_cols.clone(),
1649 rho_dim,
1650 n_obs: latent.values.n_obs(),
1651 latent_dim: latent.values.latent_dim(),
1652 id_mode: latent.values.id_mode().clone(),
1653 manifold: latent.values.manifold().clone(),
1654 retraction_registry: latent.values.retraction_registry().clone(),
1655 latent_id: latent.values.latent_id(),
1656 analytic_penalties: latent.analytic_penalties.clone(),
1657 analytic_rho_count,
1658 design_revision: 0,
1659 })
1660 }
1661
1662 fn design_revision(&self) -> u64 {
1663 self.design_revision
1664 }
1665
1666 fn design(&self) -> &TermCollectionDesign {
1667 &self.design
1668 }
1669
1670 fn latent(&self) -> Result<std::sync::Arc<gam_terms::latent::LatentCoordValues>, String> {
1671 self.current_latent
1672 .as_ref()
1673 .cloned()
1674 .ok_or_else(|| "latent-coordinate cache has not been realized".to_string())
1675 }
1676
1677 fn analytic_penalties(&self) -> Option<std::sync::Arc<gam_terms::AnalyticPenaltyRegistry>> {
1678 self.analytic_penalties.clone()
1679 }
1680
1681 fn analytic_penalty_rho_count(&self) -> usize {
1682 self.analytic_rho_count
1683 }
1684
1685 fn hyper_dirs(&self) -> Result<Vec<gam_solve::estimate::reml::DirectionalHyperParam>, String> {
1686 self.current_hyper_dirs
1687 .as_ref()
1688 .cloned()
1689 .ok_or_else(|| "latent-coordinate hyper_dirs cache has not been realized".to_string())
1690 }
1691
1692 fn latent_basis_kind(&self) -> Result<gam_solve::latent_cache::LatentBasisKind, String> {
1693 let smooth_term = self
1694 .design
1695 .smooth
1696 .terms
1697 .get(self.term_index.get())
1698 .ok_or_else(|| {
1699 SmoothError::dimension_mismatch(format!(
1700 "LatentCoord term index {} out of bounds for realized smooth design",
1701 self.term_index
1702 ))
1703 })?;
1704 let termspec = self
1705 .spec
1706 .smooth_terms
1707 .get(self.term_index.get())
1708 .ok_or_else(|| {
1709 SmoothError::dimension_mismatch(format!(
1710 "LatentCoord term index {} out of bounds for resolved smooth spec",
1711 self.term_index
1712 ))
1713 })?;
1714 match (&termspec.basis, &smooth_term.metadata) {
1715 (
1716 SmoothBasisSpec::Matern { .. },
1717 BasisMetadata::Matern {
1718 centers,
1719 length_scale,
1720 nu,
1721 aniso_log_scales,
1722 input_scale,
1723 ..
1724 },
1725 ) => Ok(gam_solve::latent_cache::LatentBasisKind::Matern {
1726 centers: centers.clone(),
1727 input_scale: *input_scale,
1730 length_scale: *length_scale,
1731 nu: *nu,
1732 aniso_log_scales: aniso_log_scales
1733 .clone()
1734 .unwrap_or_else(|| vec![0.0; centers.ncols()]),
1735 chunk_size: gam_terms::basis::auto_streaming_chunk_size_for_dense(
1736 self.n_obs,
1737 centers.nrows(),
1738 ),
1739 }),
1740 (
1741 SmoothBasisSpec::Duchon { .. },
1742 BasisMetadata::Duchon {
1743 centers,
1744 length_scale,
1745 power,
1746 nullspace_order,
1747 aniso_log_scales,
1748 input_scale,
1749 ..
1750 },
1751 ) => Ok(gam_solve::latent_cache::LatentBasisKind::Duchon {
1752 centers: centers.clone(),
1753 input_scale: *input_scale,
1755 length_scale: *length_scale,
1756 power: *power,
1757 nullspace_order: *nullspace_order,
1758 aniso_log_scales: aniso_log_scales
1759 .clone()
1760 .unwrap_or_else(|| vec![0.0; centers.ncols()]),
1761 }),
1762 (
1763 SmoothBasisSpec::Sphere { .. },
1764 BasisMetadata::Sphere {
1765 centers,
1766 penalty_order,
1767 method,
1768 ..
1769 },
1770 ) if matches!(*method, gam_terms::basis::SphereMethod::Wahba) => {
1771 Ok(gam_solve::latent_cache::LatentBasisKind::Sphere {
1772 centers: centers.clone(),
1773 penalty_order: *penalty_order,
1774 chunk_size: gam_terms::basis::auto_streaming_chunk_size_for_dense(
1775 self.n_obs,
1776 centers.nrows(),
1777 ),
1778 })
1779 }
1780 (
1781 SmoothBasisSpec::BSpline1D { spec, .. },
1782 BasisMetadata::BSpline1D {
1783 knots,
1784 periodic,
1785 degree: meta_degree,
1786 ..
1787 },
1788 ) => {
1789 let effective_degree = meta_degree.unwrap_or(spec.degree);
1793 if let Some((domain_start, period, num_basis)) = periodic {
1794 Ok(gam_solve::latent_cache::LatentBasisKind::PeriodicBspline {
1795 domain_start: *domain_start,
1796 period: *period,
1797 degree: effective_degree,
1798 num_basis: *num_basis,
1799 chunk_size: gam_terms::basis::auto_streaming_chunk_size_for_dense(
1800 self.n_obs, *num_basis,
1801 ),
1802 })
1803 } else {
1804 let num_basis_est = knots.len().saturating_sub(effective_degree + 1);
1805 Ok(gam_solve::latent_cache::LatentBasisKind::TensorBspline {
1806 knots: vec![knots.clone()],
1807 degrees: vec![effective_degree],
1808 chunk_size: gam_terms::basis::auto_streaming_chunk_size_for_dense(
1809 self.n_obs,
1810 num_basis_est,
1811 ),
1812 })
1813 }
1814 }
1815 (
1816 SmoothBasisSpec::TensorBSpline { .. },
1817 BasisMetadata::TensorBSpline { knots, degrees, .. },
1818 ) => Ok(gam_solve::latent_cache::LatentBasisKind::TensorBspline {
1819 knots: knots.clone(),
1820 degrees: degrees.clone(),
1821 chunk_size: None,
1822 }),
1823 (
1824 SmoothBasisSpec::Pca { .. },
1825 BasisMetadata::Pca {
1826 basis_matrix,
1827 centered,
1828 smooth_penalty,
1829 center_mean,
1830 pca_basis_path,
1831 chunk_size,
1832 ..
1833 },
1834 ) => {
1835 let center_mean_fingerprint = if *centered && pca_basis_path.is_none() {
1836 let mean = center_mean.as_ref().ok_or_else(|| {
1837 SmoothError::invalid_config(
1838 "latent-coordinate Pca cache key requires center_mean when centered",
1839 )
1840 })?;
1841 Some(gam_solve::latent_cache::pca_center_mean_fingerprint(mean))
1842 } else {
1843 None
1844 };
1845 Ok(gam_solve::latent_cache::LatentBasisKind::Pca {
1846 basis_matrix: basis_matrix.clone(),
1847 centered: *centered,
1848 center_mean_fingerprint,
1849 smooth_penalty: *smooth_penalty,
1850 pca_basis_path: pca_basis_path.clone(),
1851 chunk_size: *chunk_size,
1852 })
1853 }
1854 _ => Err(SmoothError::invalid_config(
1855 "latent-coordinate design cache could not key the realized latent smooth basis"
1856 .to_string(),
1857 )
1858 .into()),
1859 }
1860 }
1861
1862 fn ensure_theta(&mut self, theta: &Array1<f64>) -> Result<(), String> {
1863 if self
1864 .current_theta
1865 .as_ref()
1866 .is_some_and(|cached| theta_values_match(cached, theta))
1867 {
1868 return Ok(());
1869 }
1870 let latent_flat_len = self.n_obs * self.latent_dim;
1871 let direct_hyper_count = latent_coord_direct_hyper_count(&self.id_mode, self.latent_dim);
1872 let expected =
1873 self.rho_dim + latent_flat_len + self.analytic_rho_count + direct_hyper_count;
1874 if theta.len() != expected {
1875 return Err(SmoothError::dimension_mismatch(format!(
1876 "latent-coordinate theta length mismatch: got {}, expected {} (rho_dim={}, n={}, d={}, analytic_rhos={}, direct_hypers={})",
1877 theta.len(),
1878 expected,
1879 self.rho_dim,
1880 self.n_obs,
1881 self.latent_dim,
1882 self.analytic_rho_count,
1883 direct_hyper_count
1884 ))
1885 .into());
1886 }
1887 let flat = theta
1888 .slice(s![self.rho_dim..self.rho_dim + latent_flat_len])
1889 .to_owned();
1890 let latent = std::sync::Arc::new(
1891 gam_terms::latent::LatentCoordValues::from_flat_with_manifold_and_retraction_and_id(
1892 flat,
1893 self.n_obs,
1894 self.latent_dim,
1895 self.id_mode.clone(),
1896 self.manifold.clone(),
1897 self.retraction_registry.clone(),
1898 self.latent_id,
1899 ),
1900 );
1901 let latent_values_changed = self
1902 .current_latent
1903 .as_ref()
1904 .map(|cached| !latent_values_match(cached.as_flat(), latent.as_flat()))
1905 .unwrap_or(true);
1906 if latent_values_changed {
1907 self.latent_design_cache.invalidate_all();
1908 self.current_design_cache_id = None;
1909 self.design_revision = self.design_revision.wrapping_add(1);
1910 }
1911 for n in 0..self.n_obs {
1912 for axis in 0..self.latent_dim {
1913 let col = self.feature_cols[axis];
1914 self.data[[n, col]] = latent.as_flat()[n * self.latent_dim + axis];
1915 }
1916 }
1917
1918 let basis_kind = self.latent_basis_kind()?;
1919 let rebuilt_width = self.design.design.ncols();
1920 let spec = self.spec.clone();
1921 let term_index = self.term_index;
1922 let analytic_rho_count = self.analytic_rho_count;
1923 let data = self.data.view();
1924 let design_context_digest = gam_solve::latent_cache::latent_design_context_cache_digest(
1925 data,
1926 &spec,
1927 term_index,
1928 analytic_rho_count,
1929 &self.feature_cols,
1930 )
1931 .map_err(|e| e.to_string())?;
1932 let lookup = self
1933 .latent_design_cache
1934 .lookup_or_compute(latent.clone(), basis_kind, design_context_digest, || {
1935 let rebuilt = build_term_collection_design(data, &spec).map_err(|e| {
1936 EstimationError::InvalidInput(format!(
1937 "failed to rebuild latent-coordinate design: {e}"
1938 ))
1939 })?;
1940 if rebuilt.design.ncols() != rebuilt_width {
1941 crate::bail_invalid_estim!(
1942 "latent-coordinate design topology changed: rebuilt p={}, cached p={}",
1943 rebuilt.design.ncols(),
1944 rebuilt_width
1945 );
1946 }
1947 let hyper_dirs = try_build_latent_coord_hyper_dirs(
1948 latent.clone(),
1949 &spec,
1950 &rebuilt,
1951 &[term_index],
1952 analytic_rho_count,
1953 )?
1954 .ok_or_else(|| {
1955 EstimationError::InvalidInput(
1956 "failed to build latent-coordinate hyper_dirs".to_string(),
1957 )
1958 })?;
1959 Ok(gam_solve::latent_cache::ComputedLatentDesign {
1960 design: rebuilt,
1961 hyper_dirs,
1962 })
1963 })
1964 .map_err(|e| e.to_string())?;
1965 if lookup.cached.design.design.ncols() != self.design.design.ncols() {
1966 return Err(SmoothError::dimension_mismatch(format!(
1967 "latent-coordinate design topology changed: rebuilt p={}, cached p={}",
1968 lookup.cached.design.design.ncols(),
1969 self.design.design.ncols()
1970 ))
1971 .into());
1972 }
1973 self.design = lookup.cached.design.clone();
1974 self.current_hyper_dirs = Some(lookup.cached.hyper_dirs.clone());
1975 self.current_latent = Some(latent);
1976 self.current_theta = Some(theta.clone());
1977 self.last_cost = None;
1978 self.last_eval = None;
1979 if !latent_values_changed && self.current_design_cache_id != Some(lookup.entry_id) {
1980 self.design_revision = self.design_revision.wrapping_add(1);
1981 }
1982 self.current_design_cache_id = Some(lookup.entry_id);
1983 Ok(())
1984 }
1985
1986 fn memoized_cost(&self, theta: &Array1<f64>) -> Option<f64> {
1987 if self
1988 .current_theta
1989 .as_ref()
1990 .is_some_and(|cached| theta_values_match(cached, theta))
1991 {
1992 self.last_eval
1993 .as_ref()
1994 .map(|cached| cached.0)
1995 .or(self.last_cost)
1996 } else {
1997 None
1998 }
1999 }
2000
2001 fn memoized_eval(
2002 &self,
2003 theta: &Array1<f64>,
2004 ) -> Option<(f64, Array1<f64>, gam_problem::HessianValue)> {
2005 if self
2006 .current_theta
2007 .as_ref()
2008 .is_some_and(|cached| theta_values_match(cached, theta))
2009 {
2010 self.last_eval.clone()
2011 } else {
2012 None
2013 }
2014 }
2015
2016 fn store_eval(&mut self, eval: (f64, Array1<f64>, gam_problem::HessianValue)) {
2017 self.last_cost = Some(eval.0);
2018 self.last_eval = Some(eval);
2019 }
2020
2021 fn store_cost(&mut self, cost: f64) {
2022 self.last_cost = Some(cost);
2023 }
2024
2025 fn reset(&mut self) {
2026 self.current_theta = None;
2027 self.current_latent = None;
2028 self.current_hyper_dirs = None;
2029 self.current_design_cache_id = None;
2030 self.latent_design_cache.invalidate();
2031 self.last_cost = None;
2032 self.last_eval = None;
2033 }
2034}
2035
2036pub fn fixed_kappa_profiled_reml_score(
2056 data: ArrayView2<'_, f64>,
2057 y: ArrayView1<'_, f64>,
2058 weights: ArrayView1<'_, f64>,
2059 offset: ArrayView1<'_, f64>,
2060 resolvedspec: &TermCollectionSpec,
2061 term_idx: usize,
2062 kappa: f64,
2063 family: LikelihoodSpec,
2064 options: &FitOptions,
2065) -> Result<f64, EstimationError> {
2066 if !kappa.is_finite() {
2067 crate::bail_invalid_estim!("fixed-κ profiled score probed a non-finite κ = {kappa}");
2068 }
2069 if y.len() != data.nrows() || weights.len() != data.nrows() || offset.len() != data.nrows() {
2070 crate::bail_invalid_estim!(
2071 "fixed-κ profiled score row mismatch: data={}, y={}, weights={}, offset={}",
2072 data.nrows(),
2073 y.len(),
2074 weights.len(),
2075 offset.len(),
2076 );
2077 }
2078 let mut probe_spec = resolvedspec.clone();
2085 match probe_spec
2086 .smooth_terms
2087 .get_mut(term_idx)
2088 .map(|t| &mut t.basis)
2089 {
2090 Some(SmoothBasisSpec::ConstantCurvature { spec, .. }) => spec.kappa = kappa,
2091 _ => {
2092 crate::bail_invalid_estim!(
2093 "fixed-κ profiled score: term {term_idx} is not a constant-curvature smooth"
2094 )
2095 }
2096 }
2097 let fixed_kappa_options = SpatialLengthScaleOptimizationOptions {
2098 enabled: false,
2099 ..SpatialLengthScaleOptimizationOptions::default()
2100 };
2101 let fit = fit_term_collectionwith_spatial_length_scale_optimization(
2102 data,
2103 y.to_owned(),
2104 weights.to_owned(),
2105 offset.to_owned(),
2106 &probe_spec,
2107 family,
2108 options,
2109 &fixed_kappa_options,
2110 )?;
2111 let Some(score) = fit.fit.reml_score() else {
2112 crate::bail_invalid_estim!(
2113 "fixed-κ profiled fit at κ={kappa} has no REML/LAML score to profile against"
2114 );
2115 };
2116 if !score.is_finite() {
2117 crate::bail_invalid_estim!(
2118 "fixed-κ profiled fit at κ={kappa} returned a non-finite REML/LAML score"
2119 );
2120 }
2121 Ok(score)
2122}
2123
2124pub(crate) const JOINT_RHO_BOUND: f64 = 12.0;
2132
2133pub(crate) fn joint_rho_search_box(
2195 rho_seed: ArrayView1<'_, f64>,
2196 rho_upper_bound: f64,
2197) -> (Array1<f64>, Array1<f64>) {
2198 let rho_dim = rho_seed.len();
2199 let lower = Array1::<f64>::from_shape_fn(rho_dim, |k| {
2200 let seed = rho_seed[k];
2201 if seed.is_finite() && seed <= -JOINT_RHO_BOUND {
2202 -gam_solve::estimate::RHO_BOUND
2203 } else {
2204 -JOINT_RHO_BOUND
2205 }
2206 });
2207 let upper = Array1::<f64>::from_shape_fn(rho_dim, |k| {
2208 let seed = rho_seed[k];
2209 if seed.is_finite() && seed >= rho_upper_bound {
2210 gam_solve::estimate::RHO_BOUND
2211 } else {
2212 rho_upper_bound
2213 }
2214 });
2215 (lower, upper)
2216}
2217
2218enum JointSpatialKappaOutcome {
2244 Optimized(Box<FittedTermCollectionWithSpec>),
2246 DeclinedKeepIncumbent {
2249 baseline_score: f64,
2250 optimized_score: f64,
2251 },
2252 Unavailable,
2254}
2255
2256fn try_exact_joint_spatial_length_scale_optimization(
2257 data: ArrayView2<'_, f64>,
2258 y: ArrayView1<'_, f64>,
2259 weights: ArrayView1<'_, f64>,
2260 offset: ArrayView1<'_, f64>,
2261 resolvedspec: &TermCollectionSpec,
2262 best: &FittedTermCollection,
2263 family: LikelihoodSpec,
2264 options: &FitOptions,
2265 kappa_options: &SpatialLengthScaleOptimizationOptions,
2266 spatial_terms: &[usize],
2267) -> Result<JointSpatialKappaOutcome, EstimationError> {
2268 if spatial_terms.is_empty() {
2269 return Ok(JointSpatialKappaOutcome::Unavailable);
2270 }
2271 kappa_options
2276 .validate()
2277 .map_err(EstimationError::InvalidInput)?;
2278
2279 if try_build_spatial_log_kappa_hyper_dirs(data, resolvedspec, &best.design, spatial_terms)?
2280 .is_none()
2281 {
2282 if !constant_curvature_term_indices(resolvedspec).is_empty() {
2283 log::info!(
2284 "[#1464-trace] try_exact_joint RETURNED None (hyper_dirs unavailable); \
2285 κ̂ comes from a NON-joint path"
2286 );
2287 }
2288 return Ok(JointSpatialKappaOutcome::Unavailable);
2289 }
2290 if !constant_curvature_term_indices(resolvedspec).is_empty() {
2291 log::info!(
2292 "[#1464-trace] try_exact_joint ENTERED for {} spatial term(s); CC present",
2293 spatial_terms.len()
2294 );
2295 }
2296
2297 let rho_dim = best.fit.lambdas.len();
2298
2299 let has_constant_curvature_term = !constant_curvature_term_indices(resolvedspec).is_empty();
2313 let rho_upper_bound = if has_constant_curvature_term {
2314 gam_solve::estimate::RHO_BOUND
2315 } else {
2316 JOINT_RHO_BOUND
2317 };
2318
2319 let dims_per_term = spatial_dims_per_term(resolvedspec, spatial_terms);
2321 let use_aniso = has_aniso_terms(resolvedspec, spatial_terms);
2322
2323 let log_kappa0 = if use_aniso {
2328 SpatialLogKappaCoords::from_length_scales_aniso(resolvedspec, spatial_terms, kappa_options)
2329 } else {
2330 SpatialLogKappaCoords::from_length_scales(resolvedspec, spatial_terms, kappa_options)
2331 };
2332 let mut log_kappa0 = log_kappa0
2335 .reseed_from_data(data, resolvedspec, spatial_terms, kappa_options)
2336 .map_err(EstimationError::BasisError)?;
2337 let mut cc_profiled_values: Vec<(usize, f64)> = Vec::new();
2342 if has_constant_curvature_term {
2343 for (slot, &term_idx) in spatial_terms.iter().enumerate() {
2344 if constant_curvature_term_spec(resolvedspec, term_idx).is_none() {
2345 continue;
2346 }
2347 let kappa = get_constant_curvature_kappa(resolvedspec, term_idx)
2348 .expect("constant-curvature term exposes its kappa");
2349 log_kappa0.set_scalar_slot(slot, kappa);
2350 cc_profiled_values.push((slot, kappa));
2351 }
2352 }
2353 let log_kappa_lower = if use_aniso {
2354 SpatialLogKappaCoords::lower_bounds_aniso_from_data(
2355 data,
2356 resolvedspec,
2357 spatial_terms,
2358 &dims_per_term,
2359 kappa_options,
2360 )
2361 } else {
2362 SpatialLogKappaCoords::lower_bounds_from_data(
2363 data,
2364 resolvedspec,
2365 spatial_terms,
2366 kappa_options,
2367 )
2368 }
2369 .map_err(EstimationError::BasisError)?;
2370 let log_kappa_upper = if use_aniso {
2371 SpatialLogKappaCoords::upper_bounds_aniso_from_data(
2372 data,
2373 resolvedspec,
2374 spatial_terms,
2375 &dims_per_term,
2376 kappa_options,
2377 )
2378 } else {
2379 SpatialLogKappaCoords::upper_bounds_from_data(
2380 data,
2381 resolvedspec,
2382 spatial_terms,
2383 kappa_options,
2384 )
2385 }
2386 .map_err(EstimationError::BasisError)?;
2387 let mut log_kappa_lower = log_kappa_lower;
2388 let mut log_kappa_upper = log_kappa_upper;
2389 for &(slot, kappa) in &cc_profiled_values {
2390 log_kappa_lower.set_scalar_slot(slot, kappa);
2391 log_kappa_upper.set_scalar_slot(slot, kappa);
2392 log::info!("[spatial-kappa] slot {slot}: profiling rho at certified kappa={kappa}");
2393 }
2394 let log_kappa0 = log_kappa0.clamp_to_bounds(&log_kappa_lower, &log_kappa_upper);
2400
2401 for (slot, &term_idx) in spatial_terms.iter().enumerate() {
2414 if constant_curvature_term_spec(resolvedspec, term_idx).is_some()
2415 || measure_jet_term_spec(resolvedspec, term_idx).is_some()
2416 {
2417 continue;
2418 }
2419 let Some(incumbent) = get_spatial_length_scale(resolvedspec, term_idx) else {
2420 continue;
2423 };
2424 if !(incumbent.is_finite() && incumbent > 0.0) {
2425 continue;
2426 }
2427 let psi_incumbent = -incumbent.ln();
2428 let axes = log_kappa0.term_slice(slot);
2429 if axes.is_empty() {
2430 continue;
2431 }
2432 let psi_bar = axes.iter().sum::<f64>() / axes.len() as f64;
2433 let max_abs_axis = axes.iter().fold(0.0_f64, |acc, v| acc.max(v.abs()));
2441 let mean_roundoff =
2442 (axes.len() as f64 + 1.0) * f64::EPSILON * (max_abs_axis + psi_incumbent.abs());
2443 if (psi_bar - psi_incumbent).abs() > mean_roundoff {
2444 return Err(EstimationError::RemlOptimizationFailed(format!(
2445 "exact joint spatial optimization would grade its criterion at a psi the \
2446 scalar-rho incumbent was never realized at (term {term_idx}): \
2447 seed_psi_bar={psi_bar:.17e}, incumbent_psi={psi_incumbent:.17e}, \
2448 delta={:.6e}, incumbent_length_scale={incumbent:.17e}, \
2449 window=[{:.6e}, {:.6e}]. theta0 is not shared, so the monotonicity \
2450 certificate below would compare two different functions (#2726).",
2451 psi_bar - psi_incumbent,
2452 kappa_options.min_length_scale,
2453 kappa_options.max_length_scale,
2454 )));
2455 }
2456 }
2457
2458 let rho_seed = best.fit.lambdas.mapv(f64::ln);
2459 let (rho_lower, rho_upper) = joint_rho_search_box(rho_seed.view(), rho_upper_bound);
2460 let widened: Vec<usize> = (0..rho_dim)
2461 .filter(|&k| rho_lower[k] < -JOINT_RHO_BOUND || rho_upper[k] > rho_upper_bound)
2462 .collect();
2463 if !widened.is_empty() {
2464 log::info!(
2465 "[spatial-kappa] joint rho box fell back to the engine's own +/-RHO_BOUND on \
2466 coordinate(s) {widened:?}: their incumbent is not strictly inside the joint \
2467 +/-{JOINT_RHO_BOUND} prior, so the prior is falsified there and the search \
2468 region becomes the one the incumbent was found in (gam#2760). \
2469 seed={:?} box=[{:?}, {:?}]",
2470 rho_seed.to_vec(),
2471 rho_lower.to_vec(),
2472 rho_upper.to_vec(),
2473 );
2474 }
2475 let setup = ExactJointHyperSetup::new(rho_seed, rho_lower, rho_upper, log_kappa0, log_kappa_lower, log_kappa_upper);
2476
2477 let theta0 = setup.theta0();
2478 let lower = setup.lower();
2479 let upper = setup.upper();
2480
2481 let kind = if use_aniso {
2493 SpatialHyperKind::Anisotropic
2494 } else {
2495 SpatialHyperKind::Isotropic
2496 };
2497 let (theta_star, joint_final_value, joint_seed_value, kappa_timing) = run_exact_joint_spatial_optimization(
2498 kind,
2499 data,
2500 y,
2501 weights,
2502 offset,
2503 resolvedspec,
2504 &best.design,
2505 family.clone(),
2506 options,
2507 spatial_terms,
2508 &dims_per_term,
2509 &theta0,
2510 &lower,
2511 &upper,
2512 rho_dim,
2513 kappa_options,
2514 )?;
2515
2516 let baseline_score = fit_score(&best.fit);
2517
2518 let accept_tol = options.tol.max(1e-8 * baseline_score.abs()).max(1e-12);
2523 if !joint_seed_value.is_finite() {
2559 return Err(EstimationError::RemlOptimizationFailed(format!(
2560 "exact joint spatial optimization could not evaluate its own criterion at the \
2561 seed (seed_value={joint_seed_value:.6e}), so neither its descent nor its \
2562 agreement with the scalar-rho route is checkable; baseline={baseline_score:.6e}"
2563 )));
2564 }
2565 log::info!(
2576 "[spatial-kappa] route agreement at theta0: joint_seed={joint_seed_value:.12e} \
2577 baseline={baseline_score:.12e} gap={:.6e} ({:.6e} relative) \
2578 agreement_tolerance={accept_tol:.6e} ({}) sqrt_eps_scale={:.6e}",
2579 joint_seed_value - baseline_score,
2580 (joint_seed_value - baseline_score) / baseline_score.abs().max(f64::MIN_POSITIVE),
2581 if (joint_seed_value - baseline_score).abs() > accept_tol {
2582 "REFUSES"
2583 } else {
2584 "admits"
2585 },
2586 baseline_score.abs() * f64::EPSILON.sqrt(),
2587 );
2588 if (joint_seed_value - baseline_score).abs() > accept_tol {
2619 log::warn!(
2620 "[spatial-kappa] the joint and scalar-rho routes disagree about the criterion AT \
2621 THE SAME POINT theta0: joint_seed={joint_seed_value:.12e}, \
2622 baseline={baseline_score:.12e}, gap={:.3e} ({:.3e} relative) against a \
2623 {accept_tol:.3e} agreement tolerance. Two independent assemblies of one \
2624 criterion; their forward error is O(eps*kappa) in the penalized Hessian, so this \
2625 is only evidence of a formula difference when it exceeds what the conditioning \
2626 explains. The joint result is graded on the SHIPPED scalar-route score below, \
2627 which is a like-for-like comparison; this line is the record that the two \
2628 assemblies parted company (joint_final={joint_final_value:.12e}, \
2629 theta_checkpoint={:?}).",
2630 joint_seed_value - baseline_score,
2631 (joint_seed_value - baseline_score) / baseline_score.abs().max(f64::MIN_POSITIVE),
2632 theta_star.to_vec(),
2633 );
2634 }
2635 let (theta_star, joint_final_value) = if joint_final_value > joint_seed_value + accept_tol {
2656 log::warn!(
2657 "[spatial-kappa] the exact joint search terminated ABOVE its own seed \
2658 (seed={joint_seed_value:.12e}, final={joint_final_value:.12e}, \
2659 regression={:.3e}, acceptance_tolerance={accept_tol:.3e}); its terminal \
2660 certificate is local/boundary at theta={:?} and does not dominate the seed, \
2661 so the seed is kept and joint kappa optimization is a no-op for this fit. \
2662 A descent method returning a point worse than its start is a solver defect \
2663 in its own right and this line is the record of it.",
2664 joint_final_value - joint_seed_value,
2665 theta_star.to_vec(),
2666 );
2667 (theta0.clone(), joint_seed_value)
2668 } else {
2669 (theta_star, joint_final_value)
2670 };
2671
2672 let selected_lambdas = Array1::from_vec(
2673 gam_problem::checked_exp_log_strengths(
2674 theta_star.slice(s![..rho_dim]).iter().copied(),
2675 )
2676 .map_err(|error| {
2677 EstimationError::InvalidInput(format!(
2678 "selected joint spatial smoothing coordinate is outside the canonical log-strength domain: {error}"
2679 ))
2680 })?,
2681 );
2682 let log_kappa_star =
2683 SpatialLogKappaCoords::from_theta_tail_with_dims(&theta_star, rho_dim, dims_per_term);
2684 if has_constant_curvature_term {
2690 let star = log_kappa_star.as_array();
2691 let dims = log_kappa_star.dims_per_term();
2692 for (slot, &term_idx) in spatial_terms.iter().enumerate() {
2693 if constant_curvature_term_spec(resolvedspec, term_idx).is_some() {
2694 let off: usize = dims[..slot].iter().sum();
2695 log::info!(
2696 "[#1464-trace] term {term_idx}: joint solver CONVERGED ψ-tail κ = {} \
2697 (this is the optimised candidate; joint_final_value={joint_final_value})",
2698 star[off]
2699 );
2700 }
2701 }
2702 }
2703 let optimized_spec = log_kappa_star.apply_tospec(resolvedspec, spatial_terms)?;
2704 let optimized = fit_term_collection_forspecwith_heuristic_lambdas(
2705 data,
2706 y,
2707 weights,
2708 offset,
2709 &optimized_spec,
2710 selected_lambdas.as_slice(),
2711 family.clone(),
2712 options,
2713 )?;
2714
2715 let optimized_score = fit_score(&optimized.fit);
2726 if optimized_score > baseline_score + accept_tol {
2727 log::warn!(
2728 "[spatial-kappa] joint kappa optimization did not improve the SHIPPED scalar-route \
2729 score (baseline={baseline_score:.12e}, at theta_star={optimized_score:.12e}, \
2730 regression={:.3e}, acceptance_tolerance={accept_tol:.3e}); keeping the incumbent \
2731 fit and treating joint kappa optimization as a no-op for this fit. Both numbers \
2732 are `fit_score` of a scalar-route fit, so unlike the theta0 cross-route line this \
2733 comparison is like-for-like and a regression here is a real one.",
2734 optimized_score - baseline_score,
2735 );
2736 return Ok(JointSpatialKappaOutcome::DeclinedKeepIncumbent {
2737 baseline_score,
2738 optimized_score,
2739 });
2740 }
2741
2742 let mut fit = optimized.fit;
2746 fit.set_criterion(Some(joint_final_value));
2747 let optimized_result = FittedTermCollectionWithSpec {
2748 fit,
2749 design: optimized.design,
2750 resolvedspec: optimized_spec,
2751 adaptive_diagnostics: optimized.adaptive_diagnostics,
2752 kappa_timing: Some(kappa_timing),
2753 };
2754
2755 Ok(JointSpatialKappaOutcome::Optimized(Box::new(
2756 optimized_result,
2757 )))
2758}
2759
2760#[derive(Clone, Copy, PartialEq, Eq, Debug)]
2772enum SpatialHyperKind {
2773 Anisotropic,
2774 Isotropic,
2775}
2776
2777impl SpatialHyperKind {
2778 fn label(self) -> &'static str {
2781 match self {
2782 SpatialHyperKind::Anisotropic => "spatial-aniso-joint",
2783 SpatialHyperKind::Isotropic => "spatial-iso-joint",
2784 }
2785 }
2786
2787 fn adjective(self) -> &'static str {
2789 match self {
2790 SpatialHyperKind::Anisotropic => "anisotropic",
2791 SpatialHyperKind::Isotropic => "isotropic",
2792 }
2793 }
2794
2795 fn coord_name(self) -> &'static str {
2798 match self {
2799 SpatialHyperKind::Anisotropic => "psi",
2800 SpatialHyperKind::Isotropic => "kappa",
2801 }
2802 }
2803}
2804
2805struct SpatialFrozenGlmInputs {
2811 y: Array1<f64>,
2812 weights: Array1<f64>,
2813 offset: Array1<f64>,
2814 family: LikelihoodSpec,
2815}
2816
2817fn frozen_glm_tensor_eligible_family(family: &LikelihoodSpec) -> bool {
2834 !family.is_gaussian_identity()
2835 && matches!(
2836 &family.response,
2837 ResponseFamily::Binomial
2838 | ResponseFamily::Poisson
2839 | ResponseFamily::Gamma
2840 | ResponseFamily::NegativeBinomial { .. }
2841 )
2842}
2843
2844struct SpatialJointContext<'d> {
2845 data: ArrayView2<'d, f64>,
2846 rho_dim: usize,
2847 kind: SpatialHyperKind,
2848 cache: SingleBlockExactJointDesignCache<'d>,
2849 evaluator: gam_solve::estimate::ExternalJointHyperEvaluator<'d>,
2850 frozen_glm_inputs: Option<SpatialFrozenGlmInputs>,
2851 frozen_glm_psi_bounds: Option<(f64, f64)>,
2852 frozen_glm_tensor: Option<gam_solve::glm_sufficient_lane::FrozenWeightGramTensor>,
2853 frozen_glm_tensor_attempted: bool,
2854 frozen_glm_weight_memo: Option<(Array1<f64>, Array1<f64>)>,
2866 value_realization_failures: usize,
2872 value_evaluation_failures: usize,
2873 nfree_polish_boundary: Option<(u64, u64)>,
2888}
2889
2890#[derive(Clone, Copy, Debug, Default)]
2891struct NfreeSkipGateStatus {
2892 shape: bool,
2893 value: bool,
2894 gradient: bool,
2895 penalty: bool,
2896 revision: bool,
2897 second_order: bool,
2898}
2899
2900impl NfreeSkipGateStatus {
2901 fn would_skip(self, require_gradient: bool) -> bool {
2902 self.shape
2903 && self.value
2904 && (!require_gradient || self.gradient)
2905 && self.penalty
2906 && self.revision
2907 && !self.second_order
2908 }
2909}
2910
2911fn nfree_skip_gate_status_from_parts(
2912 shape: bool,
2913 covers_value: bool,
2914 covers_skip: bool,
2915 covers_gradient: bool,
2916 penalty: bool,
2917 revision: bool,
2918 allow_second_order: bool,
2919 require_gradient: bool,
2920) -> NfreeSkipGateStatus {
2921 NfreeSkipGateStatus {
2922 shape,
2923 value: shape && covers_value && (!require_gradient || covers_skip),
2931 gradient: shape && (!require_gradient || covers_gradient),
2932 penalty,
2933 revision,
2934 second_order: allow_second_order,
2935 }
2936}
2937
2938fn classify_spatial_value_probe_failure(
2943 error: EstimationError,
2944) -> Result<f64, EstimationError> {
2945 if is_recoverable_trial_point_error(&error) {
2946 Ok(f64::INFINITY)
2947 } else {
2948 Err(error)
2949 }
2950}
2951
2952impl<'d> SpatialJointContext<'d> {
2953 fn nfree_skip_gate_status(
2954 &self,
2955 theta: &Array1<f64>,
2956 allow_second_order: bool,
2957 require_gradient: bool,
2958 ) -> NfreeSkipGateStatus {
2959 let shape = theta.len() == self.rho_dim + 1;
2960 let (covers_value, covers_skip, covers_gradient) = if shape {
2961 let psi = theta[self.rho_dim];
2962 (
2963 self.evaluator.psi_gram_tensor_covers(psi),
2964 self.evaluator.psi_gram_tensor_covers_skip(psi),
2965 self.evaluator.psi_gram_tensor_covers_gradient(psi),
2966 )
2967 } else {
2968 (false, false, false)
2969 };
2970 nfree_skip_gate_status_from_parts(
2971 shape,
2972 covers_value,
2973 covers_skip,
2974 covers_gradient,
2975 self.evaluator.supports_nfree_penalty_rekey(),
2976 self.evaluator.nfree_fast_path_revision().is_some(),
2977 allow_second_order,
2978 require_gradient,
2979 )
2980 }
2981
2982 fn frozen_glm_working_state(
2983 &self,
2984 beta: &Array1<f64>,
2985 ) -> Result<Option<(Array1<f64>, Array1<f64>)>, EstimationError> {
2986 let Some(inputs) = self.frozen_glm_inputs.as_ref() else {
2987 return Ok(None);
2988 };
2989 if beta.len() != self.cache.design().design.ncols() {
2990 return Ok(None);
2991 }
2992 let mut eta = self.cache.design().design.matrixvectormultiply(beta);
2993 if eta.len() != inputs.offset.len() {
2994 crate::bail_invalid_estim!(
2995 "frozen GLM tensor warm-state row mismatch: eta={}, offset={}",
2996 eta.len(),
2997 inputs.offset.len()
2998 );
2999 }
3000 eta += &inputs.offset;
3001 let obs = evaluate_standard_familyobservations(
3002 inputs.family.clone(),
3003 None,
3004 None,
3005 None,
3006 &inputs.y,
3007 &inputs.weights,
3008 &eta,
3009 )?;
3010 let mut working_response = obs.eta.clone();
3011 for i in 0..working_response.len() {
3012 let wi = obs.fisherweight[i].max(1e-12);
3013 working_response[i] += obs.score[i] / wi;
3014 }
3015 Ok(Some((obs.fisherweight, working_response)))
3016 }
3017
3018 fn frozen_glm_trial_weights(
3027 &mut self,
3028 beta: &Array1<f64>,
3029 ) -> Result<Option<Array1<f64>>, EstimationError> {
3030 if let Some((memo_beta, memo_w)) = self.frozen_glm_weight_memo.as_ref()
3031 && memo_beta.len() == beta.len()
3032 && memo_beta
3033 .iter()
3034 .zip(beta.iter())
3035 .all(|(a, b)| a.to_bits() == b.to_bits())
3036 {
3037 return Ok(Some(memo_w.clone()));
3038 }
3039 match self.frozen_glm_working_state(beta)? {
3040 Some((current_w, _)) => {
3041 self.frozen_glm_weight_memo = Some((beta.clone(), current_w.clone()));
3042 Ok(Some(current_w))
3043 }
3044 None => Ok(None),
3045 }
3046 }
3047
3048 fn ensure_frozen_glm_tensor(
3049 &mut self,
3050 theta: &Array1<f64>,
3051 warm_beta: Option<&Array1<f64>>,
3052 ) -> Result<(), EstimationError> {
3053 if self.frozen_glm_tensor.is_some() || self.frozen_glm_tensor_attempted {
3054 return Ok(());
3055 }
3056 let Some((psi_lo, psi_hi)) = self.frozen_glm_psi_bounds else {
3057 return Ok(());
3058 };
3059 if theta.len() != self.rho_dim + 1 {
3060 self.frozen_glm_tensor_attempted = true;
3061 return Ok(());
3062 }
3063 let Some(beta) = warm_beta else {
3064 return Ok(());
3065 };
3066 let Some((frozen_w, working_z)) = self.frozen_glm_working_state(beta)? else {
3067 self.frozen_glm_tensor_attempted = true;
3068 return Ok(());
3069 };
3070 let theta_probe_base = theta.clone();
3071 let rho_dim = self.rho_dim;
3072 let Self {
3079 cache, evaluator, ..
3080 } = self;
3081 let tensor = evaluator.build_frozen_glm_gram_tensor(
3082 |psi| {
3083 let mut theta_probe = theta_probe_base.clone();
3084 theta_probe[rho_dim] = psi;
3085 cache.ensure_theta(&theta_probe).map_err(|e| e.to_string())?;
3086 Ok(cache.design().design.clone())
3087 },
3088 frozen_w.view(),
3089 working_z.view(),
3090 psi_lo,
3091 psi_hi,
3092 );
3093 self.cache
3094 .ensure_theta(theta)?;
3095 self.frozen_glm_tensor_attempted = true;
3096 if let Some(tensor) = tensor {
3097 self.frozen_glm_tensor = Some(tensor);
3098 log::info!(
3099 "[STAGE] {} certified frozen-W GLM ψ tensor over [{psi_lo:.3}, {psi_hi:.3}]",
3100 self.kind.label(),
3101 );
3102 } else {
3103 log::info!(
3104 "[STAGE] {} frozen-W GLM ψ tensor did not certify over [{psi_lo:.3}, {psi_hi:.3}]",
3105 self.kind.label(),
3106 );
3107 }
3108 Ok(())
3109 }
3110
3111 fn stage_frozen_glm_trial_statistics(
3112 &mut self,
3113 theta: &Array1<f64>,
3114 warm_beta: Option<&Array1<f64>>,
3115 allow_gradient: bool,
3116 ) -> Result<(), EstimationError> {
3117 let kind = self.kind;
3118 let mut staged_gram: Option<Array2<f64>> = None;
3119 let mut staged_deriv: Option<(Array2<f64>, Array1<f64>)> = None;
3120 if theta.len() == self.rho_dim + 1 {
3121 let psi = theta[self.rho_dim];
3122 let tensor_covers = self
3129 .frozen_glm_tensor
3130 .as_ref()
3131 .is_some_and(|t| t.contains(psi));
3132 let current_w = if tensor_covers {
3133 match warm_beta {
3134 Some(beta) => self.frozen_glm_trial_weights(beta)?,
3135 None => None,
3136 }
3137 } else {
3138 None
3139 };
3140 if let (Some(tensor), Some(current_w)) =
3141 (self.frozen_glm_tensor.as_ref(), current_w.as_ref())
3142 {
3143 const FROZEN_GLM_WEIGHT_DRIFT_RTOL: f64 = 1e-3;
3144 if tensor.weight_drift_within(current_w.view(), FROZEN_GLM_WEIGHT_DRIFT_RTOL) {
3145 staged_gram = Some(tensor.gram_at(psi));
3146 log::debug!(
3147 "[STAGE] {} trial at psi={psi:.6}: serving frozen-W GLM \
3148 first-Fisher-step XᵀWX n-free (weight drift within tol)",
3149 kind.label(),
3150 );
3151 }
3152 if allow_gradient
3153 && tensor.contains_for_gradient(psi)
3154 && let Some((dgram_dpsi, drhs_dpsi)) =
3155 tensor.gradient_pair_if_sound(psi, current_w.view())
3156 {
3157 staged_deriv = Some((dgram_dpsi, drhs_dpsi));
3158 log::debug!(
3159 "[STAGE] {} trial at psi={psi:.6}: serving frozen-W GLM \
3160 ψ-gradient (∂G/∂ψ, ∂b/∂ψ) n-free (gradient weight drift within \
3161 tight tol); B_j stays exact",
3162 kind.label(),
3163 );
3164 }
3165 }
3166 }
3167 self.evaluator.stage_glm_first_step_gram(staged_gram);
3168 self.evaluator.stage_glm_psi_gram_deriv(staged_deriv);
3169 Ok(())
3170 }
3171
3172 fn eval_full(
3174 &mut self,
3175 theta: &Array1<f64>,
3176 order: gam_solve::rho_optimizer::OuterEvalOrder,
3177 analytic_outer_hessian_available: bool,
3178 ) -> Result<(f64, Array1<f64>, gam_problem::HessianValue), EstimationError> {
3179 use gam_solve::rho_optimizer::OuterEvalOrder;
3180 let allow_second_order = matches!(order, OuterEvalOrder::ValueGradientHessian)
3181 && analytic_outer_hessian_available;
3182 if let Some(eval) = self.cache.memoized_eval(theta) {
3183 let cached_satisfies_order = !allow_second_order || eval.2.is_analytic();
3184 if cached_satisfies_order {
3185 return Ok(eval);
3186 }
3187 }
3188 let kind = self.kind;
3189 let nfree_fast_path_revision = self.evaluator.nfree_fast_path_revision();
3225 let skip_design_realization = !allow_second_order && theta.len() == self.rho_dim + 1 && {
3226 let psi = theta[self.rho_dim];
3227 self.evaluator.psi_gram_tensor_covers(psi)
3228 && self.evaluator.psi_gram_tensor_covers_gradient(psi)
3235 && self.evaluator.psi_gram_tensor_covers_skip(psi)
3252 && self.evaluator.supports_nfree_penalty_rekey()
3257 && nfree_fast_path_revision.is_some()
3258 };
3259 if skip_design_realization {
3271 log::debug!(
3272 "[STAGE] {} eval_full at psi={:.6}: skipping n×k design re-realization \
3273 + reconditioning — criterion/gradient/inner-solve served n-free from \
3274 the certified ψ-gram tensor (GaussianFixedCache + k-space ψ-derivatives)",
3275 kind.label(),
3276 theta[self.rho_dim],
3277 );
3278 } else {
3279 self.cache
3280 .ensure_theta(theta)?;
3281 }
3282 let warm_beta = self.evaluator.current_beta();
3283 self.ensure_frozen_glm_tensor(theta, warm_beta.as_ref())?;
3284 self.stage_frozen_glm_trial_statistics(theta, warm_beta.as_ref(), !allow_second_order)?;
3292 let hyper_dirs = if skip_design_realization {
3299 self.cache.nfree_tensor_gradient_hyper_dirs(theta)?
3300 } else {
3301 self.cache.hyper_dirs_for_current_design(self.data, kind)?
3302 };
3303
3304 let design_revision = if skip_design_realization {
3305 nfree_fast_path_revision
3306 } else {
3307 Some(self.cache.design_revision())
3308 };
3309 if self.evaluator.supports_nfree_penalty_rekey() {
3323 match self.cache.canonical_penalties_at(theta) {
3324 Ok(penalty) => self.evaluator.stage_fast_path_penalty(Some(penalty)),
3325 Err(e) => {
3326 log::warn!(
3327 "[STAGE] {} eval_full at psi={:.6}: exact n-free S(ψ) rebuild failed \
3328 ({e}); clearing stage (eval falls to slow path)",
3329 kind.label(),
3330 theta[self.rho_dim],
3331 );
3332 self.evaluator.stage_fast_path_penalty(None);
3333 }
3334 }
3335 }
3336 let eval = evaluate_joint_reml_outer_eval_at_theta(
3343 &mut self.evaluator,
3344 self.cache.design(),
3345 theta,
3346 self.rho_dim,
3347 hyper_dirs,
3348 warm_beta.as_ref().map(|b: &Array1<f64>| b.view()),
3349 if allow_second_order {
3350 order
3351 } else {
3352 OuterEvalOrder::ValueAndGradient
3353 },
3354 design_revision,
3355 );
3356 if let Ok(ref value) = eval {
3357 self.cache.store_eval_at(theta, value.clone());
3358 }
3359 eval
3360 }
3361
3362 fn eval_efs(&mut self, theta: &Array1<f64>) -> Result<gam_problem::EfsEval, EstimationError> {
3363 self.cache
3364 .ensure_theta(theta)?;
3365 let kind = self.kind;
3366 let hyper_dirs = try_build_spatial_log_kappa_hyper_dirs(
3367 self.data,
3368 self.cache.spec(),
3369 self.cache.design(),
3370 &self.cache.spatial_terms,
3371 )?
3372 .ok_or_else(|| {
3373 EstimationError::InvalidInput(format!(
3374 "failed to build {} hyper_dirs for exact-joint EFS",
3375 kind.adjective(),
3376 ))
3377 })?;
3378 let design_revision = Some(self.cache.design_revision());
3379 let warm_beta = self.evaluator.current_beta();
3380 evaluate_joint_reml_efs_at_theta(
3381 &mut self.evaluator,
3382 self.cache.design(),
3383 theta,
3384 self.rho_dim,
3385 hyper_dirs,
3386 warm_beta.as_ref().map(|b: &Array1<f64>| b.view()),
3387 design_revision,
3388 )
3389 }
3390
3391 fn eval_cost(&mut self, theta: &Array1<f64>) -> Result<f64, EstimationError> {
3397 if let Some(cost) = self.cache.memoized_cost(theta) {
3398 return Ok(cost);
3399 }
3400 let probe_start = std::time::Instant::now();
3415 let psi_distance = self
3416 .cache
3417 .current_theta
3418 .as_ref()
3419 .filter(|reference| reference.len() == theta.len())
3420 .map(|reference| {
3421 reference
3422 .iter()
3423 .zip(theta.iter())
3424 .map(|(a, b)| (a - b) * (a - b))
3425 .sum::<f64>()
3426 .sqrt()
3427 })
3428 .unwrap_or(f64::NAN);
3429 let nfree_fast_path_revision = self.evaluator.nfree_fast_path_revision();
3443 let skip_value_realization = theta.len() == self.rho_dim + 1 && {
3444 let psi = theta[self.rho_dim];
3445 self.evaluator.psi_gram_tensor_covers(psi)
3446 && self.evaluator.supports_nfree_penalty_rekey()
3480 && nfree_fast_path_revision.is_some()
3481 };
3482 if theta.len() == self.rho_dim + 1
3483 && self.evaluator.has_psi_gram_tensor()
3484 && !self.evaluator.psi_gram_tensor_covers(theta[self.rho_dim])
3485 {
3486 self.cache.store_cost_at(theta, f64::INFINITY);
3487 return Ok(f64::INFINITY);
3488 }
3489 if !skip_value_realization && let Err(error) = self.cache.ensure_theta(theta) {
3493 self.value_realization_failures += 1;
3494 let (theta_norm, log_kappa_norm) = kphase_log_norms(theta, self.rho_dim);
3495 if is_recoverable_trial_point_error(&error) {
3496 log::debug!(
3497 "[STAGE] {} value-probe: design realization makes this trial infeasible at theta_norm={:.4e} log_kappa_norm={:.4e} ({error}); retreating",
3498 self.kind.label(), theta_norm, log_kappa_norm,
3499 );
3500 } else {
3501 log::warn!(
3502 "[STAGE] {} value-probe: design realization FAILED fatally at theta_norm={:.4e} log_kappa_norm={:.4e} ({error}); propagating",
3503 self.kind.label(), theta_norm, log_kappa_norm,
3504 );
3505 }
3506 return classify_spatial_value_probe_failure(error);
3507 }
3508 if self.evaluator.supports_nfree_penalty_rekey() {
3514 match self.cache.canonical_penalties_at(theta) {
3515 Ok(penalty) => self.evaluator.stage_fast_path_penalty(Some(penalty)),
3516 Err(_) => self.evaluator.stage_fast_path_penalty(None),
3517 }
3518 }
3519 let warm_beta = self.evaluator.current_beta();
3520 if let Err(err) = self.ensure_frozen_glm_tensor(theta, warm_beta.as_ref()) {
3521 log::warn!(
3522 "[STAGE] {} value-probe at psi={:.6}: frozen-W GLM tensor setup failed ({err}); \
3523 falling back to exact streamed Gram",
3524 self.kind.label(),
3525 if theta.len() > self.rho_dim {
3526 theta[self.rho_dim]
3527 } else {
3528 f64::NAN
3529 },
3530 );
3531 self.evaluator.stage_glm_first_step_gram(None);
3532 self.evaluator.stage_glm_psi_gram_deriv(None);
3533 } else if let Err(err) =
3534 self.stage_frozen_glm_trial_statistics(theta, warm_beta.as_ref(), false)
3535 {
3536 log::warn!(
3537 "[STAGE] {} value-probe at psi={:.6}: frozen-W GLM staging failed ({err}); \
3538 falling back to exact streamed Gram",
3539 self.kind.label(),
3540 if theta.len() > self.rho_dim {
3541 theta[self.rho_dim]
3542 } else {
3543 f64::NAN
3544 },
3545 );
3546 self.evaluator.stage_glm_first_step_gram(None);
3547 self.evaluator.stage_glm_psi_gram_deriv(None);
3548 }
3549 let design_revision = if skip_value_realization {
3550 nfree_fast_path_revision
3551 } else {
3552 Some(self.cache.design_revision())
3553 };
3554 let cost_label = self.kind.label();
3555 let result = {
3556 let design = self.cache.design();
3557 self.evaluator.evaluate_cost_only(
3558 &design.design,
3559 &design.penalties,
3560 &design.nullspace_dims,
3561 design.linear_constraints.clone(),
3562 theta,
3563 self.rho_dim,
3564 warm_beta.as_ref().map(|b: &Array1<f64>| b.view()),
3565 cost_label,
3566 design_revision,
3567 )
3568 };
3569 match result {
3570 Ok(cost) => {
3571 log::debug!(
3572 "[STAGE] {cost_label} value-probe (order=Value): elapsed={:.3}s \
3573 cost={cost:.6e} trial_theta_distance={psi_distance:.3e}",
3574 probe_start.elapsed().as_secs_f64(),
3575 );
3576 self.cache.store_cost_at(theta, cost);
3577 Ok(cost)
3578 }
3579 Err(error) => {
3582 self.value_evaluation_failures += 1;
3583 let (theta_norm, log_kappa_norm) = kphase_log_norms(theta, self.rho_dim);
3584 if is_recoverable_trial_point_error(&error) {
3585 log::debug!(
3586 "[STAGE] {cost_label} value-probe: cost evaluator makes this trial infeasible at theta_norm={theta_norm:.4e} log_kappa_norm={log_kappa_norm:.4e} ({error}); retreating",
3587 );
3588 } else {
3589 log::warn!(
3590 "[STAGE] {cost_label} value-probe: cost evaluation FAILED fatally at theta_norm={theta_norm:.4e} log_kappa_norm={log_kappa_norm:.4e} ({error}); propagating",
3591 );
3592 }
3593 classify_spatial_value_probe_failure(error)
3594 }
3595 }
3596 }
3597
3598 fn reset(&mut self) {
3599 self.cache.current_theta = None;
3600 self.cache.last_eval_theta = None;
3601 self.cache.last_cost = None;
3602 self.cache.last_eval = None;
3603 }
3604}
3605
3606fn kphase_psi_display(theta: &Array1<f64>, rho_dim: usize) -> String {
3639 let mut out = String::from("[");
3640 for (offset, value) in theta.iter().skip(rho_dim).enumerate() {
3641 if offset > 0 {
3642 out.push(',');
3643 }
3644 out.push_str(&format!("{value:+.4e}"));
3645 }
3646 out.push(']');
3647 out
3648}
3649
3650fn kphase_log_norms(theta: &Array1<f64>, rho_dim: usize) -> (f64, f64) {
3651 let theta_norm = theta.iter().map(|v| v * v).sum::<f64>().sqrt();
3652 let log_kappa_norm = theta
3653 .iter()
3654 .skip(rho_dim)
3655 .map(|v| v * v)
3656 .sum::<f64>()
3657 .sqrt();
3658 (theta_norm, log_kappa_norm)
3659}
3660
3661fn run_exact_joint_spatial_optimization(
3662 kind: SpatialHyperKind,
3663 data: ArrayView2<'_, f64>,
3664 y: ArrayView1<'_, f64>,
3665 weights: ArrayView1<'_, f64>,
3666 offset: ArrayView1<'_, f64>,
3667 resolvedspec: &TermCollectionSpec,
3668 baseline_design: &TermCollectionDesign,
3669 family: LikelihoodSpec,
3670 options: &FitOptions,
3671 spatial_terms: &[usize],
3672 dims_per_term: &[usize],
3673 theta0: &Array1<f64>,
3674 lower: &Array1<f64>,
3675 upper: &Array1<f64>,
3676 rho_dim: usize,
3677 kappa_options: &SpatialLengthScaleOptimizationOptions,
3678) -> Result<(Array1<f64>, f64, f64, SpatialLengthScaleOptimizationTiming), EstimationError> {
3679 let label = kind.label();
3680 let effective_offset = baseline_design
3681 .compose_offset(offset, "spatial joint fit")
3682 .map_err(EstimationError::BasisError)?;
3683 let offset = effective_offset.view();
3684 let external_opts = external_opts_for_design(&family, baseline_design, options);
3685 let joint_conditioned_y = gam_solve::estimate::gaussian_identity_outer_response_conditioning(
3719 &baseline_design.design,
3720 &baseline_design.penalties,
3721 &external_opts,
3722 y,
3723 weights,
3724 offset,
3725 )?;
3726 if joint_conditioned_y.is_some() {
3727 log::info!(
3728 "[{label}] outer response conditioned for the joint [rho, psi] search (#2671): the \
3729 criterion is now formed in the same coordinates as the scalar-rho route it is \
3730 graded against"
3731 );
3732 }
3733 let y = joint_conditioned_y
3734 .as_ref()
3735 .map_or(y, |conditioned| conditioned.view());
3736 assert!(
3738 lower.len() == theta0.len() && upper.len() == theta0.len(),
3739 "spatial hyperparameter bounds must match theta length: lower_len={}, upper_len={}, theta_len={}",
3740 lower.len(),
3741 upper.len(),
3742 theta0.len()
3743 );
3744 assert!(
3745 baseline_design.smooth.terms.len() >= spatial_terms.len(),
3746 "baseline design must have at least one smooth term per spatial term: baseline_terms={}, spatial_terms={}",
3747 baseline_design.smooth.terms.len(),
3748 spatial_terms.len()
3749 );
3750 use gam_problem::{DeclaredHessianForm, Derivative, OuterEval};
3751 use gam_solve::rho_optimizer::OuterEvalOrder;
3752
3753 let theta_dim = theta0.len();
3754 let coord_dim = theta_dim - rho_dim;
3757 let analytic_outer_hessian_available =
3762 exact_joint_spatial_outer_hessian_available(&family, baseline_design);
3763 if !analytic_outer_hessian_available {
3764 log::info!(
3765 "[{label}] analytic outer Hessian unavailable for family/design; routing without second-order geometry (coord_dim={coord_dim})"
3766 );
3767 }
3768 let mut suppress_outer_hessian_for_nfree = false;
3783
3784 log::trace!(
3785 "[{}] starting analytic optimization: rho_dim={}, coord_dim={}, dims_per_term={:?}",
3786 label,
3787 rho_dim,
3788 coord_dim,
3789 dims_per_term,
3790 );
3791
3792 let mut ctx = SpatialJointContext {
3793 data,
3794 rho_dim,
3795 kind,
3796 value_realization_failures: 0,
3797 value_evaluation_failures: 0,
3798 nfree_polish_boundary: None,
3799 cache: SingleBlockExactJointDesignCache::new_with_policy(
3800 data,
3801 resolvedspec.clone(),
3802 baseline_design.clone(),
3803 spatial_terms.to_vec(),
3804 rho_dim,
3805 dims_per_term.to_vec(),
3806 &options.resource_policy,
3807 )
3808 .map_err(EstimationError::InvalidInput)?,
3809 evaluator: gam_solve::estimate::ExternalJointHyperEvaluator::new(
3810 y,
3811 weights,
3812 &baseline_design.design,
3813 offset,
3814 &baseline_design.penalties,
3815 &external_opts,
3816 label,
3817 )?,
3818 frozen_glm_inputs: if coord_dim == 1 && frozen_glm_tensor_eligible_family(&family) {
3819 Some(SpatialFrozenGlmInputs {
3820 y: y.to_owned(),
3821 weights: weights.to_owned(),
3822 offset: offset.to_owned(),
3823 family: family.clone(),
3824 })
3825 } else {
3826 None
3827 },
3828 frozen_glm_psi_bounds: if coord_dim == 1 && frozen_glm_tensor_eligible_family(&family) {
3829 Some((lower[rho_dim], upper[rho_dim]))
3830 } else {
3831 None
3832 },
3833 frozen_glm_tensor: None,
3834 frozen_glm_tensor_attempted: false,
3835 frozen_glm_weight_memo: None,
3836 };
3837
3838 let mut psi_rank_stable_floor: Option<f64> = None;
3861 let mut psi_rank_stable_ceiling: Option<f64> = None;
3870 let nfree_penalty_capable =
3871 coord_dim == 1 && family.is_gaussian_identity() && ctx.cache.supports_nfree_penalty_rekey();
3872 if nfree_penalty_capable {
3873 let psi_lo = lower[rho_dim];
3874 let psi_hi = upper[rho_dim];
3875 let z = Array1::from_iter(y.iter().zip(offset.iter()).map(|(yi, oi)| yi - oi));
3876 let theta_probe_base = theta0.clone();
3877 let SpatialJointContext {
3880 cache, evaluator, ..
3881 } = &mut ctx;
3882 let attached = evaluator.build_and_set_psi_gram_tensor(
3883 |psi| {
3884 let mut theta_probe = theta_probe_base.clone();
3885 theta_probe[rho_dim] = psi;
3886 cache.ensure_theta(&theta_probe).map_err(|e| e.to_string())?;
3887 Ok(cache.design().design.clone())
3888 },
3889 weights,
3890 z.view(),
3891 psi_lo,
3892 psi_hi,
3893 );
3894 if attached {
3895 log::info!(
3896 "[{label}] certified ψ-gram tensor over [{psi_lo:.3}, {psi_hi:.3}]: \
3897 in-window trials assemble Gaussian sufficient statistics n-free"
3898 );
3899 let psi_anchor = theta0[rho_dim];
3904 let psi_projector_bar = evaluator.psi_gram_projector_error_bar(psi_anchor);
3910 let psi_rank_stable_floor_raw = evaluator.psi_gram_rank_stable_floor(psi_anchor);
3913 psi_rank_stable_floor = psi_rank_stable_floor_raw
3914 .filter(|&f| f.is_finite() && f > psi_lo && f < psi_anchor);
3915 log::info!(
3916 "[KAPPA-PHASE-FLOOR] n_rows={} psi_lo={psi_lo:.6} psi_anchor={psi_anchor:.6} \
3917 rank_stable_floor={psi_rank_stable_floor_raw:?} lifted={} \
3918 projector_error_bar={psi_projector_bar:?}",
3919 data.nrows(),
3920 psi_rank_stable_floor.is_some(),
3921 );
3922 if let Some(floor) = psi_rank_stable_floor {
3923 log::info!(
3924 "[{label}] rank-stable κ-floor ψ_floor={floor:.6} > window floor \
3925 ψ_lo={psi_lo:.6}: lifting the optimizer lower bound to keep every \
3926 in-window trial on the n-free design-realization skip (#1033). The \
3927 conditioned Gram is rank-deficient below ψ_floor (longest-length-scale \
3928 radial mode collapses into the nullspace), where the skip is soundly \
3929 refused. The SEARCH is n-free — O(iters·k³) off the k-space tensor, \
3930 zero row access — but the EDGE IS NOT AN n-INVARIANT CONSTANT of the \
3931 design (#2408): the tensor is built from n rows, so its Gram is an \
3932 O(1/n) relative perturbation of the continuum Gram, which moves the \
3933 rank margin additively and displaces this root by \
3934 sup|δ margin| / inf|d margin/dψ|. A steep cliff pins it to machine \
3935 precision; a grazing crossing does not. Treat it as a clamp carrying \
3936 that transport bound, not as the n-independent answer."
3937 );
3938 }
3939 let psi_rank_stable_ceiling_raw = evaluator.psi_gram_rank_stable_ceiling(psi_anchor);
3948 psi_rank_stable_ceiling = psi_rank_stable_ceiling_raw
3949 .filter(|&c| c.is_finite() && c < psi_hi && c > psi_anchor);
3950 log::info!(
3951 "[KAPPA-PHASE-CEIL] n_rows={} psi_hi={psi_hi:.6} psi_anchor={psi_anchor:.6} \
3952 rank_stable_ceiling={psi_rank_stable_ceiling_raw:?} clamped={} \
3953 projector_error_bar={psi_projector_bar:?}",
3954 data.nrows(),
3955 psi_rank_stable_ceiling.is_some(),
3956 );
3957 if let Some(ceiling) = psi_rank_stable_ceiling {
3958 log::info!(
3959 "[{label}] rank-stable κ-ceiling ψ_ceil={ceiling:.6} < window ceiling \
3960 ψ_hi={psi_hi:.6}: clamping the optimizer upper bound to keep every \
3961 in-window trial on the n-free design-realization skip (#1033). The \
3962 conditioned Gram is rank-deficient above ψ_ceil (longest-frequency \
3963 radial mode goes collinear), where the skip is soundly refused; a \
3964 line-search overshoot there trips the O(n) reset_surface lane (and the \
3965 deficient pinning ψ it records resets the next in-band trial too)."
3966 );
3967 }
3968 if let Some(bar) = psi_projector_bar
3977 && bar > gam_solve::psi_gram_tensor::PSI_GRAM_SKIP_PROJ_ATOL
3978 {
3979 log::warn!(
3980 "[{label}] ψ-gram range projector at the anchor ψ={psi_anchor:.6} is \
3981 UNRESOLVED: Davis–Kahan bar {bar:.3e} exceeds the {:.3e} subspace \
3982 tolerance the design-revision skip gates on (#2448). The conditioned \
3983 Gram has no kept/dropped eigen-gap wide enough to decide subspace \
3984 identity at double precision here — its spectrum decays smoothly \
3985 through the rank cutoff instead of cliffing — so the skip witness \
3986 soundly refuses every trial and the n-free fast path will not fire \
3987 at all. Results are unaffected (the exact O(n) path runs); the cost \
3988 is the fast path. The lever is the geometry (basis size / centers) \
3989 or the rank cutoff, not this clamp.",
3990 gam_solve::psi_gram_tensor::PSI_GRAM_SKIP_PROJ_ATOL
3991 );
3992 }
3993 let gradient_covers_full_window = evaluator.psi_gram_tensor_covers_gradient(psi_lo)
3994 && evaluator.psi_gram_tensor_covers_gradient(psi_hi);
3995 if gradient_covers_full_window {
3996 log::info!(
3997 "[{label}] certified ψ-gram tensor gradient lane covers the full \
3998 optimizer window [{psi_lo:.3}, {psi_hi:.3}]"
3999 );
4000 } else {
4001 log::info!(
4002 "[{label}] ψ-gram tensor value lane certified, but the gradient lane \
4003 does not cover the full optimizer window [{psi_lo:.3}, {psi_hi:.3}]; \
4004 keeping exact streamed kappa routing"
4005 );
4006 }
4007 evaluator.set_supports_nfree_penalty_rekey(true);
4027 log::info!(
4028 "[{label}] exact n-free ψ-penalty re-key enabled over [{psi_lo:.3}, \
4029 {psi_hi:.3}]: in-window fast-path trials rebuild S(ψ) n-free from frozen \
4030 geometry (no reset_surface)"
4031 );
4032 } else {
4033 log::info!(
4034 "[{label}] ψ-gram tensor did not certify over [{psi_lo:.3}, {psi_hi:.3}]; \
4035 keeping the exact per-trial path"
4036 );
4037 }
4038 if attached
4059 && evaluator.psi_gram_tensor_covers_gradient(psi_lo)
4060 && evaluator.psi_gram_tensor_covers_gradient(psi_hi)
4061 && evaluator.supports_nfree_penalty_rekey()
4062 && cache.supports_nfree_gradient_only_routing()
4063 {
4064 suppress_outer_hessian_for_nfree = true;
4065 log::info!(
4066 "[{label}] n-free Gaussian ψ-lane armed; routing the SEARCH gradient-only \
4067 (BFGS, fixed-point lane off) so no in-window κ-trial realizes the O(n) \
4068 second-order slab — n-independent outer loop (#1033). The terminal \
4069 certificate keeps its one exact curvature evaluation (gam#2760)."
4070 );
4071 }
4072 } else if coord_dim == 1 && family.is_gaussian_identity() {
4073 log::info!(
4074 "[{label}] exact n-free ψ-penalty re-key unavailable; skipping ψ-gram tensor \
4075 attachment so value, gradient, and Hessian remain on the same exact streamed \
4076 objective"
4077 );
4078 }
4079
4080 let kphase_prime_order = OuterEvalOrder::ValueAndGradient;
4083 let kphase_prime_start = std::time::Instant::now();
4084 let seed_value = ctx
4093 .eval_full(theta0, kphase_prime_order, analytic_outer_hessian_available)?
4094 .0;
4095 log::info!(
4096 "[KAPPA-PHASE-PRIME] n_rows={} order={:?} seed_value={seed_value:.12e} elapsed_s={:.4} slow_path_resets_total={} design_revision={}",
4097 data.nrows(),
4098 kphase_prime_order,
4099 kphase_prime_start.elapsed().as_secs_f64(),
4100 ctx.evaluator.slow_path_reset_count(),
4101 ctx.cache.design_revision(),
4102 );
4103
4104 let kphase_cost_calls = std::cell::Cell::new(0usize);
4105 let kphase_eval_calls = std::cell::Cell::new(0usize);
4106 let kphase_efs_calls = std::cell::Cell::new(0usize);
4107 let kphase_cost_total_s = std::cell::Cell::new(0.0);
4108 let kphase_eval_total_s = std::cell::Cell::new(0.0);
4109 let kphase_efs_total_s = std::cell::Cell::new(0.0);
4110 let kphase_nfree_miss_shape = std::cell::Cell::new(0u64);
4111 let kphase_nfree_miss_value = std::cell::Cell::new(0u64);
4112 let kphase_nfree_miss_gradient = std::cell::Cell::new(0u64);
4113 let kphase_nfree_miss_penalty = std::cell::Cell::new(0u64);
4114 let kphase_nfree_miss_revision = std::cell::Cell::new(0u64);
4115 let kphase_nfree_miss_second_order = std::cell::Cell::new(0u64);
4116 let kphase_nfree_miss_other = std::cell::Cell::new(0u64);
4117 let kphase_optim_start = std::time::Instant::now();
4118 let kphase_log_kappa_dim = coord_dim;
4119 let kphase_slow_resets_start = ctx.evaluator.slow_path_reset_count();
4120 let kphase_design_revision_start = ctx.cache.design_revision();
4121 let kphase_nfree_skip_touches_start = gam_solve::pirls::nfree_skip_row_element_touches();
4125
4126 let lower_effective: std::borrow::Cow<'_, Array1<f64>> = match psi_rank_stable_floor {
4133 Some(floor) if coord_dim == 1 && floor > lower[rho_dim] => {
4134 let mut lifted = lower.clone();
4135 lifted[rho_dim] = floor;
4136 std::borrow::Cow::Owned(lifted)
4137 }
4138 _ => std::borrow::Cow::Borrowed(lower),
4139 };
4140 let lower = lower_effective.as_ref();
4141
4142 let upper_effective: std::borrow::Cow<'_, Array1<f64>> = match psi_rank_stable_ceiling {
4150 Some(ceiling) if coord_dim == 1 && ceiling < upper[rho_dim] => {
4151 let mut clamped = upper.clone();
4152 clamped[rho_dim] = ceiling;
4153 std::borrow::Cow::Owned(clamped)
4154 }
4155 _ => std::borrow::Cow::Borrowed(upper),
4156 };
4157 let upper = upper_effective.as_ref();
4158
4159 let problem = exact_joint_multistart_outer_problem(
4160 theta0,
4161 lower,
4162 upper,
4163 rho_dim,
4164 coord_dim,
4165 theta_dim,
4166 Derivative::Analytic,
4167 if analytic_outer_hessian_available && !suppress_outer_hessian_for_nfree {
4168 DeclaredHessianForm::Either
4239 } else {
4240 DeclaredHessianForm::Unavailable
4241 },
4242 suppress_outer_hessian_for_nfree,
4253 seed_risk_profile_for_likelihood_family(&family),
4254 kappa_options.rel_tol.max(1e-6),
4255 kappa_options.max_outer_iter.max(1),
4256 Some(5.0),
4259 Some(kappa_options.log_step.clamp(0.25, 1.0)),
4261 None,
4262 Some((data.nrows(), baseline_design.design.ncols())),
4267 !constant_curvature_term_indices(resolvedspec).is_empty(),
4271 kind == SpatialHyperKind::Isotropic
4276 && constant_curvature_term_indices(resolvedspec).is_empty()
4277 && spatial_terms.iter().any(|&term_idx| {
4278 matches!(
4279 resolvedspec
4280 .smooth_terms
4281 .get(term_idx)
4282 .map(|term| &term.basis),
4283 Some(SmoothBasisSpec::Matern { .. })
4284 )
4285 }),
4286 )?;
4287
4288 let eval_outer = |ctx: &mut &mut SpatialJointContext<'_>,
4289 theta: &Array1<f64>,
4290 order: OuterEvalOrder|
4291 -> Result<OuterEval, EstimationError> {
4292 let t0 = std::time::Instant::now();
4293 let allow_second_order_for_call = matches!(order, OuterEvalOrder::ValueGradientHessian)
4294 && analytic_outer_hessian_available;
4295 let gate = ctx.nfree_skip_gate_status(theta, allow_second_order_for_call, true);
4296 let resets_before = ctx.evaluator.slow_path_reset_count();
4297 let raw = ctx.eval_full(theta, order, analytic_outer_hessian_available);
4298 let reset_delta = ctx
4299 .evaluator
4300 .slow_path_reset_count()
4301 .saturating_sub(resets_before);
4302 if reset_delta > 0 {
4303 if !gate.shape {
4304 kphase_nfree_miss_shape.set(kphase_nfree_miss_shape.get() + reset_delta);
4305 }
4306 if gate.shape && !gate.value {
4307 kphase_nfree_miss_value.set(kphase_nfree_miss_value.get() + reset_delta);
4308 }
4309 if gate.shape && gate.value && !gate.gradient {
4310 kphase_nfree_miss_gradient.set(kphase_nfree_miss_gradient.get() + reset_delta);
4311 }
4312 if gate.shape && gate.value && gate.gradient && !gate.penalty {
4313 kphase_nfree_miss_penalty.set(kphase_nfree_miss_penalty.get() + reset_delta);
4314 }
4315 if gate.shape && gate.value && gate.gradient && gate.penalty && !gate.revision {
4316 kphase_nfree_miss_revision.set(kphase_nfree_miss_revision.get() + reset_delta);
4317 }
4318 if gate.shape
4319 && gate.value
4320 && gate.gradient
4321 && gate.penalty
4322 && gate.revision
4323 && gate.second_order
4324 {
4325 kphase_nfree_miss_second_order
4326 .set(kphase_nfree_miss_second_order.get() + reset_delta);
4327 }
4328 if gate.would_skip(true) {
4329 kphase_nfree_miss_other.set(kphase_nfree_miss_other.get() + reset_delta);
4330 }
4331 }
4332 let elapsed_s = t0.elapsed().as_secs_f64();
4333 kphase_eval_calls.set(kphase_eval_calls.get() + 1);
4334 kphase_eval_total_s.set(kphase_eval_total_s.get() + elapsed_s);
4335 let (theta_norm, log_kappa_norm) = kphase_log_norms(theta, rho_dim);
4336 log::info!(
4337 "[KAPPA-PHASE] phase=eval_outer call={} order={:?} design_revision={:?} theta_norm={:.4e} log_kappa_norm={:.4e} psi={} elapsed_s={:.4}",
4338 kphase_eval_calls.get(),
4339 order,
4340 Some(ctx.cache.design_revision()),
4341 theta_norm,
4342 log_kappa_norm,
4343 kphase_psi_display(theta, rho_dim),
4344 elapsed_s,
4345 );
4346 match raw {
4347 Ok((cost, grad, hess)) => Ok(OuterEval {
4348 cost,
4349 gradient: grad,
4350 hessian: hess,
4351 inner_beta_hint: None,
4352 }),
4353 Err(err) if is_recoverable_trial_point_error(&err) => {
4361 log::debug!(
4362 "[{label}] trial point infeasible (kernel design \
4363 not constructible at theta={theta:?}): {err}; retreating",
4364 );
4365 Ok(OuterEval::infeasible(theta_dim))
4366 }
4367 Err(err) => Err(err),
4368 }
4369 };
4370
4371 let obj = problem.build_objective_with_eval_order(
4372 &mut ctx,
4373 |ctx: &mut &mut SpatialJointContext<'_>, theta: &Array1<f64>| {
4374 let t0 = std::time::Instant::now();
4375 let gate = ctx.nfree_skip_gate_status(theta, false, false);
4376 let resets_before = ctx.evaluator.slow_path_reset_count();
4377 let cost = ctx.eval_cost(theta);
4378 let reset_delta = ctx
4379 .evaluator
4380 .slow_path_reset_count()
4381 .saturating_sub(resets_before);
4382 if reset_delta > 0 {
4383 if !gate.shape {
4384 kphase_nfree_miss_shape.set(kphase_nfree_miss_shape.get() + reset_delta);
4385 }
4386 if gate.shape && !gate.value {
4387 kphase_nfree_miss_value.set(kphase_nfree_miss_value.get() + reset_delta);
4388 }
4389 if gate.shape && gate.value && !gate.penalty {
4390 kphase_nfree_miss_penalty.set(kphase_nfree_miss_penalty.get() + reset_delta);
4391 }
4392 if gate.shape && gate.value && gate.penalty && !gate.revision {
4393 kphase_nfree_miss_revision.set(kphase_nfree_miss_revision.get() + reset_delta);
4394 }
4395 if gate.would_skip(false) {
4396 kphase_nfree_miss_other.set(kphase_nfree_miss_other.get() + reset_delta);
4397 }
4398 }
4399 let elapsed_s = t0.elapsed().as_secs_f64();
4400 kphase_cost_calls.set(kphase_cost_calls.get() + 1);
4401 kphase_cost_total_s.set(kphase_cost_total_s.get() + elapsed_s);
4402 let (theta_norm, log_kappa_norm) = kphase_log_norms(theta, rho_dim);
4403 log::info!(
4404 "[KAPPA-PHASE] phase=cost call={} design_revision={:?} theta_norm={:.4e} log_kappa_norm={:.4e} elapsed_s={:.4}",
4405 kphase_cost_calls.get(),
4406 Some(ctx.cache.design_revision()),
4407 theta_norm,
4408 log_kappa_norm,
4409 elapsed_s,
4410 );
4411 cost
4412 },
4413 |ctx: &mut &mut SpatialJointContext<'_>, theta: &Array1<f64>| {
4414 eval_outer(
4415 ctx,
4416 theta,
4417 OuterEvalOrder::ValueAndGradient,
4421 )
4422 },
4423 |ctx: &mut &mut SpatialJointContext<'_>, theta: &Array1<f64>, order: OuterEvalOrder| {
4424 eval_outer(ctx, theta, order)
4425 },
4426 Some(|ctx: &mut &mut SpatialJointContext<'_>| {
4427 ctx.reset();
4428 }),
4429 Some(|ctx: &mut &mut SpatialJointContext<'_>, theta: &Array1<f64>| {
4430 let t0 = std::time::Instant::now();
4431 let eval = ctx.eval_efs(theta);
4432 let elapsed_s = t0.elapsed().as_secs_f64();
4433 kphase_efs_calls.set(kphase_efs_calls.get() + 1);
4434 kphase_efs_total_s.set(kphase_efs_total_s.get() + elapsed_s);
4435 let (theta_norm, log_kappa_norm) = kphase_log_norms(theta, rho_dim);
4436 log::info!(
4437 "[KAPPA-PHASE] phase=efs call={} design_revision={:?} theta_norm={:.4e} log_kappa_norm={:.4e} elapsed_s={:.4}",
4438 kphase_efs_calls.get(),
4439 Some(ctx.cache.design_revision()),
4440 theta_norm,
4441 log_kappa_norm,
4442 elapsed_s,
4443 );
4444 eval
4445 }),
4446 );
4447
4448 let mut obj = obj
4459 .with_criterion_invariance(
4460 |ctx: &mut &mut SpatialJointContext<'_>, rho: &Array1<f64>| {
4461 ctx.evaluator.criterion_invariant_directions(rho)
4462 },
4463 )
4464 .with_exact_polish(|ctx: &mut &mut SpatialJointContext<'_>| {
4486 if !ctx.evaluator.retire_psi_gram_tensor() {
4487 return false;
4488 }
4489 ctx.cache.forget_eval_memo();
4492 ctx.nfree_polish_boundary = Some((
4495 ctx.evaluator.slow_path_reset_count(),
4496 gam_solve::pirls::nfree_skip_row_element_touches(),
4497 ));
4498 log::info!(
4499 "[KAPPA-PHASE-POLISH] the certified n-free psi-Gram surrogate is retired at \
4500 the search checkpoint; the optimizer continues and certifies on the exact \
4501 streamed criterion (gam#2760)"
4502 );
4503 true
4504 });
4505
4506 let run_label = match kind {
4507 SpatialHyperKind::Anisotropic => "aniso-psi joint REML",
4508 SpatialHyperKind::Isotropic => "iso-kappa joint REML",
4509 };
4510 let result = problem.run(&mut obj, run_label)?;
4511 if !result.converged() {
4512 crate::bail_invalid_estim!(
4513 "{} did not converge after {} iterations (final_objective={:.6e}, final_grad_norm={})",
4514 run_label,
4515 result.iterations,
4516 result.final_value,
4517 result.final_grad_norm_report(),
4518 );
4519 }
4520 drop(obj);
4521 let kphase_total_s = kphase_optim_start.elapsed().as_secs_f64();
4522 let slow_resets_end = ctx.evaluator.slow_path_reset_count();
4523 let skip_touches_end = gam_solve::pirls::nfree_skip_row_element_touches();
4524 let (search_slow_resets_end, search_skip_touches_end) =
4530 ctx.nfree_polish_boundary.unwrap_or((slow_resets_end, skip_touches_end));
4531 let kphase_slow_resets = search_slow_resets_end.saturating_sub(kphase_slow_resets_start);
4532 let kphase_polish_slow_resets = slow_resets_end.saturating_sub(search_slow_resets_end);
4533 let kphase_design_revision_delta = ctx
4534 .cache
4535 .design_revision()
4536 .saturating_sub(kphase_design_revision_start);
4537 let kphase_nfree_skip_touches =
4538 search_skip_touches_end.saturating_sub(kphase_nfree_skip_touches_start);
4539 let kphase_polish_skip_touches = skip_touches_end.saturating_sub(search_skip_touches_end);
4540 log::info!(
4541 "[KAPPA-PHASE-POLISH-SUMMARY] n_rows={} exact_polish_ran={} polish_slow_path_resets={} polish_nfree_skip_row_touches={}",
4542 data.nrows(),
4543 ctx.nfree_polish_boundary.is_some(),
4544 kphase_polish_slow_resets,
4545 kphase_polish_skip_touches,
4546 );
4547 log::info!(
4548 "[KAPPA-PHASE-SUMMARY] n_rows={} log_kappa_dim={} n_cost={} cost_total_s={:.4} n_eval={} eval_total_s={:.4} n_efs={} efs_total_s={:.4} value_realization_failures={} value_evaluation_failures={} slow_path_resets={} design_revision_delta={} nfree_skip_row_touches={} nfree_miss_shape={} nfree_miss_value={} nfree_miss_gradient={} nfree_miss_penalty={} nfree_miss_revision={} nfree_miss_second_order={} nfree_miss_other={} optim_total_s={:.4}",
4549 data.nrows(),
4550 kphase_log_kappa_dim,
4551 kphase_cost_calls.get(),
4552 kphase_cost_total_s.get(),
4553 kphase_eval_calls.get(),
4554 kphase_eval_total_s.get(),
4555 kphase_efs_calls.get(),
4556 kphase_efs_total_s.get(),
4557 ctx.value_realization_failures,
4558 ctx.value_evaluation_failures,
4559 kphase_slow_resets,
4560 kphase_design_revision_delta,
4561 kphase_nfree_skip_touches,
4562 kphase_nfree_miss_shape.get(),
4563 kphase_nfree_miss_value.get(),
4564 kphase_nfree_miss_gradient.get(),
4565 kphase_nfree_miss_penalty.get(),
4566 kphase_nfree_miss_revision.get(),
4567 kphase_nfree_miss_second_order.get(),
4568 kphase_nfree_miss_other.get(),
4569 kphase_total_s,
4570 );
4571 let timing = SpatialLengthScaleOptimizationTiming {
4572 log_kappa_dim: kphase_log_kappa_dim,
4573 cost_calls: kphase_cost_calls.get(),
4574 cost_total_s: kphase_cost_total_s.get(),
4575 eval_calls: kphase_eval_calls.get(),
4576 eval_total_s: kphase_eval_total_s.get(),
4577 efs_calls: kphase_efs_calls.get(),
4578 efs_total_s: kphase_efs_total_s.get(),
4579 slow_path_resets: kphase_slow_resets,
4580 design_revision_delta: kphase_design_revision_delta,
4581 nfree_skip_row_touches: kphase_nfree_skip_touches,
4582 nfree_miss_shape: kphase_nfree_miss_shape.get(),
4583 nfree_miss_value: kphase_nfree_miss_value.get(),
4584 nfree_miss_gradient: kphase_nfree_miss_gradient.get(),
4585 nfree_miss_penalty: kphase_nfree_miss_penalty.get(),
4586 nfree_miss_revision: kphase_nfree_miss_revision.get(),
4587 nfree_miss_second_order: kphase_nfree_miss_second_order.get(),
4588 nfree_miss_other: kphase_nfree_miss_other.get(),
4589 exact_polish_ran: ctx.nfree_polish_boundary.is_some(),
4590 polish_slow_path_resets: kphase_polish_slow_resets,
4591 polish_nfree_skip_row_touches: kphase_polish_skip_touches,
4592 optim_total_s: kphase_total_s,
4593 };
4594 log::trace!(
4595 "[{}] converged in {} iterations, final_value={:.6e}, grad_norm={}",
4596 label,
4597 result.iterations,
4598 result.final_value,
4599 result.final_grad_norm_report(),
4600 );
4601 let theta_star = result.rho;
4605 Ok((theta_star, result.final_value, seed_value, timing))
4606}
4607
4608fn set_single_term_spatial_length_scale(
4612 term: &mut SmoothTermSpec,
4613 length_scale: f64,
4614) -> Result<(), EstimationError> {
4615 match &mut term.basis {
4616 SmoothBasisSpec::ThinPlate { spec, .. } => {
4617 spec.length_scale = length_scale;
4618 Ok(())
4619 }
4620 SmoothBasisSpec::Matern { spec, .. } => {
4621 spec.length_scale.set_resolved(length_scale);
4622 Ok(())
4623 }
4624 SmoothBasisSpec::Duchon { spec, .. } => {
4625 spec.length_scale = Some(length_scale);
4626 Ok(())
4627 }
4628 _ => Err(EstimationError::InvalidInput(format!(
4629 "term '{}' does not expose a spatial length scale",
4630 term.name
4631 ))),
4632 }
4633}
4634
4635fn set_single_term_spatial_aniso_log_scales(
4639 term: &mut SmoothTermSpec,
4640 eta: Vec<f64>,
4641) -> Result<(), EstimationError> {
4642 let eta = center_aniso_log_scales(&eta);
4643 match &mut term.basis {
4644 SmoothBasisSpec::Matern { spec, .. } => {
4645 spec.aniso_log_scales = Some(eta);
4646 Ok(())
4647 }
4648 SmoothBasisSpec::Duchon { spec, .. } => {
4649 spec.aniso_log_scales = Some(eta);
4650 Ok(())
4651 }
4652 _ => Err(EstimationError::InvalidInput(format!(
4653 "term '{}' does not support aniso_log_scales",
4654 term.name
4655 ))),
4656 }
4657}
4658
4659pub fn get_constant_curvature_kappa(spec: &TermCollectionSpec, term_idx: usize) -> Option<f64> {
4678 constant_curvature_term_spec(spec, term_idx).map(|cc| cc.kappa)
4679}
4680
4681pub fn constant_curvature_kappa_is_fixed(spec: &TermCollectionSpec, term_idx: usize) -> bool {
4688 constant_curvature_term_spec(spec, term_idx).is_some_and(|cc| cc.kappa_fixed)
4689}
4690
4691pub fn constant_curvature_length_scale_is_fixed(
4699 spec: &TermCollectionSpec,
4700 term_idx: usize,
4701) -> bool {
4702 constant_curvature_term_spec(spec, term_idx).is_some_and(|cc| cc.length_scale_fixed)
4703}
4704
4705pub fn constant_curvature_term_indices(spec: &TermCollectionSpec) -> Vec<usize> {
4707 (0..spec.smooth_terms.len())
4708 .filter(|&idx| constant_curvature_term_spec(spec, idx).is_some())
4709 .collect()
4710}
4711
4712#[derive(Debug, Clone)]
4713struct SingleSmoothTermRealization {
4714 design_local: DesignMatrix,
4715 term: SmoothTerm,
4716}
4717
4718fn wrap_local_build_as_realization(
4725 mut local: LocalSmoothTermBuild,
4726 termspec: &SmoothTermSpec,
4727) -> Result<SingleSmoothTermRealization, String> {
4728 let p_local = local.dim;
4729 let lb_local = if local.box_reparam {
4730 shape_lower_bounds_local(termspec.shape, p_local)
4731 } else {
4732 None
4733 };
4734
4735 let applied_rotation: Option<gam_terms::basis::JointNullRotation> = match (
4739 local.joint_null_rotation.take(),
4740 lb_local.is_some(),
4741 local.linear_constraints.is_some(),
4742 ) {
4743 (Some(rot), false, false) => {
4744 let q = &rot.rotation;
4745 local.design =
4746 apply_smooth_transform_to_design(local.design.clone(), q, &termspec.name).map_err(
4747 |e| {
4748 format!(
4749 "joint-null absorption rotation failed for term '{}': {}",
4750 termspec.name, e
4751 )
4752 },
4753 )?;
4754 for penalty in &mut local.active_penalties {
4755 let qt_s = gam_linalg::faer_ndarray::fast_atb(q, &penalty.matrix);
4756 penalty.matrix = gam_linalg::faer_ndarray::fast_ab(&qt_s, q);
4757 penalty.null_eigenvectors = penalty
4758 .null_eigenvectors
4759 .as_ref()
4760 .map(|basis| gam_linalg::faer_ndarray::fast_atb(q, basis));
4761 penalty.info.structural_null_frame = penalty
4767 .info
4768 .structural_null_frame
4769 .as_ref()
4770 .map(|frame| gam_linalg::faer_ndarray::fast_atb(q, frame));
4771 penalty.op = None;
4772 penalty.info.kronecker_factors = None;
4773 }
4774 local.kronecker_factored = None;
4775 Some(rot)
4776 }
4777 (Some(_), _, _) => None,
4778 (None, _, _) => None,
4779 };
4780
4781 let smooth_term = SmoothTerm {
4782 parametric_residualization: None,
4783 collection_gauge: None,
4787 name: termspec.name.clone(),
4788 coeff_range: 0..p_local,
4789 shape: termspec.shape,
4790 active_penalties: local.active_penalties.clone(),
4791 dropped_penalties: local.dropped_penalties.clone(),
4792 metadata: local.metadata.clone(),
4793 lower_bounds_local: lb_local,
4794 linear_constraints_local: local.linear_constraints.clone(),
4795 kronecker_factored: local.kronecker_factored.take(),
4796 joint_null_rotation: applied_rotation,
4797 unabsorbed_global_orthogonality: None,
4800 };
4801
4802 Ok(SingleSmoothTermRealization {
4803 design_local: local.design,
4804 term: smooth_term,
4805 })
4806}
4807
4808fn freeze_geometry_from_metadata(
4819 termspec: &SmoothTermSpec,
4820 metadata: &BasisMetadata,
4821) -> Option<SmoothTermSpec> {
4822 let mut frozen = termspec.clone();
4823 match (&mut frozen.basis, metadata) {
4824 (
4825 SmoothBasisSpec::Matern {
4826 spec,
4827 input_scale: spec_scale,
4828 ..
4829 },
4830 BasisMetadata::Matern {
4831 centers,
4832 input_scale: metadata_scale,
4833 identifiability_transform,
4834 ..
4835 },
4836 ) => {
4837 spec.center_strategy = CenterStrategy::UserProvided(centers.clone());
4838 *spec_scale = Some(*metadata_scale);
4839 if let Some(transform) = identifiability_transform.clone() {
4843 spec.identifiability = MaternIdentifiability::FrozenTransform { transform };
4844 }
4845 Some(frozen)
4846 }
4847 (
4848 SmoothBasisSpec::Duchon {
4849 spec,
4850 input_scale: spec_scale,
4851 ..
4852 },
4853 BasisMetadata::Duchon {
4854 centers,
4855 input_scale: metadata_scale,
4856 ..
4857 },
4858 ) => {
4859 spec.center_strategy = CenterStrategy::UserProvided(centers.clone());
4860 *spec_scale = Some(*metadata_scale);
4861 Some(frozen)
4872 }
4873 (
4874 SmoothBasisSpec::ThinPlate {
4875 spec,
4876 input_scale: spec_scale,
4877 ..
4878 },
4879 BasisMetadata::ThinPlate {
4880 centers,
4881 input_scale: metadata_scale,
4882 ..
4883 },
4884 ) => {
4885 spec.center_strategy = CenterStrategy::UserProvided(centers.clone());
4886 *spec_scale = Some(*metadata_scale);
4887 Some(frozen)
4888 }
4889 _ => None,
4892 }
4893}
4894
4895fn restore_local_identifiability_chart(
4910 replay: &mut SmoothBasisSpec,
4911 local_chart: Option<&Array2<f64>>,
4912) {
4913 let spatial = |chart: Option<&Array2<f64>>| match chart {
4914 Some(transform) => SpatialIdentifiability::FrozenTransform {
4915 transform: transform.clone(),
4916 },
4917 None => SpatialIdentifiability::None,
4918 };
4919 if let SmoothBasisSpec::Duchon { spec, .. } = &mut *replay {
4920 spec.identifiability = spatial(local_chart);
4921 }
4922 if let SmoothBasisSpec::ThinPlate { spec, .. } = &mut *replay {
4923 spec.identifiability = spatial(local_chart);
4924 }
4925 if let SmoothBasisSpec::Matern { spec, .. } = &mut *replay {
4926 spec.identifiability = match local_chart {
4927 Some(transform) => MaternIdentifiability::FrozenTransform {
4928 transform: transform.clone(),
4929 },
4930 None => MaternIdentifiability::None,
4931 };
4932 }
4933 if let (SmoothBasisSpec::ConstantCurvature { spec, .. }, Some(transform)) =
4939 (&mut *replay, local_chart)
4940 {
4941 spec.identifiability = gam_terms::basis::ConstantCurvatureIdentifiability::FrozenTransform {
4942 transform: transform.clone(),
4943 };
4944 }
4945 if let (SmoothBasisSpec::MeasureJet { spec, .. }, Some(transform)) = (&mut *replay, local_chart)
4946 {
4947 spec.identifiability = gam_terms::basis::MeasureJetIdentifiability::FrozenTransform {
4948 transform: transform.clone(),
4949 };
4950 }
4951}
4952
4953fn spatial_frozen_radial_chart_shape(termspec: &SmoothTermSpec) -> Option<(usize, usize)> {
4955 match &termspec.basis {
4956 SmoothBasisSpec::Duchon { spec, .. } => spec.radial_reparam.as_ref().map(|v| v.dim()),
4957 SmoothBasisSpec::ThinPlate { spec, .. } => spec.radial_reparam.as_ref().map(|v| v.dim()),
4958 _ => None,
4959 }
4960}
4961
4962fn spatial_realized_radial_chart_shape(metadata: &BasisMetadata) -> Option<(usize, usize)> {
4964 match metadata {
4965 BasisMetadata::Duchon { radial_reparam, .. } => radial_reparam.as_ref().map(|v| v.dim()),
4966 BasisMetadata::ThinPlate { radial_reparam, .. } => radial_reparam.as_ref().map(|v| v.dim()),
4967 _ => None,
4968 }
4969}
4970
4971fn rebuild_smooth_auxiliary_state(
4972 smooth: &mut SmoothDesign,
4973 dropped_penaltyinfo_by_term: &[Vec<DroppedPenaltyBlockInfo>],
4974) -> Result<(), String> {
4975 if dropped_penaltyinfo_by_term.len() != smooth.terms.len() {
4976 return Err(SmoothError::dimension_mismatch(format!(
4977 "smooth dropped-penalty cache mismatch: terms={}, dropped_sets={}",
4978 smooth.terms.len(),
4979 dropped_penaltyinfo_by_term.len()
4980 ))
4981 .into());
4982 }
4983
4984 let total_p = smooth.total_smooth_cols();
4985 let mut coefficient_lower_bounds = Array1::<f64>::from_elem(total_p, f64::NEG_INFINITY);
4986 let mut any_bounds = false;
4987 let mut linear_constraintrows: Vec<Array1<f64>> = Vec::new();
4988 let mut linear_constraint_b: Vec<f64> = Vec::new();
4989
4990 for term in &smooth.terms {
4991 let range = term.coeff_range.clone();
4992 if let Some(lb_local) = term.lower_bounds_local.as_ref() {
4993 if lb_local.len() != range.len() {
4994 return Err(SmoothError::dimension_mismatch(format!(
4995 "smooth lower-bound cache mismatch for term '{}': bounds={}, coeffs={}",
4996 term.name,
4997 lb_local.len(),
4998 range.len()
4999 ))
5000 .into());
5001 }
5002 coefficient_lower_bounds
5003 .slice_mut(s![range.clone()])
5004 .assign(lb_local);
5005 any_bounds = true;
5006 }
5007 if let Some(lin_local) = term.linear_constraints_local.as_ref() {
5008 if lin_local.a.ncols() != range.len() {
5009 return Err(SmoothError::dimension_mismatch(format!(
5010 "smooth linear-constraint cache mismatch for term '{}': cols={}, coeffs={}",
5011 term.name,
5012 lin_local.a.ncols(),
5013 range.len()
5014 ))
5015 .into());
5016 }
5017 for r in 0..lin_local.a.nrows() {
5018 let mut row = Array1::<f64>::zeros(total_p);
5019 row.slice_mut(s![range.clone()]).assign(&lin_local.a.row(r));
5020 linear_constraintrows.push(row);
5021 linear_constraint_b.push(lin_local.b[r]);
5022 }
5023 }
5024 }
5025
5026 smooth.coefficient_lower_bounds = if any_bounds {
5027 Some(coefficient_lower_bounds)
5028 } else {
5029 None
5030 };
5031 smooth.linear_constraints = if linear_constraintrows.is_empty() {
5032 None
5033 } else {
5034 let mut a = Array2::<f64>::zeros((linear_constraintrows.len(), total_p));
5035 for (i, row) in linear_constraintrows.iter().enumerate() {
5036 a.row_mut(i).assign(row);
5037 }
5038 Some(LinearInequalityConstraints {
5039 a,
5040 b: Array1::from_vec(linear_constraint_b),
5041 })
5042 };
5043 smooth.dropped_penaltyinfo = dropped_penaltyinfo_by_term
5044 .iter()
5045 .flat_map(|infos| infos.iter().cloned())
5046 .collect();
5047 Ok(())
5048}
5049
5050fn rebuild_term_collection_auxiliary_state(
5051 spec: &TermCollectionSpec,
5052 design: &mut TermCollectionDesign,
5053) -> Result<(), String> {
5054 if spec.linear_terms.len() != design.linear_ranges.len() {
5055 return Err(SmoothError::dimension_mismatch(format!(
5056 "term-collection linear bookkeeping mismatch: spec_terms={}, design_ranges={}",
5057 spec.linear_terms.len(),
5058 design.linear_ranges.len()
5059 ))
5060 .into());
5061 }
5062
5063 let p_total = design.design.ncols();
5064 let smooth_start = p_total.saturating_sub(design.smooth.total_smooth_cols());
5065 let mut coefficient_lower_bounds = Array1::<f64>::from_elem(p_total, f64::NEG_INFINITY);
5066 let mut any_bounds = false;
5067 let mut linear_constraintrows: Vec<Array1<f64>> = Vec::new();
5068 let mut linear_constraint_b: Vec<f64> = Vec::new();
5069
5070 for (linear, (_, range)) in spec.linear_terms.iter().zip(design.linear_ranges.iter()) {
5071 if range.len() != 1 {
5072 return Err(SmoothError::dimension_mismatch(format!(
5073 "linear term '{}' expected one coefficient column, found {}",
5074 linear.name,
5075 range.len()
5076 ))
5077 .into());
5078 }
5079 let col = range.start;
5080 if let Some(lb) = linear.coefficient_min {
5081 let mut row = Array1::<f64>::zeros(p_total);
5082 row[col] = 1.0;
5083 linear_constraintrows.push(row);
5084 linear_constraint_b.push(lb);
5085 }
5086 if let Some(ub) = linear.coefficient_max {
5087 let mut row = Array1::<f64>::zeros(p_total);
5088 row[col] = -1.0;
5089 linear_constraintrows.push(row);
5090 linear_constraint_b.push(-ub);
5091 }
5092 }
5093
5094 if let Some(lb_smooth) = design.smooth.coefficient_lower_bounds.as_ref() {
5095 if lb_smooth.len() != design.smooth.total_smooth_cols() {
5096 return Err(SmoothError::dimension_mismatch(format!(
5097 "smooth lower-bound width mismatch: bounds={}, smooth_cols={}",
5098 lb_smooth.len(),
5099 design.smooth.total_smooth_cols()
5100 ))
5101 .into());
5102 }
5103 coefficient_lower_bounds
5104 .slice_mut(s![
5105 smooth_start..(smooth_start + design.smooth.total_smooth_cols())
5106 ])
5107 .assign(lb_smooth);
5108 any_bounds = true;
5109 }
5110 if let Some(lin_smooth) = design.smooth.linear_constraints.as_ref() {
5111 if lin_smooth.a.ncols() != design.smooth.total_smooth_cols() {
5112 return Err(SmoothError::dimension_mismatch(format!(
5113 "smooth linear-constraint width mismatch: cols={}, smooth_cols={}",
5114 lin_smooth.a.ncols(),
5115 design.smooth.total_smooth_cols()
5116 ))
5117 .into());
5118 }
5119 let mut a_global = Array2::<f64>::zeros((lin_smooth.a.nrows(), p_total));
5120 a_global
5121 .slice_mut(s![
5122 ..,
5123 smooth_start..(smooth_start + design.smooth.total_smooth_cols())
5124 ])
5125 .assign(&lin_smooth.a);
5126 for r in 0..a_global.nrows() {
5127 linear_constraintrows.push(a_global.row(r).to_owned());
5128 linear_constraint_b.push(lin_smooth.b[r]);
5129 }
5130 }
5131
5132 let lower_bound_constraints = if any_bounds {
5133 linear_constraints_from_lower_bounds_global(&coefficient_lower_bounds)
5134 } else {
5135 None
5136 };
5137 let explicit_linear_constraints = if linear_constraintrows.is_empty() {
5138 None
5139 } else {
5140 let mut a = Array2::<f64>::zeros((linear_constraintrows.len(), p_total));
5141 for (i, row) in linear_constraintrows.iter().enumerate() {
5142 a.row_mut(i).assign(row);
5143 }
5144 Some(LinearInequalityConstraints {
5145 a,
5146 b: Array1::from_vec(linear_constraint_b),
5147 })
5148 };
5149
5150 design.coefficient_lower_bounds = if any_bounds {
5151 Some(coefficient_lower_bounds)
5152 } else {
5153 None
5154 };
5155 design.linear_constraints =
5156 merge_linear_constraints_global(explicit_linear_constraints, lower_bound_constraints)
5157 .map_err(|error| error.to_string())?;
5158 design.dropped_penaltyinfo = design.smooth.dropped_penaltyinfo.clone();
5159 Ok(())
5160}
5161
5162fn theta_values_match(left: &Array1<f64>, right: &Array1<f64>) -> bool {
5163 left.len() == right.len()
5164 && left
5165 .iter()
5166 .zip(right.iter())
5167 .all(|(&l, &r)| l.to_bits() == r.to_bits())
5168}
5169
5170fn latent_values_match(left: &Array1<f64>, right: &Array1<f64>) -> bool {
5171 theta_values_match(left, right)
5172}
5173
5174fn spatial_aniso_matches(left: Option<&[f64]>, right: Option<&[f64]>) -> bool {
5175 match (left, right) {
5176 (None, None) => true,
5177 (Some(a), Some(b)) => {
5178 a.len() == b.len()
5179 && a.iter()
5180 .zip(b.iter())
5181 .all(|(&x, &y)| x.to_bits() == y.to_bits())
5182 }
5183 _ => false,
5184 }
5185}
5186
5187fn spatial_length_scale_matches(left: Option<f64>, right: Option<f64>) -> bool {
5188 match (left, right) {
5189 (None, None) => true,
5190 (Some(a), Some(b)) => a.to_bits() == b.to_bits(),
5191 _ => false,
5192 }
5193}
5194
5195struct FrozenTermCollectionIncrementalRealizer<'d> {
5196 data: ArrayView2<'d, f64>,
5197 spec: TermCollectionSpec,
5198 design: TermCollectionDesign,
5199 fixed_blocks: Vec<DesignBlock>,
5200 dropped_penaltyinfo_by_term: Vec<Vec<DroppedPenaltyBlockInfo>>,
5201 smooth_penalty_ranges: Vec<Range<usize>>,
5202 full_penalty_ranges: Vec<Range<usize>>,
5203 basisworkspace: gam_terms::basis::BasisWorkspace,
5207 spatial_realization_geometry: Vec<Option<SmoothTermSpec>>,
5220 design_revision: u64,
5226}
5227
5228impl<'d> std::fmt::Debug for FrozenTermCollectionIncrementalRealizer<'d> {
5229 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
5230 f.debug_struct("FrozenTermCollectionIncrementalRealizer")
5231 .field("data_shape", &(self.data.nrows(), self.data.ncols()))
5232 .field("fixed_blocks", &self.fixed_blocks.len())
5233 .finish_non_exhaustive()
5234 }
5235}
5236
5237fn emitted_smooth_penalty_ranges(
5246 design: &TermCollectionDesign,
5247) -> Result<(Vec<Range<usize>>, Vec<Range<usize>>), String> {
5248 let leading = design.leading_penalty_blocks_before_smooth();
5249 let mut smooth_penalty_ranges = Vec::with_capacity(design.smooth.terms.len());
5250 let mut full_penalty_ranges = Vec::with_capacity(design.smooth.terms.len());
5251 let mut smooth_cursor = 0usize;
5252 for term_idx in 0..design.smooth.terms.len() {
5253 let full_range = design.smooth_term_penalty_range(term_idx)?;
5254 match full_range {
5255 Some(full_range) => {
5256 let local_start = full_range.start.checked_sub(leading).ok_or_else(|| {
5257 "incremental realizer smooth penalty range precedes the emitted smooth prefix"
5258 .to_string()
5259 })?;
5260 let local_end = full_range.end.checked_sub(leading).ok_or_else(|| {
5261 "incremental realizer smooth penalty range precedes the emitted smooth prefix"
5262 .to_string()
5263 })?;
5264 if local_start != smooth_cursor {
5265 return Err(format!(
5266 "incremental realizer non-contiguous emitted smooth layout at term {term_idx}: expected local start {smooth_cursor}, got {local_start}"
5267 ));
5268 }
5269 smooth_cursor = local_end;
5270 smooth_penalty_ranges.push(local_start..local_end);
5271 full_penalty_ranges.push(full_range);
5272 }
5273 None => {
5274 smooth_penalty_ranges.push(smooth_cursor..smooth_cursor);
5275 let global_cursor = leading.checked_add(smooth_cursor).ok_or_else(|| {
5276 "incremental realizer empty smooth penalty range overflow".to_string()
5277 })?;
5278 full_penalty_ranges.push(global_cursor..global_cursor);
5279 }
5280 }
5281 }
5282 if smooth_cursor != design.smooth.penalties.len() {
5283 return Err(format!(
5284 "incremental realizer smooth penalty mismatch: ranged={}, actual={}",
5285 smooth_cursor,
5286 design.smooth.penalties.len()
5287 ));
5288 }
5289 Ok((smooth_penalty_ranges, full_penalty_ranges))
5290}
5291
5292impl<'d> FrozenTermCollectionIncrementalRealizer<'d> {
5293 fn new(
5294 data: ArrayView2<'d, f64>,
5295 spec: TermCollectionSpec,
5296 design: TermCollectionDesign,
5297 ) -> Result<Self, String> {
5298 let policy = gam_runtime::resource::ResourcePolicy::default_library();
5299 Self::new_with_policy(data, spec, design, &policy)
5300 }
5301
5302 fn new_with_policy(
5303 data: ArrayView2<'d, f64>,
5304 spec: TermCollectionSpec,
5305 design: TermCollectionDesign,
5306 policy: &gam_runtime::resource::ResourcePolicy,
5307 ) -> Result<Self, String> {
5308 if spec.smooth_terms.len() != design.smooth.terms.len() {
5309 return Err(SmoothError::dimension_mismatch(format!(
5310 "incremental realizer smooth term mismatch: spec_terms={}, design_terms={}",
5311 spec.smooth_terms.len(),
5312 design.smooth.terms.len()
5313 ))
5314 .into());
5315 }
5316
5317 let (smooth_penalty_ranges, full_penalty_ranges) = emitted_smooth_penalty_ranges(&design)?;
5322 let mut spec = freeze_term_collection_from_design(&spec, &design)
5371 .map_err(|e| format!("failed to freeze incremental replay specification: {e}"))?;
5372 for (term_idx, term) in design.smooth.terms.iter().enumerate() {
5373 let Some(gauge) = term.collection_gauge.as_ref() else {
5374 continue;
5375 };
5376 let Some(replay) = spec.smooth_terms.get_mut(term_idx) else {
5377 continue;
5378 };
5379 restore_local_identifiability_chart(
5380 &mut replay.basis,
5381 gauge.local_identifiability_transform.as_ref(),
5382 );
5383 }
5384 let spec = spec;
5385 let fixed_blocks = build_term_collection_fixed_blocks(data, &spec)
5386 .map_err(|e| format!("failed to cache fixed term-collection blocks: {e}"))?;
5387
5388 let dropped_penaltyinfo_by_term: Vec<Vec<DroppedPenaltyBlockInfo>> = design
5402 .smooth
5403 .terms
5404 .iter()
5405 .map(|term| {
5406 term.dropped_penalties
5407 .iter()
5408 .cloned()
5409 .map(|penalty| DroppedPenaltyBlockInfo {
5410 termname: Some(term.name.clone()),
5411 penalty,
5412 })
5413 .collect()
5414 })
5415 .collect();
5416
5417 let geometry_slots = spec.smooth_terms.len();
5418 Ok(Self {
5419 data,
5420 spec,
5421 design,
5422 fixed_blocks,
5423 dropped_penaltyinfo_by_term,
5424 smooth_penalty_ranges,
5425 full_penalty_ranges,
5426 basisworkspace: gam_terms::basis::BasisWorkspace::with_policy(policy.clone()),
5427 spatial_realization_geometry: vec![None; geometry_slots],
5428 design_revision: 0,
5429 })
5430 }
5431
5432 fn design_revision(&self) -> u64 {
5433 self.design_revision
5434 }
5435
5436 fn spec(&self) -> &TermCollectionSpec {
5437 &self.spec
5438 }
5439
5440 fn design(&self) -> &TermCollectionDesign {
5441 &self.design
5442 }
5443
5444 fn supports_nfree_penalty_rekey(&self, spatial_terms: &[usize]) -> bool {
5485 if spatial_terms.len() != 1 {
5486 return false;
5487 }
5488 let term_idx = spatial_terms[0];
5489 matches!(
5490 self.design.smooth.terms.get(term_idx).map(|t| &t.metadata),
5491 Some(BasisMetadata::Duchon { .. } | BasisMetadata::ThinPlate { .. })
5492 )
5493 }
5494
5495 fn supports_nfree_gradient_only_routing(&self, spatial_terms: &[usize]) -> bool {
5504 if spatial_terms.len() != 1 {
5505 return false;
5506 }
5507 let term_idx = spatial_terms[0];
5508 matches!(
5509 self.design.smooth.terms.get(term_idx).map(|t| &t.metadata),
5510 Some(BasisMetadata::Duchon { .. } | BasisMetadata::ThinPlate { .. })
5511 )
5512 }
5513
5514 fn canonical_penalties_at_psi(
5527 &mut self,
5528 spatial_terms: &[usize],
5529 psi: &[f64],
5530 ) -> Result<(Vec<gam_terms::construction::CanonicalPenalty>, Vec<usize>), String> {
5531 if spatial_terms.len() != 1 {
5532 return Err(format!(
5533 "n-free penalty re-key requires exactly one spatial term, found {}",
5534 spatial_terms.len()
5535 ));
5536 }
5537 let term_idx = spatial_terms[0];
5538 let (ls_opt, aniso_from_psi) = spatial_term_psi_to_length_scale_and_aniso(psi);
5544 let termspec =
5547 self.spec.smooth_terms.get(term_idx).ok_or_else(|| {
5548 format!("spatial term {term_idx} out of range for n-free penalty")
5549 })?;
5550 let term = self
5551 .design
5552 .smooth
5553 .terms
5554 .get(term_idx)
5555 .ok_or_else(|| format!("realized smooth term {term_idx} out of range"))?;
5556 let p_total = self.design.design.ncols();
5559 let (locals, nullspace_dims): (Vec<Array2<f64>>, Vec<usize>) = match &term.metadata {
5560 BasisMetadata::Duchon {
5561 centers,
5562 identifiability_transform,
5563 operator_collocation_points,
5564 power,
5565 nullspace_order,
5566 aniso_log_scales,
5567 input_scale,
5568 radial_reparam,
5569 ..
5570 } => {
5571 let operator_penalties = match &termspec.basis {
5572 SmoothBasisSpec::Duchon { spec, .. } => spec.operator_penalties.clone(),
5573 _ => gam_terms::basis::DuchonOperatorPenaltySpec::default(),
5574 };
5575 let effective_ls = ls_opt.map(|length| {
5582 input_scale
5583 .to_standardized_units(gam_terms::OriginalUnits::new(length))
5584 .standardized_value()
5585 });
5586 gam_terms::basis::duchon_penalties_at_length_scale(
5587 centers.view(),
5588 identifiability_transform.as_ref(),
5589 operator_collocation_points.as_ref().map(|p| p.view()),
5590 &operator_penalties,
5591 *power,
5592 *nullspace_order,
5593 aniso_log_scales.as_deref(),
5594 radial_reparam.as_ref(),
5595 effective_ls,
5596 &mut self.basisworkspace,
5597 )
5598 .map_err(|e| e.to_string())?
5599 }
5600 BasisMetadata::Matern {
5601 centers,
5602 periodic,
5603 nu,
5604 include_intercept,
5605 identifiability_transform,
5606 aniso_log_scales,
5607 input_scale,
5608 ..
5609 } => {
5610 let ls = ls_opt.ok_or_else(|| {
5617 "Matérn n-free penalty re-key requires a finite length-scale".to_string()
5618 })?;
5619 let effective_ls = input_scale
5620 .to_standardized_units(gam_terms::OriginalUnits::new(ls))
5621 .standardized_value();
5622 let aniso_for_penalty = aniso_from_psi.as_deref().or(aniso_log_scales.as_deref());
5623 let filtered = matern_operator_penalty_triplet_at_length_scale(
5634 centers.view(),
5635 periodic.as_deref(),
5636 identifiability_transform.as_ref(),
5637 *nu,
5638 *include_intercept,
5639 aniso_for_penalty,
5640 effective_ls,
5641 )
5642 .map_err(|e| e.to_string())?;
5643 let locals = filtered
5644 .active
5645 .iter()
5646 .map(|penalty| penalty.matrix.clone())
5647 .collect();
5648 let nullspace_dims = filtered
5649 .active
5650 .iter()
5651 .map(|penalty| penalty.nullity)
5652 .collect();
5653 (locals, nullspace_dims)
5654 }
5655 BasisMetadata::ThinPlate {
5656 centers,
5657 identifiability_transform,
5658 radial_reparam,
5659 ..
5660 } => {
5661 let ls = ls_opt.ok_or_else(|| {
5662 "thin-plate n-free penalty re-key requires a finite length-scale".to_string()
5663 })?;
5664 let double_penalty = match &termspec.basis {
5665 SmoothBasisSpec::ThinPlate { spec, .. } => spec.double_penalty,
5666 _ => false,
5667 };
5668 gam_terms::basis::thin_plate_penalties_at_length_scale(
5669 centers.view(),
5670 identifiability_transform.as_ref(),
5671 radial_reparam.as_ref(),
5672 ls,
5673 double_penalty,
5674 &mut self.basisworkspace,
5675 )
5676 .map_err(|e| e.to_string())?
5677 }
5678 other => {
5679 return Err(format!(
5680 "n-free penalty re-key unsupported for basis metadata {:?}",
5681 std::mem::discriminant(other)
5682 ));
5683 }
5684 };
5685 let templates = &self.design.penalties;
5690 if templates.len() != locals.len() {
5691 return Err(format!(
5692 "n-free penalty re-key produced {} blocks but the frozen design carries {} \
5693 — penalty topology is not ψ-stable",
5694 locals.len(),
5695 templates.len()
5696 ));
5697 }
5698 let specs: Vec<gam_solve::estimate::PenaltySpec> = templates
5699 .iter()
5700 .zip(locals.into_iter())
5701 .map(|(tmpl, local)| gam_solve::estimate::PenaltySpec::Block {
5702 local,
5703 col_range: tmpl.col_range.clone(),
5704 prior_mean: tmpl.prior_mean.clone(),
5705 structure_hint: tmpl.structure_hint.clone(),
5706 op: tmpl.op.clone(),
5707 })
5708 .collect();
5709 gam_terms::construction::canonicalize_penalty_specs(
5710 &specs,
5711 &nullspace_dims,
5712 p_total,
5713 "nfree-psi-penalty",
5714 )
5715 .map_err(|e| e.to_string())
5716 }
5717
5718 fn canonical_penalty_derivatives_at_psi(
5719 &mut self,
5720 spatial_terms: &[usize],
5721 psi: &[f64],
5722 ) -> Result<(Range<usize>, usize, Vec<Array2<f64>>), String> {
5723 if spatial_terms.len() != 1 {
5724 return Err(format!(
5725 "n-free penalty derivative re-key requires exactly one spatial term, found {}",
5726 spatial_terms.len()
5727 ));
5728 }
5729 let term_idx = spatial_terms[0];
5730 let (ls_opt, aniso_from_psi) = spatial_term_psi_to_length_scale_and_aniso(psi);
5731 let termspec = self.spec.smooth_terms.get(term_idx).ok_or_else(|| {
5732 format!("spatial term {term_idx} out of range for n-free penalty derivative")
5733 })?;
5734 let term = self
5735 .design
5736 .smooth
5737 .terms
5738 .get(term_idx)
5739 .ok_or_else(|| format!("realized smooth term {term_idx} out of range"))?;
5740 let p_total = self.design.design.ncols();
5741 let smooth_start = p_total.saturating_sub(self.design.smooth.total_smooth_cols());
5742 let global_range =
5743 (smooth_start + term.coeff_range.start)..(smooth_start + term.coeff_range.end);
5744
5745 let locals = match &term.metadata {
5746 BasisMetadata::Duchon {
5747 centers,
5748 identifiability_transform,
5749 operator_collocation_points,
5750 power,
5751 nullspace_order,
5752 aniso_log_scales,
5753 input_scale,
5754 radial_reparam,
5755 ..
5756 } => {
5757 let mut spec = match &termspec.basis {
5758 SmoothBasisSpec::Duchon { spec, .. } => spec.clone(),
5759 _ => {
5760 return Err(
5761 "Duchon n-free penalty derivative requires a Duchon term spec"
5762 .to_string(),
5763 );
5764 }
5765 };
5766 let effective_ls = ls_opt.map(|length| {
5767 input_scale
5768 .to_standardized_units(gam_terms::OriginalUnits::new(length))
5769 .standardized_value()
5770 });
5771 spec.length_scale = effective_ls;
5772 spec.power = *power;
5773 spec.nullspace_order = *nullspace_order;
5774 spec.aniso_log_scales = aniso_log_scales.clone();
5775 spec.radial_reparam = radial_reparam.clone();
5778 if spec.length_scale.is_none() {
5779 return Err(
5780 "Duchon n-free penalty derivative requires a hybrid length-scale"
5781 .to_string(),
5782 );
5783 }
5784 let collocation = operator_collocation_points
5785 .as_ref()
5786 .map(|points| points.view())
5787 .unwrap_or_else(|| centers.view());
5788 let (_native_sources, mut first, _native_second) =
5789 gam_terms::basis::build_duchon_native_penalty_psi_derivatives(
5790 centers.view(),
5791 &spec,
5792 identifiability_transform.as_ref(),
5793 &mut self.basisworkspace,
5794 )
5795 .map_err(|e| e.to_string())?;
5796 let (_operator_sources, operator_first, _operator_second) =
5797 gam_terms::basis::build_duchon_operator_penalty_psi_derivatives(
5798 collocation,
5799 centers.view(),
5800 &spec,
5801 identifiability_transform.as_ref(),
5802 &mut self.basisworkspace,
5803 )
5804 .map_err(|e| e.to_string())?;
5805 first.extend(operator_first);
5806 first
5807 }
5808 BasisMetadata::Matern {
5809 centers,
5810 periodic,
5811 nu,
5812 include_intercept,
5813 identifiability_transform,
5814 aniso_log_scales,
5815 input_scale,
5816 ..
5817 } => {
5818 let ls = ls_opt.ok_or_else(|| {
5819 "Matérn n-free penalty derivative requires a finite length-scale".to_string()
5820 })?;
5821 let effective_ls = input_scale
5822 .to_standardized_units(gam_terms::OriginalUnits::new(ls))
5823 .standardized_value();
5824 let penalty_centers = gam_terms::basis::expand_periodic_centers(
5825 ¢ers.to_owned(),
5826 periodic.as_deref(),
5827 )
5828 .map_err(|e| e.to_string())?;
5829 let aniso_for_penalty = aniso_from_psi.as_deref().or(aniso_log_scales.as_deref());
5830 let (first, _second) =
5831 gam_terms::basis::build_matern_operator_penalty_psi_derivatives(
5832 penalty_centers.view(),
5833 effective_ls,
5834 *nu,
5835 *include_intercept,
5836 identifiability_transform.as_ref(),
5837 aniso_for_penalty,
5838 )
5839 .map_err(|e| e.to_string())?;
5840 first
5841 }
5842 BasisMetadata::ThinPlate {
5843 centers,
5844 identifiability_transform,
5845 radial_reparam,
5846 ..
5847 } => {
5848 let ls = ls_opt.ok_or_else(|| {
5849 "thin-plate n-free penalty derivative requires a finite length-scale"
5850 .to_string()
5851 })?;
5852 let mut spec = match &termspec.basis {
5853 SmoothBasisSpec::ThinPlate { spec, .. } => spec.clone(),
5854 _ => {
5855 return Err(
5856 "thin-plate n-free penalty derivative requires a ThinPlate term spec"
5857 .to_string(),
5858 );
5859 }
5860 };
5861 spec.length_scale = ls;
5862 if spec.radial_reparam.is_none() {
5863 spec.radial_reparam = radial_reparam.clone();
5864 }
5865 let (primary, _primary_second, nullspace, _nullspace_second) =
5866 gam_terms::basis::build_thin_plate_penalty_psi_derivativeswithworkspace(
5867 centers.view(),
5868 &spec,
5869 identifiability_transform.as_ref(),
5870 &mut self.basisworkspace,
5871 )
5872 .map_err(|e| e.to_string())?;
5873 if self.design.penalties.len() > 1 {
5874 vec![primary, nullspace]
5875 } else {
5876 vec![primary]
5877 }
5878 }
5879 other => {
5880 return Err(format!(
5881 "n-free penalty derivative re-key unsupported for basis metadata {:?}",
5882 std::mem::discriminant(other)
5883 ));
5884 }
5885 };
5886 if locals.len() != self.design.penalties.len() {
5887 return Err(format!(
5888 "n-free penalty derivative re-key produced {} blocks but the frozen design carries {} \
5889 — penalty topology is not ψ-stable",
5890 locals.len(),
5891 self.design.penalties.len()
5892 ));
5893 }
5894 Ok((global_range, p_total, locals))
5895 }
5896
5897 fn apply_log_kappa(
5907 &mut self,
5908 log_kappa: &SpatialLogKappaCoords,
5909 term_indices: &[usize],
5910 ) -> Result<(), EstimationError> {
5911 if term_indices.len() != log_kappa.dims_per_term().len() {
5912 return Err(EstimationError::InvalidInput(
5913 SmoothError::dimension_mismatch(format!(
5914 "incremental realizer log-kappa term mismatch: term_indices={}, dims_per_term={}",
5915 term_indices.len(),
5916 log_kappa.dims_per_term().len()
5917 ))
5918 .to_string(),
5919 ));
5920 }
5921
5922 let mut any_changed = false;
5923 for (slot, &term_idx) in term_indices.iter().enumerate() {
5924 any_changed |= self.apply_log_kappa_to_term(term_idx, log_kappa.term_slice(slot))?;
5925 }
5926
5927 if any_changed {
5928 self.refresh_full_design_operator()
5929 .map_err(EstimationError::InvalidInput)?;
5930 rebuild_smooth_auxiliary_state(
5931 &mut self.design.smooth,
5932 &self.dropped_penaltyinfo_by_term,
5933 )
5934 .map_err(EstimationError::InvalidInput)?;
5935 rebuild_term_collection_auxiliary_state(&self.spec, &mut self.design)
5936 .map_err(EstimationError::InvalidInput)?;
5937 self.design_revision = self.design_revision.wrapping_add(1);
5938 }
5939 Ok(())
5940 }
5941
5942 fn apply_log_kappa_to_term(
5943 &mut self,
5944 term_idx: usize,
5945 psi: &[f64],
5946 ) -> Result<bool, EstimationError> {
5947 if !spatial_term_supports_hyper_optimization(&self.spec, term_idx) {
5948 return Err(EstimationError::InvalidInput(
5949 SmoothError::invalid_config(format!(
5950 "incremental realizer term {term_idx} does not expose spatial hyperparameters"
5951 ))
5952 .to_string(),
5953 ));
5954 }
5955 let measure_jet_term = measure_jet_term_spec(&self.spec, term_idx).is_some();
5959 let constant_curvature_term = constant_curvature_term_spec(&self.spec, term_idx).is_some();
5963 let mut next_length_scale = None;
5964 let mut next_aniso: Option<Vec<f64>> = None;
5965 if measure_jet_term {
5966 if !set_measure_jet_psi_dials(&mut self.spec, term_idx, psi)
5967 ?
5968 {
5969 return Ok(false);
5970 }
5971 } else if constant_curvature_term {
5972 if !set_constant_curvature_kappa(&mut self.spec, term_idx, psi)
5973 ?
5974 {
5975 return Ok(false);
5976 }
5977 } else {
5978 let current_length_scale = get_spatial_length_scale(&self.spec, term_idx);
5979 let current_aniso = get_spatial_aniso_log_scales(&self.spec, term_idx);
5980 let (ls, eta) = spatial_term_psi_to_length_scale_and_aniso(psi);
5981 next_length_scale = ls;
5982 next_aniso = eta;
5983 let same_length = spatial_length_scale_matches(current_length_scale, next_length_scale);
5984 let same_aniso = spatial_aniso_matches(current_aniso.as_deref(), next_aniso.as_deref());
5985 if same_length && same_aniso {
5986 return Ok(false);
5987 }
5988 if let Some(length_scale) = next_length_scale {
5989 set_spatial_length_scale(&mut self.spec, term_idx, length_scale)
5990 ?;
5991 }
5992 if let Some(eta) = next_aniso.clone() {
5993 set_spatial_aniso_log_scales(&mut self.spec, term_idx, eta)
5994 ?;
5995 }
5996 }
5997
5998 let geometry_slot = self
6009 .spatial_realization_geometry
6010 .get(term_idx)
6011 .ok_or_else(|| EstimationError::InvalidInput(format!("incremental realizer geometry slot {term_idx} out of range")))?;
6012 let geometry_cached = geometry_slot.is_some();
6013 let mut build_spec = match geometry_slot {
6014 Some(cached) => cached.clone(),
6015 None => self
6016 .spec
6017 .smooth_terms
6018 .get(term_idx)
6019 .ok_or_else(|| EstimationError::InvalidInput(format!("incremental realizer smooth term {term_idx} out of range")))?
6020 .clone(),
6021 };
6022 if measure_jet_term {
6023 set_single_term_measure_jet_psi_dials(&mut build_spec, psi)
6027 ?;
6028 } else if constant_curvature_term {
6029 set_single_term_constant_curvature_kappa(&mut build_spec, psi)
6034 ?;
6035 } else {
6036 if let Some(length_scale) = next_length_scale {
6037 set_single_term_spatial_length_scale(&mut build_spec, length_scale)
6038 ?;
6039 }
6040 if let Some(eta) = next_aniso {
6041 set_single_term_spatial_aniso_log_scales(&mut build_spec, eta)
6042 ?;
6043 }
6044 }
6045
6046 let termname = build_spec.name.clone();
6047 let local = build_single_local_smooth_term(
6048 self.data,
6049 &build_spec,
6050 &mut self.basisworkspace,
6051 )
6052 .map_err(|e| {
6053 EstimationError::InvalidInput(format!(
6054 "failed to rebuild smooth term '{termname}' during incremental κ realization: {e}"
6055 ))
6056 })?;
6057
6058 if self.spatial_realization_geometry[term_idx].is_none()
6063 && let Some(frozen) = freeze_geometry_from_metadata(&build_spec, &local.metadata)
6064 {
6065 if let (
6077 SmoothBasisSpec::Matern {
6078 spec: frozen_spec, ..
6079 },
6080 Some(SmoothBasisSpec::Matern {
6081 spec: live_spec, ..
6082 }),
6083 ) = (
6084 &frozen.basis,
6085 self.spec
6086 .smooth_terms
6087 .get_mut(term_idx)
6088 .map(|t| &mut t.basis),
6089 ) {
6090 live_spec.identifiability = frozen_spec.identifiability.clone();
6091 live_spec.center_strategy = frozen_spec.center_strategy.clone();
6092 }
6093 self.spatial_realization_geometry[term_idx] = Some(frozen);
6094 }
6095
6096 let trial_report = format!(
6099 "psi={psi:?}, length_scale={next_length_scale:?}, geometry_cached={geometry_cached}, \
6100 frozen_radial_chart={:?}, realized_radial_chart={:?}, local_cols={}",
6101 spatial_frozen_radial_chart_shape(&build_spec),
6102 spatial_realized_radial_chart_shape(&local.metadata),
6103 local.design.ncols(),
6104 );
6105 let realization = wrap_local_build_as_realization(local, &build_spec)
6106 .map_err(EstimationError::InvalidInput)?;
6107 self.replace_term_realization(term_idx, realization, &trial_report)?;
6108 Ok(true)
6109 }
6110
6111 fn replace_term_realization(
6112 &mut self,
6113 term_idx: usize,
6114 realization: SingleSmoothTermRealization,
6115 trial_report: &str,
6116 ) -> Result<(), EstimationError> {
6117 let t_replace = std::time::Instant::now();
6118 let SingleSmoothTermRealization { design_local, term } = realization;
6119 let SmoothTerm {
6120 name,
6121 active_penalties,
6122 dropped_penalties,
6123 metadata,
6124 lower_bounds_local,
6125 linear_constraints_local,
6126 joint_null_rotation,
6127 ..
6128 } = term;
6129 let collection_gauge = self
6149 .design
6150 .smooth
6151 .terms
6152 .get(term_idx)
6153 .and_then(|target| target.collection_gauge.clone());
6154 let pre_gauge_cols = design_local.ncols();
6160 let gauge_report = match collection_gauge.as_ref() {
6161 Some(gauge) => format!(
6162 "arm={:?}, constraint_block={}x{}, owner_terms={:?}, local_columns={}",
6163 gauge.arm,
6164 gauge.constraint_block.nrows(),
6165 gauge.constraint_block.ncols(),
6166 gauge.owner_terms,
6167 gauge.local_columns,
6168 ),
6169 None => "none".to_string(),
6170 };
6171 let collection_gauge_local_columns = collection_gauge
6172 .as_ref()
6173 .map(|gauge| gauge.local_columns);
6174 let (
6175 design_local,
6176 metadata,
6177 active_penalties,
6178 dropped_penalties,
6179 linear_constraints_local,
6180 joint_null_rotation,
6181 regauged_residualization,
6182 ) = match collection_gauge {
6183 Some(gauge) => {
6184 let placed = gam_terms::smooth::place_term_in_collection_gauge(
6185 &gauge,
6186 gam_terms::smooth::LocalTermRealization {
6187 design: design_local,
6188 metadata: &metadata,
6189 active_penalties: &active_penalties,
6190 dropped_penalties,
6191 linear_constraints_local: linear_constraints_local.as_ref(),
6192 joint_null_rotation: joint_null_rotation.as_ref(),
6193 termname: &name,
6194 },
6195 )
6196 .map_err(|e| {
6197 EstimationError::InvalidInput(format!(
6198 "term '{name}' could not be returned to its collection's identifiability \
6199 gauge after an incremental rebuild: {e}"
6200 ))
6201 })?;
6202 (
6203 placed.design,
6204 placed.metadata,
6205 placed.active_penalties,
6206 placed.dropped_penalties,
6207 placed.linear_constraints_local,
6208 None,
6211 Some(placed.parametric_residualization),
6212 )
6213 }
6214 None => (
6215 design_local,
6216 metadata,
6217 active_penalties,
6218 dropped_penalties,
6219 linear_constraints_local,
6220 joint_null_rotation,
6221 None,
6222 ),
6223 };
6224 let dropped_penaltyinfo: Vec<DroppedPenaltyBlockInfo> = dropped_penalties
6228 .iter()
6229 .map(|info| DroppedPenaltyBlockInfo {
6230 termname: Some(name.clone()),
6231 penalty: info.clone(),
6232 })
6233 .collect();
6234 let coeff_range = self
6235 .design
6236 .smooth
6237 .terms
6238 .get(term_idx)
6239 .ok_or_else(|| EstimationError::InvalidInput(format!("incremental realizer smooth term {term_idx} out of range")))?
6240 .coeff_range
6241 .clone();
6242 if design_local.ncols() != coeff_range.len() {
6243 let local_width_moved = collection_gauge_local_columns
6259 .is_some_and(|expected| pre_gauge_cols != expected);
6260 let reason = format!(
6261 "incremental realizer width mismatch for term {term_idx} ('{name}'): rebuilt_cols={}, \
6262 cached_cols={}; the local rebuild produced {pre_gauge_cols} column(s) before the \
6263 collection gauge ({gauge_report}) and {} after it. Trial: {trial_report}",
6264 design_local.ncols(),
6265 coeff_range.len(),
6266 design_local.ncols(),
6267 );
6268 if local_width_moved {
6269 return Err(EstimationError::InvalidInput(format!(
6270 "{reason}. The LOCAL width moved, so this rebuild is not the basis the \
6271 collection gauged (gam#2760)"
6272 )));
6273 }
6274 return Err(EstimationError::TrialPointRefused {
6275 reason: format!(
6276 "{reason}. The local width is unchanged, so the realized design loses rank in \
6277 the collection gauge's chart at this psi and the model the collection \
6278 specified does not exist here (gam#2760)"
6279 ),
6280 });
6281 }
6282 if design_local.nrows() != self.design.design.nrows() {
6283 return Err(EstimationError::InvalidInput(SmoothError::dimension_mismatch(format!(
6284 "incremental realizer row mismatch for term {}: rebuilt_rows={}, design_rows={}",
6285 term_idx,
6286 design_local.nrows(),
6287 self.design.design.nrows()
6288 )).to_string()));
6289 }
6290
6291 let smooth_penalty_range = self
6292 .smooth_penalty_ranges
6293 .get(term_idx)
6294 .ok_or_else(|| {
6295 EstimationError::InvalidInput(format!("incremental realizer missing smooth penalty range for term {term_idx}"))
6296 })?
6297 .clone();
6298 let full_penalty_range = self
6299 .full_penalty_ranges
6300 .get(term_idx)
6301 .ok_or_else(|| {
6302 EstimationError::InvalidInput(format!("incremental realizer missing full penalty range for term {term_idx}"))
6303 })?
6304 .clone();
6305 let cached_originals: Vec<usize> = self
6322 .design
6323 .smooth
6324 .terms
6325 .get(term_idx)
6326 .map(|term| {
6327 term.active_penalties
6328 .iter()
6329 .map(|active| active.info.original_index)
6330 .collect()
6331 })
6332 .unwrap_or_default();
6333 let (active_penalties, dropped_penalties) = if cached_originals.len()
6334 == smooth_penalty_range.len()
6335 && active_penalties.len() != smooth_penalty_range.len()
6336 {
6337 let mut slots: Vec<Option<gam_terms::basis::ActivePenalty>> =
6338 active_penalties.into_iter().map(Some).collect();
6339 let mut kept = Vec::with_capacity(cached_originals.len());
6340 for original in &cached_originals {
6341 let Some(found) = slots
6342 .iter_mut()
6343 .find(|slot| {
6344 slot.as_ref()
6345 .is_some_and(|active| active.info.original_index == *original)
6346 })
6347 .and_then(Option::take)
6348 else {
6349 return Err(EstimationError::InvalidInput(SmoothError::dimension_mismatch(format!(
6350 "incremental realizer lost cached penalty {original} for term \
6351 '{name}': the rebuild produced {:?}",
6352 slots
6353 .iter()
6354 .flatten()
6355 .map(|active| active.info.original_index)
6356 .collect::<Vec<_>>()
6357 )).to_string()));
6358 };
6359 kept.push(found);
6360 }
6361 let mut dropped = dropped_penalties;
6362 dropped.extend(slots.into_iter().flatten().map(|active| {
6363 gam_terms::basis::DroppedPenaltyInfo {
6364 source: active.info.source.clone(),
6365 original_index: active.info.original_index,
6366 reason: gam_terms::basis::PenaltyDropReason::ZeroMatrix,
6367 normalization_scale: active.info.normalization_scale,
6368 }
6369 }));
6370 (kept, dropped)
6371 } else {
6372 (active_penalties, dropped_penalties)
6373 };
6374 if active_penalties.len() != smooth_penalty_range.len() {
6375 return Err(EstimationError::InvalidInput(SmoothError::dimension_mismatch(format!(
6376 "incremental realizer topology changed for term '{}': active_penalties={}, cached_penalties={}",
6377 name,
6378 active_penalties.len(),
6379 smooth_penalty_range.len()
6380 )).to_string()));
6381 }
6382
6383 self.design.smooth.term_designs[term_idx] = design_local;
6384
6385 for (offset, active_penalty) in active_penalties.iter().enumerate() {
6386 let smooth_penalty_idx = smooth_penalty_range.start + offset;
6387 let full_penalty_idx = full_penalty_range.start + offset;
6388 let penalty_local = &active_penalty.matrix;
6389
6390 if penalty_local.nrows() != coeff_range.len()
6391 || penalty_local.ncols() != coeff_range.len()
6392 {
6393 return Err(EstimationError::InvalidInput(
6394 SmoothError::dimension_mismatch(format!(
6395 "incremental realizer penalty shape mismatch for term '{}' penalty {}: \
6396 penalty is {}x{} but coeff_range has {} columns",
6397 name,
6398 offset,
6399 penalty_local.nrows(),
6400 penalty_local.ncols(),
6401 coeff_range.len()
6402 ))
6403 .to_string(),
6404 ));
6405 }
6406
6407 let smooth_penalty = self
6408 .design
6409 .smooth
6410 .penalties
6411 .get_mut(smooth_penalty_idx)
6412 .ok_or_else(|| {
6413 EstimationError::InvalidInput(format!(
6414 "incremental realizer smooth penalty {} out of range for term {}",
6415 smooth_penalty_idx, term_idx
6416 ))
6417 })?;
6418 smooth_penalty.local.assign(penalty_local);
6421 smooth_penalty.op = active_penalty.op.clone();
6422
6423 let full_bp = self
6424 .design
6425 .penalties
6426 .get_mut(full_penalty_idx)
6427 .ok_or_else(|| {
6428 EstimationError::InvalidInput(format!(
6429 "incremental realizer full penalty {} out of range for term {}",
6430 full_penalty_idx, term_idx
6431 ))
6432 })?;
6433 full_bp.local.assign(penalty_local);
6436 full_bp.op = active_penalty.op.clone();
6437
6438 self.design.smooth.nullspace_dims[smooth_penalty_idx] = active_penalty.nullity;
6439 self.design.nullspace_dims[full_penalty_idx] = active_penalty.nullity;
6440
6441 self.design.smooth.penaltyinfo[smooth_penalty_idx].global_index = smooth_penalty_idx;
6442 self.design.smooth.penaltyinfo[smooth_penalty_idx].termname = Some(name.clone());
6443 self.design.smooth.penaltyinfo[smooth_penalty_idx].penalty =
6444 active_penalty.info.clone();
6445
6446 self.design.penaltyinfo[full_penalty_idx].global_index = full_penalty_idx;
6447 self.design.penaltyinfo[full_penalty_idx].termname = Some(name.clone());
6448 self.design.penaltyinfo[full_penalty_idx].penalty = active_penalty.info.clone();
6449 }
6450
6451 let target_term = self.design.smooth.terms.get_mut(term_idx).ok_or_else(|| {
6452 EstimationError::InvalidInput(format!("incremental realizer smooth term {term_idx} disappeared during replacement"))
6453 })?;
6454 target_term.active_penalties = active_penalties;
6455 target_term.dropped_penalties = dropped_penalties;
6456 target_term.metadata = metadata;
6457 target_term.lower_bounds_local = lower_bounds_local;
6458 target_term.linear_constraints_local = linear_constraints_local;
6459 target_term.joint_null_rotation = joint_null_rotation;
6460 if let Some(chart) = regauged_residualization {
6464 target_term.parametric_residualization = chart;
6465 }
6466 self.dropped_penaltyinfo_by_term[term_idx] = dropped_penaltyinfo;
6467 log::info!(
6468 "[STAGE] smooth basis rebuild (term {}, '{}', cols={}): {:.3}s",
6469 term_idx,
6470 target_term.name,
6471 coeff_range.len(),
6472 t_replace.elapsed().as_secs_f64(),
6473 );
6474 Ok(())
6475 }
6476
6477 fn refresh_full_design_operator(&mut self) -> Result<(), String> {
6478 let mut blocks = Vec::<DesignBlock>::with_capacity(
6479 self.fixed_blocks.len() + self.design.smooth.term_designs.len(),
6480 );
6481 blocks.extend(self.fixed_blocks.iter().cloned());
6482 for term_design in &self.design.smooth.term_designs {
6483 blocks.push(DesignBlock::from(term_design));
6484 }
6485 self.design.design = assemble_term_collection_design_matrix(blocks)
6486 .map_err(|e| format!("failed to refresh term-collection design: {e}"))?;
6487 Ok(())
6488 }
6489}
6490
6491fn build_term_collection_fixed_blocks(
6492 data: ArrayView2<'_, f64>,
6493 spec: &TermCollectionSpec,
6494) -> Result<Vec<DesignBlock>, BasisError> {
6495 let mut blocks = Vec::<DesignBlock>::new();
6496 if !term_collection_has_anchored_bspline(spec) {
6497 blocks.push(DesignBlock::Intercept(data.nrows()));
6498 }
6499
6500 if !spec.linear_terms.is_empty() {
6501 let mut linear_block = Array2::<f64>::zeros((data.nrows(), spec.linear_terms.len()));
6502 for (j, linear) in spec.linear_terms.iter().enumerate() {
6503 let column = linear
6507 .realized_design_column(data)
6508 .map_err(BasisError::InvalidInput)?;
6509 linear_block.column_mut(j).assign(&column);
6510 }
6511 blocks.push(DesignBlock::Dense(
6512 gam_linalg::matrix::DenseDesignMatrix::from(linear_block),
6513 ));
6514 }
6515
6516 for term in &spec.random_effect_terms {
6517 let block = build_random_effect_block(data, term)?;
6518 let re_op = RandomEffectOperator::new(block.group_ids, block.num_groups);
6519 blocks.push(DesignBlock::RandomEffect(Arc::new(re_op)));
6520 }
6521
6522 Ok(blocks)
6523}
6524
6525pub struct SpatialLengthScaleOptimizationResult<FitOut> {
6530 pub resolved_specs: Vec<TermCollectionSpec>,
6531 pub designs: Vec<TermCollectionDesign>,
6532 pub fit: FitOut,
6533 pub certified_outer: Option<gam_solve::rho_optimizer::CertifiedOuterResult>,
6534 pub timing: Option<SpatialLengthScaleOptimizationTiming>,
6535}
6536
6537pub struct ExactJointEvaluation<M> {
6545 pub objective: f64,
6546 pub gradient: Array1<f64>,
6547 pub hessian: gam_problem::HessianValue,
6548 pub mode: M,
6549}
6550
6551pub struct ExactJointEfsEvaluation<M> {
6554 pub evaluation: gam_problem::EfsEval,
6555 pub mode: M,
6556}
6557
6558pub enum SpatialFitProvenance<'a, M> {
6559 NoOuterOptimization,
6560 Certified {
6561 outer: &'a gam_solve::rho_optimizer::CertifiedOuterResult,
6562 mode: M,
6563 },
6564}
6565
6566#[derive(Debug, Clone)]
6568pub struct ExactJointHyperSetup {
6569 rho0: Array1<f64>,
6570 rho_lower: Array1<f64>,
6571 rho_upper: Array1<f64>,
6572 log_kappa0: SpatialLogKappaCoords,
6573 log_kappa_lower: SpatialLogKappaCoords,
6574 log_kappa_upper: SpatialLogKappaCoords,
6575 auxiliary0: Array1<f64>,
6576 auxiliary_lower: Array1<f64>,
6577 auxiliary_upper: Array1<f64>,
6578}
6579
6580impl ExactJointHyperSetup {
6581 fn sanitize_rho_seed(
6582 rho0: Array1<f64>,
6583 rho_lower: &Array1<f64>,
6584 rho_upper: &Array1<f64>,
6585 ) -> Array1<f64> {
6586 Array1::from_iter(rho0.iter().enumerate().map(|(idx, &value)| {
6587 let lo = rho_lower[idx];
6588 let hi = rho_upper[idx];
6589 let fallback = 0.0_f64.clamp(lo, hi);
6590 if value.is_finite() {
6591 value.clamp(lo, hi)
6592 } else {
6593 fallback
6594 }
6595 }))
6596 }
6597
6598 pub(crate) fn new(
6599 rho0: Array1<f64>,
6600 rho_lower: Array1<f64>,
6601 rho_upper: Array1<f64>,
6602 log_kappa0: SpatialLogKappaCoords,
6603 log_kappa_lower: SpatialLogKappaCoords,
6604 log_kappa_upper: SpatialLogKappaCoords,
6605 ) -> Self {
6606 let rho0 = Self::sanitize_rho_seed(rho0, &rho_lower, &rho_upper);
6607 Self {
6608 rho0,
6609 rho_lower,
6610 rho_upper,
6611 log_kappa0,
6612 log_kappa_lower,
6613 log_kappa_upper,
6614 auxiliary0: Array1::zeros(0),
6615 auxiliary_lower: Array1::zeros(0),
6616 auxiliary_upper: Array1::zeros(0),
6617 }
6618 }
6619
6620 pub(crate) fn with_auxiliary(
6621 mut self,
6622 auxiliary0: Array1<f64>,
6623 auxiliary_lower: Array1<f64>,
6624 auxiliary_upper: Array1<f64>,
6625 ) -> Self {
6626 assert_eq!(
6627 auxiliary0.len(),
6628 auxiliary_lower.len(),
6629 "auxiliary lower bound length mismatch"
6630 );
6631 assert_eq!(
6632 auxiliary0.len(),
6633 auxiliary_upper.len(),
6634 "auxiliary upper bound length mismatch"
6635 );
6636 self.auxiliary0 = Self::sanitize_rho_seed(auxiliary0, &auxiliary_lower, &auxiliary_upper);
6637 self.auxiliary_lower = auxiliary_lower;
6638 self.auxiliary_upper = auxiliary_upper;
6639 self
6640 }
6641
6642 pub(crate) fn rho_dim(&self) -> usize {
6643 self.rho0.len()
6644 }
6645
6646 pub(crate) fn log_kappa_dim(&self) -> usize {
6647 self.log_kappa0.len()
6648 }
6649
6650 pub(crate) fn auxiliary_dim(&self) -> usize {
6651 self.auxiliary0.len()
6652 }
6653
6654 pub(crate) fn theta0(&self) -> Array1<f64> {
6655 let mut out =
6656 Array1::<f64>::zeros(self.rho_dim() + self.log_kappa_dim() + self.auxiliary_dim());
6657 out.slice_mut(s![..self.rho_dim()]).assign(&self.rho0);
6658 out.slice_mut(s![self.rho_dim()..self.rho_dim() + self.log_kappa_dim()])
6659 .assign(self.log_kappa0.as_array());
6660 out.slice_mut(s![self.rho_dim() + self.log_kappa_dim()..])
6661 .assign(&self.auxiliary0);
6662 out
6663 }
6664
6665 pub(crate) fn lower(&self) -> Array1<f64> {
6666 let mut out =
6667 Array1::<f64>::zeros(self.rho_dim() + self.log_kappa_dim() + self.auxiliary_dim());
6668 out.slice_mut(s![..self.rho_dim()]).assign(&self.rho_lower);
6669 out.slice_mut(s![self.rho_dim()..self.rho_dim() + self.log_kappa_dim()])
6670 .assign(self.log_kappa_lower.as_array());
6671 out.slice_mut(s![self.rho_dim() + self.log_kappa_dim()..])
6672 .assign(&self.auxiliary_lower);
6673 out
6674 }
6675
6676 pub(crate) fn upper(&self) -> Array1<f64> {
6677 let mut out =
6678 Array1::<f64>::zeros(self.rho_dim() + self.log_kappa_dim() + self.auxiliary_dim());
6679 out.slice_mut(s![..self.rho_dim()]).assign(&self.rho_upper);
6680 out.slice_mut(s![self.rho_dim()..self.rho_dim() + self.log_kappa_dim()])
6681 .assign(self.log_kappa_upper.as_array());
6682 out.slice_mut(s![self.rho_dim() + self.log_kappa_dim()..])
6683 .assign(&self.auxiliary_upper);
6684 out
6685 }
6686
6687 pub(crate) fn log_kappa_dims_per_term(&self) -> Vec<usize> {
6689 self.log_kappa0.dims_per_term().to_vec()
6690 }
6691}
6692
6693struct ExactJointDesignCache<'d> {
6699 realizers: Vec<FrozenTermCollectionIncrementalRealizer<'d>>,
6700 block_term_indices: Vec<Vec<usize>>,
6701 current_theta: Option<Array1<f64>>,
6702 last_cost: Option<f64>,
6703 last_eval: Option<(f64, Array1<f64>, gam_problem::HessianValue)>,
6704 rho_dim: usize,
6705 all_dims: Vec<usize>,
6706 log_kappa_dim: usize,
6707 block_term_counts: Vec<usize>,
6708}
6709
6710impl<'d> ExactJointDesignCache<'d> {
6711 fn new(
6712 data: ArrayView2<'d, f64>,
6713 blocks: Vec<(TermCollectionSpec, TermCollectionDesign, Vec<usize>)>,
6714 rho_dim: usize,
6715 all_dims: Vec<usize>,
6716 ) -> Result<Self, String> {
6717 let n_blocks = blocks.len();
6718 let mut realizers = Vec::with_capacity(n_blocks);
6719 let mut block_term_indices = Vec::with_capacity(n_blocks);
6720 let mut block_term_counts = Vec::with_capacity(n_blocks);
6721
6722 for (spec, design, terms) in blocks {
6723 block_term_counts.push(terms.len());
6724 block_term_indices.push(terms);
6725 realizers.push(FrozenTermCollectionIncrementalRealizer::new(
6726 data, spec, design,
6727 )?);
6728 }
6729
6730 Ok(Self {
6731 realizers,
6732 block_term_indices,
6733 current_theta: None,
6734 last_cost: None,
6735 last_eval: None,
6736 rho_dim,
6737 log_kappa_dim: all_dims.iter().sum(),
6738 all_dims,
6739 block_term_counts,
6740 })
6741 }
6742
6743 fn ensure_theta(&mut self, theta: &Array1<f64>) -> Result<(), String> {
6744 if self
6745 .current_theta
6746 .as_ref()
6747 .is_some_and(|cached| theta_values_match(cached, theta))
6748 {
6749 return Ok(());
6750 }
6751
6752 let t_ensure = std::time::Instant::now();
6753 let kappa_theta_len = self.rho_dim + self.log_kappa_dim;
6754 if theta.len() < kappa_theta_len {
6755 return Err(SmoothError::dimension_mismatch(format!(
6756 "exact-joint theta length mismatch: got {}, expected at least {} (rho_dim={}, log_kappa_dim={})",
6757 theta.len(),
6758 kappa_theta_len,
6759 self.rho_dim,
6760 self.log_kappa_dim
6761 ))
6762 .into());
6763 }
6764 let theta_kappa = theta.slice(s![..kappa_theta_len]).to_owned();
6765 let full_log_kappa = SpatialLogKappaCoords::from_theta_tail_with_dims(
6766 &theta_kappa,
6767 self.rho_dim,
6768 self.all_dims.clone(),
6769 );
6770
6771 let n = self.realizers.len();
6775 let mut remaining = full_log_kappa;
6776 for block_idx in 0..n {
6777 let count = self.block_term_counts[block_idx];
6778 if block_idx < n - 1 {
6779 let (block_lk, rest) = remaining.split_at(count);
6780 self.realizers[block_idx]
6781 .apply_log_kappa(&block_lk, &self.block_term_indices[block_idx])
6782 .map_err(|e| e.to_string())?;
6783 remaining = rest;
6784 } else {
6785 self.realizers[block_idx]
6787 .apply_log_kappa(&remaining, &self.block_term_indices[block_idx])
6788 .map_err(|e| e.to_string())?;
6789 }
6790 }
6791
6792 log::info!(
6793 "[STAGE] ensure_theta (n-block, {} blocks, {} realizers): {:.3}s",
6794 n,
6795 self.realizers.len(),
6796 t_ensure.elapsed().as_secs_f64(),
6797 );
6798 self.current_theta = Some(theta.clone());
6799 self.last_cost = None;
6800 self.last_eval = None;
6801 Ok(())
6802 }
6803
6804 impl_exact_joint_theta_memo!();
6805
6806 fn store_cost_only(&mut self, theta: &Array1<f64>, cost: f64) {
6812 if self
6813 .current_theta
6814 .as_ref()
6815 .is_some_and(|cached| theta_values_match(cached, theta))
6816 {
6817 self.last_cost = Some(cost);
6818 }
6819 }
6820
6821 fn invalidate_objective_memo(&mut self) {
6824 self.last_cost = None;
6825 self.last_eval = None;
6826 }
6827
6828 fn specs(&self) -> Vec<&TermCollectionSpec> {
6829 self.realizers.iter().map(|r| r.spec()).collect()
6830 }
6831
6832 fn designs(&self) -> Vec<&TermCollectionDesign> {
6833 self.realizers.iter().map(|r| r.design()).collect()
6834 }
6835
6836 fn design_revision(&self) -> u64 {
6846 self.realizers
6847 .iter()
6848 .fold(0u64, |acc, r| acc.wrapping_add(r.design_revision()))
6849 }
6850}
6851
6852pub(crate) fn seed_risk_profile_for_likelihood_family(
6853 family: &LikelihoodSpec,
6854) -> gam_problem::SeedRiskProfile {
6855 match &family.response {
6856 ResponseFamily::Gaussian => gam_problem::SeedRiskProfile::Gaussian,
6857 ResponseFamily::RoystonParmar => gam_problem::SeedRiskProfile::Survival,
6858 ResponseFamily::Binomial
6859 | ResponseFamily::Poisson
6860 | ResponseFamily::Tweedie { .. }
6861 | ResponseFamily::NegativeBinomial { .. }
6862 | ResponseFamily::Beta { .. }
6863 | ResponseFamily::Gamma => gam_problem::SeedRiskProfile::GeneralizedLinear,
6864 }
6865}
6866
6867fn exact_joint_seed_config(
6868 risk_profile: gam_problem::SeedRiskProfile,
6869 auxiliary_dim: usize,
6870 initial_seed_only: bool,
6871) -> gam_problem::SeedConfig {
6872 let mut config = gam_problem::SeedConfig {
6873 risk_profile,
6874 num_auxiliary_trailing: auxiliary_dim,
6875 ..Default::default()
6876 };
6877 match risk_profile {
6878 gam_problem::SeedRiskProfile::Gaussian
6879 | gam_problem::SeedRiskProfile::GaussianLocationScale => {
6880 config.max_seeds = 4;
6881 config.seed_budget = 2;
6882 }
6883 gam_problem::SeedRiskProfile::GeneralizedLinear => {
6884 config.max_seeds = 1;
6889 config.seed_budget = 1;
6890 config.screen_max_inner_iterations = 8;
6891 }
6892 gam_problem::SeedRiskProfile::Survival => {
6893 config.max_seeds = 8;
6899 config.seed_budget = 4;
6900 config.screen_max_inner_iterations = 8;
6901 }
6902 }
6903 if initial_seed_only {
6904 config.max_seeds = 1;
6911 config.seed_budget = 1;
6912 config.over_smoothing_probe_rho = None;
6913 }
6914 config
6915}
6916
6917#[cfg(test)]
6918mod exact_joint_seed_config_tests {
6919 use super::*;
6920
6921 #[test]
6922 fn exact_joint_marginal_slope_profiles_get_deeper_startup_validation() {
6923 let bms =
6924 exact_joint_seed_config(gam_problem::SeedRiskProfile::GeneralizedLinear, 2, false);
6925 assert_eq!(bms.max_seeds, 1);
6926 assert_eq!(bms.seed_budget, 1);
6927 assert_eq!(bms.screen_max_inner_iterations, 8);
6928 assert_eq!(bms.num_auxiliary_trailing, 2);
6929
6930 let survival = exact_joint_seed_config(gam_problem::SeedRiskProfile::Survival, 3, false);
6931 assert_eq!(survival.max_seeds, 8);
6932 assert_eq!(survival.seed_budget, 4);
6933 assert_eq!(survival.screen_max_inner_iterations, 8);
6934 assert_eq!(survival.num_auxiliary_trailing, 3);
6935 }
6936
6937 #[test]
6938 fn exact_joint_gaussian_keeps_tight_historical_multistart_budget() {
6939 let gaussian = exact_joint_seed_config(gam_problem::SeedRiskProfile::Gaussian, 1, false);
6940 assert_eq!(gaussian.max_seeds, 4);
6941 assert_eq!(gaussian.seed_budget, 2);
6942 assert_eq!(
6943 gaussian.screen_max_inner_iterations,
6944 gam_problem::SeedConfig::default().screen_max_inner_iterations
6945 );
6946 assert_eq!(gaussian.num_auxiliary_trailing, 1);
6947 }
6948
6949 #[test]
6950 fn certified_matern_basin_owns_the_only_joint_start() {
6951 let gaussian = exact_joint_seed_config(gam_problem::SeedRiskProfile::Gaussian, 1, true);
6952 assert_eq!(gaussian.max_seeds, 1);
6953 assert_eq!(gaussian.seed_budget, 1);
6954 assert_eq!(gaussian.over_smoothing_probe_rho, None);
6955 assert_eq!(gaussian.num_auxiliary_trailing, 1);
6956 }
6957}
6958
6959#[cfg(test)]
6964mod joint_rho_search_box_tests {
6965 use super::*;
6966 use gam_solve::estimate::RHO_BOUND;
6967
6968 #[test]
6972 fn every_finite_incumbent_is_strictly_inside_the_box() {
6973 let seeds = Array1::from(vec![
6976 0.0,
6977 -11.9,
6978 -JOINT_RHO_BOUND,
6979 -12.347_446_785_500_143,
6980 -24.126_016_487_917_27,
6981 11.9,
6982 JOINT_RHO_BOUND,
6983 17.5,
6984 ]);
6985 let (lower, upper) = joint_rho_search_box(seeds.view(), JOINT_RHO_BOUND);
6986 for (k, &seed) in seeds.iter().enumerate() {
6987 assert!(
6988 lower[k] < seed && seed < upper[k],
6989 "coordinate {k}: incumbent {seed} is not STRICTLY inside its joint box \
6990 [{}, {}] — it starts the joint search as an active constraint, which is \
6991 exactly the #2760 defect (the pre-fix rule returned lower = seed here)",
6992 lower[k],
6993 upper[k],
6994 );
6995 }
6996 }
6997
6998 #[test]
7001 fn a_strictly_interior_incumbent_keeps_the_historical_box() {
7002 let seeds = Array1::from(vec![0.0, -11.999, 11.999, -3.0, 5.0]);
7003 let (lower, upper) = joint_rho_search_box(seeds.view(), JOINT_RHO_BOUND);
7004 for k in 0..seeds.len() {
7005 assert_eq!(lower[k], -JOINT_RHO_BOUND);
7006 assert_eq!(upper[k], JOINT_RHO_BOUND);
7007 }
7008 }
7009
7010 #[test]
7013 fn the_fallback_is_per_coordinate() {
7014 let seeds = Array1::from(vec![-30.0, 0.0, 20.0]);
7015 let (lower, upper) = joint_rho_search_box(seeds.view(), JOINT_RHO_BOUND);
7016 assert_eq!((lower[0], upper[0]), (-RHO_BOUND, JOINT_RHO_BOUND));
7017 assert_eq!((lower[1], upper[1]), (-JOINT_RHO_BOUND, JOINT_RHO_BOUND));
7018 assert_eq!((lower[2], upper[2]), (-JOINT_RHO_BOUND, RHO_BOUND));
7019 }
7020
7021 #[test]
7026 fn composes_with_the_constant_curvature_upper_widening() {
7027 let seeds = Array1::from(vec![-13.0, 25.0]);
7028 let (lower, upper) = joint_rho_search_box(seeds.view(), RHO_BOUND);
7029 assert_eq!((lower[0], upper[0]), (-RHO_BOUND, RHO_BOUND));
7030 assert_eq!((lower[1], upper[1]), (-JOINT_RHO_BOUND, RHO_BOUND));
7031 }
7032
7033 #[test]
7037 fn the_box_is_always_a_nonempty_subinterval_of_the_engine_rail() {
7038 let seeds = Array1::from(vec![
7039 f64::NEG_INFINITY,
7040 f64::INFINITY,
7041 f64::NAN,
7042 -RHO_BOUND,
7043 RHO_BOUND,
7044 0.0,
7045 ]);
7046 for &upper_bound in &[JOINT_RHO_BOUND, RHO_BOUND] {
7047 let (lower, upper) = joint_rho_search_box(seeds.view(), upper_bound);
7048 for k in 0..seeds.len() {
7049 assert!(lower[k] < upper[k], "coordinate {k} has an empty box");
7050 assert!(lower[k] >= -RHO_BOUND, "coordinate {k} escaped the engine rail");
7051 assert!(upper[k] <= RHO_BOUND, "coordinate {k} escaped the engine rail");
7052 }
7053 }
7054 }
7055
7056 #[test]
7061 fn a_nonfinite_incumbent_keeps_the_prior() {
7062 let seeds = Array1::from(vec![f64::NEG_INFINITY, f64::INFINITY, f64::NAN]);
7063 let (lower, upper) = joint_rho_search_box(seeds.view(), JOINT_RHO_BOUND);
7064 for k in 0..seeds.len() {
7065 assert_eq!(lower[k], -JOINT_RHO_BOUND);
7066 assert_eq!(upper[k], JOINT_RHO_BOUND);
7067 }
7068 }
7069}
7070
7071pub(crate) fn exact_joint_multistart_outer_problem(
7072 theta0: &Array1<f64>,
7073 lower: &Array1<f64>,
7074 upper: &Array1<f64>,
7075 rho_dim: usize,
7076 auxiliary_dim: usize,
7077 n_params: usize,
7078 gradient: gam_problem::Derivative,
7079 hessian: gam_problem::DeclaredHessianForm,
7080 disable_fixed_point: bool,
7081 risk_profile: gam_problem::SeedRiskProfile,
7082 tolerance: f64,
7083 max_iter: usize,
7084 bfgs_step_cap: Option<f64>,
7093 bfgs_step_cap_psi: Option<f64>,
7094 screening_cap: Option<Arc<AtomicUsize>>,
7095 profiled_objective_size: Option<(usize, usize)>,
7116 has_constant_curvature: bool,
7125 initial_seed_only: bool,
7130) -> Result<gam_solve::rho_optimizer::OuterProblem, EstimationError> {
7131 if rho_dim > theta0.len() {
7132 crate::bail_invalid_estim!(
7133 "exact joint outer problem declares {rho_dim} smoothing coordinates for theta length {}",
7134 theta0.len(),
7135 );
7136 }
7137 let mut seed_heuristic = theta0.to_vec();
7138 let initial_lambdas = gam_problem::checked_exp_log_strengths(
7139 theta0.iter().take(rho_dim).copied(),
7140 )
7141 .map_err(|error| {
7142 EstimationError::InvalidInput(format!(
7143 "exact joint initial smoothing coordinate is outside the canonical log-strength domain: {error}"
7144 ))
7145 })?;
7146 for (value, lambda) in seed_heuristic[..rho_dim].iter_mut().zip(initial_lambdas) {
7147 *value = lambda;
7148 }
7149 let rho_ceiling = if has_constant_curvature {
7154 gam_solve::estimate::RHO_BOUND
7155 } else {
7156 12.0
7157 };
7158 let mut problem = gam_solve::rho_optimizer::OuterProblem::new(n_params)
7159 .with_gradient(gradient)
7160 .with_hessian(hessian)
7161 .with_prefer_gradient_only(true)
7168 .with_require_measured_psd(!matches!(
7200 hessian,
7201 gam_problem::DeclaredHessianForm::Unavailable
7202 ))
7203 .with_disable_fixed_point(disable_fixed_point)
7204 .with_fallback_policy(gam_solve::rho_optimizer::FallbackPolicy::Automatic)
7214 .with_psi_dim(auxiliary_dim)
7215 .with_tolerance(tolerance)
7216 .with_max_iter(max_iter)
7217 .with_bounds(lower.clone(), upper.clone())
7218 .with_initial_rho(theta0.clone())
7219 .with_bfgs_step_cap(bfgs_step_cap)
7220 .with_bfgs_step_cap_psi(bfgs_step_cap_psi)
7221 .with_seed_config({
7222 let mut sc = exact_joint_seed_config(risk_profile, auxiliary_dim, initial_seed_only);
7223 if has_constant_curvature {
7224 sc.bounds = (sc.bounds.0, rho_ceiling);
7228 }
7243 sc
7244 })
7245 .with_rho_bound(rho_ceiling)
7246 .with_heuristic_lambdas(seed_heuristic);
7247 if let Some((n_obs, p_cols)) = profiled_objective_size {
7248 problem = problem
7253 .with_objective_scale(Some(n_obs as f64))
7254 .with_problem_size(n_obs, p_cols);
7255 }
7256 if let Some(screening_cap) = screening_cap {
7257 problem = problem
7258 .with_screening_cap(screening_cap)
7259 .with_screen_initial_rho(true);
7260 }
7261 Ok(problem)
7262}
7263
7264pub fn optimize_spatial_length_scale_exact_joint<FitOut, Mode, FitFn, ExactFn, ExactEfsFn, SeedFn>(
7265 data: ArrayView2<'_, f64>,
7266 block_specs: &[TermCollectionSpec],
7267 block_term_indices: &[Vec<usize>],
7268 kappa_options: &SpatialLengthScaleOptimizationOptions,
7269 joint_setup: &ExactJointHyperSetup,
7270 seed_risk_profile: gam_problem::SeedRiskProfile,
7271 analytic_joint_gradient_available: bool,
7272 analytic_joint_hessian_available: bool,
7273 disable_fixed_point: bool,
7274 screening_cap: Option<Arc<AtomicUsize>>,
7275 outer_derivative_policy: gam_model_api::families::custom_family::OuterDerivativePolicy,
7276 mut fit_fn: FitFn,
7277 mut exact_fn: ExactFn,
7278 mut exact_efs_fn: ExactEfsFn,
7279 mut seed_inner_beta_fn: SeedFn,
7280) -> Result<SpatialLengthScaleOptimizationResult<FitOut>, String>
7281where
7282 FitFn: FnMut(
7283 &Array1<f64>,
7284 &[TermCollectionSpec],
7285 &[TermCollectionDesign],
7286 SpatialFitProvenance<'_, Mode>,
7287 ) -> Result<FitOut, String>,
7288 ExactFn: FnMut(
7289 &Array1<f64>,
7290 &[TermCollectionSpec],
7291 &[TermCollectionDesign],
7292 gam_solve::estimate::reml::reml_outer_engine::EvalMode,
7293 &gam_problem::outer_subsample::RowSet,
7294 Option<Mode>,
7295 ) -> Result<ExactJointEvaluation<Mode>, String>,
7296 ExactEfsFn: FnMut(
7297 &Array1<f64>,
7298 &[TermCollectionSpec],
7299 &[TermCollectionDesign],
7300 &gam_problem::outer_subsample::RowSet,
7301 ) -> Result<ExactJointEfsEvaluation<Mode>, String>,
7302 SeedFn: FnMut(&Array1<f64>) -> Result<gam_solve::rho_optimizer::SeedOutcome, EstimationError>,
7303{
7304 let n_blocks = block_specs.len();
7305 if block_term_indices.len() != n_blocks {
7306 return Err(SmoothError::dimension_mismatch(format!(
7307 "block_specs ({}) and block_term_indices ({}) length mismatch",
7308 n_blocks,
7309 block_term_indices.len()
7310 ))
7311 .into());
7312 }
7313
7314 let log_kappa_dim = joint_setup.log_kappa_dim();
7315
7316 log::trace!(
7317 "[spatial-exact-joint] driver entry: aux_dim={} log_kappa_dim={} kappa_enabled={} rho_dim={} theta0_len={}",
7318 joint_setup.auxiliary_dim(),
7319 log_kappa_dim,
7320 kappa_options.enabled,
7321 joint_setup.rho_dim(),
7322 joint_setup.theta0().len()
7323 );
7324
7325 if joint_setup.auxiliary_dim() == 0 && (!kappa_options.enabled || log_kappa_dim == 0) {
7329 log::trace!(
7330 "[spatial-exact-joint] taking fast path (no outer theta optimization in this driver)"
7331 );
7332 let (designs, resolved_specs) = build_term_collection_designs_and_freeze_joint(
7333 data, block_specs,
7334 )
7335 .map_err(|e| {
7336 format!("failed to build and freeze joint block designs during exact joint kappa optimization: {e}")
7337 })?;
7338 let theta0 = joint_setup.theta0();
7339
7340 let spec_refs: Vec<TermCollectionSpec> = resolved_specs.clone();
7342 let design_refs: Vec<TermCollectionDesign> = designs.clone();
7343 let fit = fit_fn(
7344 &theta0,
7345 &spec_refs,
7346 &design_refs,
7347 SpatialFitProvenance::NoOuterOptimization,
7348 )?;
7349 return Ok(SpatialLengthScaleOptimizationResult {
7350 resolved_specs,
7351 designs,
7352 fit,
7353 certified_outer: None,
7354 timing: None,
7355 });
7356 }
7357
7358 let theta0 = joint_setup.theta0();
7362 let lower = joint_setup.lower();
7363 let upper = joint_setup.upper();
7364 if theta0.len() < log_kappa_dim || lower.len() != theta0.len() || upper.len() != theta0.len() {
7365 return Err(SmoothError::dimension_mismatch(format!(
7366 "invalid exact joint theta setup: theta0={}, lower={}, upper={}, required_log_kappa_dim={}",
7367 theta0.len(),
7368 lower.len(),
7369 upper.len(),
7370 log_kappa_dim
7371 ))
7372 .into());
7373 }
7374 let rho_dim = joint_setup.rho_dim();
7375 let all_dims = joint_setup.log_kappa_dims_per_term();
7376
7377 let (boot_designs, best_specs) = build_term_collection_designs_and_freeze_joint(
7379 data,
7380 block_specs,
7381 )
7382 .map_err(|e| {
7383 format!(
7384 "failed to build and freeze joint block designs during exact joint kappa bootstrap: {e}"
7385 )
7386 })?;
7387 let policy_hessian_form = outer_derivative_policy.declared_hessian_form();
7397 let analytic_outer_hessian_available = analytic_joint_hessian_available
7398 && matches!(
7399 policy_hessian_form,
7400 gam_problem::DeclaredHessianForm::Either
7401 | gam_problem::DeclaredHessianForm::Dense
7402 | gam_problem::DeclaredHessianForm::Operator { .. }
7403 );
7404 let theta_dim = theta0.len();
7405 let psi_dim = theta_dim - rho_dim;
7406
7407 let cache_blocks: Vec<(TermCollectionSpec, TermCollectionDesign, Vec<usize>)> = best_specs
7409 .iter()
7410 .zip(boot_designs.iter())
7411 .zip(block_term_indices.iter())
7412 .map(|((spec, design), terms)| (spec.clone(), design.clone(), terms.clone()))
7413 .collect();
7414
7415 struct NBlockExactJointState<'d, M> {
7416 cache: ExactJointDesignCache<'d>,
7417 row_set: gam_problem::outer_subsample::RowSet,
7418 staged_pilot_active: bool,
7419 terminal_mode: Option<(Array1<f64>, f64, M)>,
7420 }
7421
7422 impl<M> NBlockExactJointState<'_, M> {
7423 fn ensure_theta(&mut self, theta: &Array1<f64>) -> Result<(), String> {
7424 let theta_changed = !self
7425 .cache
7426 .current_theta
7427 .as_ref()
7428 .is_some_and(|current| theta_values_match(current, theta));
7429 if theta_changed {
7430 self.terminal_mode = None;
7431 }
7432 self.cache.ensure_theta(theta)
7433 }
7434
7435 fn install_terminal_mode(&mut self, theta: &Array1<f64>, objective: f64, mode: M) {
7436 self.terminal_mode = Some((theta.clone(), objective, mode));
7437 }
7438
7439 fn terminal_mode_matches(&self, theta: &Array1<f64>, objective: f64) -> bool {
7440 self.terminal_mode
7441 .as_ref()
7442 .is_some_and(|(mode_theta, mode_objective, _)| {
7443 theta_values_match(mode_theta, theta)
7444 && mode_objective.to_bits() == objective.to_bits()
7445 })
7446 }
7447
7448 fn take_terminal_mode(&mut self, theta: &Array1<f64>) -> Option<M> {
7449 if self
7450 .terminal_mode
7451 .as_ref()
7452 .is_some_and(|(mode_theta, _, _)| theta_values_match(mode_theta, theta))
7453 {
7454 self.terminal_mode.take().map(|(_, _, mode)| mode)
7455 } else {
7456 None
7457 }
7458 }
7459 }
7460
7461 let mut state = NBlockExactJointState {
7462 cache: ExactJointDesignCache::new(data, cache_blocks, rho_dim, all_dims.clone())?,
7463 row_set: gam_problem::outer_subsample::RowSet::All,
7464 staged_pilot_active: false,
7465 terminal_mode: None,
7466 };
7467
7468 const KAPPA_PILOT_K: usize = 5_000;
7496
7497 let n_total = data.nrows();
7498 let use_staged_kappa = outer_derivative_policy.should_use_staged_kappa(n_total);
7499 if use_staged_kappa {
7500 log::info!(
7501 "[KAPPA-STAGED] auto-engaging pilot+exact schedule: n={} pilot_k={}",
7502 n_total,
7503 KAPPA_PILOT_K,
7504 );
7505 }
7506
7507 fn build_uniform_pilot_subsample(
7524 n_total: usize,
7525 k_target: usize,
7526 seed: u64,
7527 ) -> gam_problem::outer_subsample::OuterScoreSubsample {
7528 use gam_problem::outer_subsample::OuterScoreSubsample;
7529 let k = k_target.min(n_total);
7530 if k == 0 || n_total == 0 {
7531 return OuterScoreSubsample::from_uniform_inclusion_mask(Vec::new(), n_total, seed);
7532 }
7533 let mut mask: Vec<usize> = Vec::with_capacity(k);
7537 let mut state = seed.wrapping_add(0x9E3779B97F4A7C15);
7539 let splitmix = |s: &mut u64| -> u64 { gam_linalg::utils::splitmix64(s) };
7540 let mut taken = std::collections::HashSet::with_capacity(k);
7541 for j in (n_total - k)..n_total {
7542 let r = (splitmix(&mut state) % (j as u64 + 1)) as usize;
7543 if !taken.insert(r) {
7544 taken.insert(j);
7545 mask.push(j);
7546 } else {
7547 mask.push(r);
7548 }
7549 }
7550 mask.sort_unstable();
7551 mask.dedup();
7552 OuterScoreSubsample::from_uniform_inclusion_mask(mask, n_total, seed)
7553 }
7554
7555 if use_staged_kappa {
7556 let pilot = build_uniform_pilot_subsample(n_total, KAPPA_PILOT_K, n_total as u64);
7557 state.row_set = gam_problem::outer_subsample::RowSet::Subsample {
7558 rows: std::sync::Arc::clone(&pilot.rows),
7559 n_full: n_total,
7560 };
7561 state.staged_pilot_active = true;
7562 }
7563
7564 let exact_fn_cell = std::cell::RefCell::new(&mut exact_fn);
7565 let exact_efs_fn_cell = std::cell::RefCell::new(&mut exact_efs_fn);
7566
7567 use std::cell::Cell;
7582 let kphase_cost_calls: Cell<usize> = Cell::new(0);
7583 let kphase_cost_total_s: Cell<f64> = Cell::new(0.0);
7584 let kphase_eval_calls: Cell<usize> = Cell::new(0);
7585 let kphase_eval_total_s: Cell<f64> = Cell::new(0.0);
7586 let kphase_efs_calls: Cell<usize> = Cell::new(0);
7587 let kphase_efs_total_s: Cell<f64> = Cell::new(0.0);
7588 let kphase_optim_start = std::time::Instant::now();
7589 let kphase_log_kappa_dim = log_kappa_dim;
7590 let kphase_log_norms = |theta: &Array1<f64>| -> (f64, f64) {
7591 let theta_norm = theta.iter().map(|v| v * v).sum::<f64>().sqrt();
7592 let log_kappa_norm = if kphase_log_kappa_dim > 0 && theta.len() >= kphase_log_kappa_dim {
7593 let start = theta.len() - kphase_log_kappa_dim;
7594 theta.iter().skip(start).map(|v| v * v).sum::<f64>().sqrt()
7595 } else {
7596 0.0
7597 };
7598 (theta_norm, log_kappa_norm)
7599 };
7600
7601 use gam_problem::{DeclaredHessianForm, Derivative, OuterEval};
7602 use gam_solve::rho_optimizer::OuterEvalOrder;
7603
7604 let joint_p_cols: usize = boot_designs
7608 .iter()
7609 .map(|d| d.design.ncols())
7610 .sum::<usize>()
7611 .max(1);
7612
7613 let problem = exact_joint_multistart_outer_problem(
7614 &theta0,
7615 &lower,
7616 &upper,
7617 rho_dim,
7618 psi_dim,
7619 theta_dim,
7620 if analytic_joint_gradient_available {
7621 Derivative::Analytic
7622 } else {
7623 Derivative::Unavailable
7624 },
7625 if analytic_outer_hessian_available {
7626 DeclaredHessianForm::Either
7627 } else {
7628 DeclaredHessianForm::Unavailable
7629 },
7630 disable_fixed_point,
7631 seed_risk_profile,
7632 kappa_options.rel_tol.max(1e-6),
7633 kappa_options.max_outer_iter.max(1),
7634 Some(5.0),
7636 Some(kappa_options.log_step.clamp(0.25, 1.0)),
7638 screening_cap.clone(),
7639 Some((n_total, joint_p_cols)),
7642 block_specs
7645 .iter()
7646 .any(|s| !constant_curvature_term_indices(s).is_empty()),
7647 false,
7650 )
7651 .map_err(|e| e.to_string())?;
7652
7653 fn collect_specs(cache: &ExactJointDesignCache<'_>) -> Vec<TermCollectionSpec> {
7655 cache.specs().into_iter().cloned().collect()
7656 }
7657 fn collect_designs(cache: &ExactJointDesignCache<'_>) -> Vec<TermCollectionDesign> {
7658 cache.designs().into_iter().cloned().collect()
7659 }
7660
7661 let result = {
7662 let eval_outer = |ctx: &mut &mut NBlockExactJointState<'_, Mode>,
7663 theta: &Array1<f64>,
7664 order: OuterEvalOrder|
7665 -> Result<OuterEval, EstimationError> {
7666 if let Some((cost, grad, hess)) = ctx.cache.memoized_eval(theta)
7667 && ctx.terminal_mode_matches(theta, cost)
7668 {
7669 let cached_satisfies_order = match order {
7670 OuterEvalOrder::Value => true,
7671 OuterEvalOrder::ValueAndGradient => grad.len() == theta.len(),
7672 OuterEvalOrder::ValueGradientHessian => {
7673 grad.len() == theta.len() && hess.is_analytic()
7674 }
7675 };
7676 if cached_satisfies_order {
7677 if !cost.is_finite() {
7678 return Ok(OuterEval::infeasible(theta.len()));
7679 }
7680 if grad.iter().any(|v| !v.is_finite()) {
7693 return Ok(OuterEval::infeasible(theta.len()));
7694 }
7695 return Ok(OuterEval {
7696 cost,
7697 gradient: grad,
7698 hessian: hess,
7699 inner_beta_hint: None,
7700 });
7701 }
7702 }
7703 ctx.ensure_theta(theta).map_err(|err| {
7704 EstimationError::InvalidInput(format!(
7705 "n-block exact-joint spatial design realization failed: {err}"
7706 ))
7707 })?;
7708 let design_revision = Some(ctx.cache.design_revision());
7709 let specs = collect_specs(&ctx.cache);
7710 let designs = collect_designs(&ctx.cache);
7711 let clamped = outer_derivative_policy.order_for_evaluation(order);
7719 let value_only = matches!(clamped, OuterEvalOrder::Value);
7720 let need_hessian = matches!(clamped, OuterEvalOrder::ValueGradientHessian)
7721 && analytic_outer_hessian_available;
7722 let eval_mode = if value_only {
7723 gam_solve::estimate::reml::reml_outer_engine::EvalMode::ValueOnly
7724 } else if need_hessian {
7725 gam_solve::estimate::reml::reml_outer_engine::EvalMode::ValueGradientHessian
7726 } else {
7727 gam_solve::estimate::reml::reml_outer_engine::EvalMode::ValueAndGradient
7728 };
7729 let owned_value_mode = if value_only {
7730 None
7731 } else {
7732 ctx.take_terminal_mode(theta)
7733 };
7734 let t0 = std::time::Instant::now();
7735 let result = (*exact_fn_cell.borrow_mut())(
7736 theta,
7737 &specs,
7738 &designs,
7739 eval_mode,
7740 &ctx.row_set,
7741 owned_value_mode,
7742 );
7743 let elapsed_s = t0.elapsed().as_secs_f64();
7744 kphase_eval_calls.set(kphase_eval_calls.get() + 1);
7745 kphase_eval_total_s.set(kphase_eval_total_s.get() + elapsed_s);
7746 let (theta_norm, log_kappa_norm) = kphase_log_norms(theta);
7747 log::info!(
7748 "[KAPPA-PHASE] phase=eval_outer call={} order={:?} design_revision={:?} theta_norm={:.4e} log_kappa_norm={:.4e} elapsed_s={:.4}",
7749 kphase_eval_calls.get(),
7750 order,
7751 design_revision,
7752 theta_norm,
7753 log_kappa_norm,
7754 elapsed_s,
7755 );
7756 match result {
7757 Ok(ExactJointEvaluation {
7758 objective: cost,
7759 gradient: grad,
7760 hessian: hess,
7761 mode,
7762 }) => {
7763 ctx.install_terminal_mode(theta, cost, mode);
7764 if value_only {
7765 ctx.cache.store_cost_only(theta, cost);
7766 } else {
7767 ctx.cache.store_eval((cost, grad.clone(), hess.clone()));
7768 }
7769 if !cost.is_finite() {
7770 return Ok(OuterEval::infeasible(theta.len()));
7771 }
7772 if grad.iter().any(|v| !v.is_finite()) {
7785 return Ok(OuterEval::infeasible(theta.len()));
7786 }
7787 Ok(OuterEval {
7788 cost,
7789 gradient: grad,
7790 hessian: hess,
7791 inner_beta_hint: None,
7792 })
7793 }
7794 Err(err) => Err(EstimationError::TrialPointRefused {
7802 reason: format!("n-block exact-joint spatial evaluation failed: {err}"),
7803 }),
7804 }
7805 };
7806
7807 let obj = problem.build_objective_with_eval_order(
7808 &mut state,
7809 |ctx: &mut &mut NBlockExactJointState<'_, Mode>, theta: &Array1<f64>| {
7810 if let Some(cost) = ctx.cache.memoized_cost(theta)
7811 && ctx.terminal_mode_matches(theta, cost)
7812 {
7813 return Ok(cost);
7814 }
7815 ctx.ensure_theta(theta).map_err(|err| {
7816 EstimationError::InvalidInput(format!(
7817 "n-block exact-joint spatial design realization failed: {err}"
7818 ))
7819 })?;
7820 let design_revision = Some(ctx.cache.design_revision());
7821 let specs = collect_specs(&ctx.cache);
7822 let designs = collect_designs(&ctx.cache);
7823 let t0 = std::time::Instant::now();
7830 let result = (*exact_fn_cell.borrow_mut())(
7831 theta,
7832 &specs,
7833 &designs,
7834 gam_solve::estimate::reml::reml_outer_engine::EvalMode::ValueOnly,
7835 &ctx.row_set,
7836 None,
7837 );
7838 let elapsed_s = t0.elapsed().as_secs_f64();
7839 kphase_cost_calls.set(kphase_cost_calls.get() + 1);
7840 kphase_cost_total_s.set(kphase_cost_total_s.get() + elapsed_s);
7841 let (theta_norm, log_kappa_norm) = kphase_log_norms(theta);
7842 log::info!(
7843 "[KAPPA-PHASE] phase=cost call={} design_revision={:?} theta_norm={:.4e} log_kappa_norm={:.4e} elapsed_s={:.4}",
7844 kphase_cost_calls.get(),
7845 design_revision,
7846 theta_norm,
7847 log_kappa_norm,
7848 elapsed_s,
7849 );
7850 match result {
7851 Ok(ExactJointEvaluation {
7852 objective: cost,
7853 mode,
7854 ..
7855 }) => {
7856 ctx.install_terminal_mode(theta, cost, mode);
7857 ctx.cache.store_cost_only(theta, cost);
7863 Ok(cost)
7864 }
7865 Err(err) => Err(EstimationError::TrialPointRefused {
7866 reason: format!(
7867 "n-block exact-joint spatial cost evaluation failed: {err}"
7868 ),
7869 }),
7870 }
7871 },
7872 |ctx: &mut &mut NBlockExactJointState<'_, Mode>, theta: &Array1<f64>| {
7873 eval_outer(ctx, theta, OuterEvalOrder::ValueAndGradient)
7876 },
7877 |ctx: &mut &mut NBlockExactJointState<'_, Mode>,
7878 theta: &Array1<f64>,
7879 order: OuterEvalOrder| { eval_outer(ctx, theta, order) },
7880 None::<fn(&mut &mut NBlockExactJointState<'_, Mode>)>,
7881 Some(
7882 |ctx: &mut &mut NBlockExactJointState<'_, Mode>, theta: &Array1<f64>| {
7883 ctx
7884 .ensure_theta(theta)
7885 .map_err(EstimationError::InvalidInput)?;
7886 let design_revision = Some(ctx.cache.design_revision());
7887 let specs = collect_specs(&ctx.cache);
7888 let designs = collect_designs(&ctx.cache);
7889 let t0 = std::time::Instant::now();
7890 let eval_result = (*exact_efs_fn_cell.borrow_mut())(
7891 theta,
7892 &specs,
7893 &designs,
7894 &ctx.row_set,
7895 );
7896 let elapsed_s = t0.elapsed().as_secs_f64();
7897 kphase_efs_calls.set(kphase_efs_calls.get() + 1);
7898 kphase_efs_total_s.set(kphase_efs_total_s.get() + elapsed_s);
7899 let (theta_norm, log_kappa_norm) = kphase_log_norms(theta);
7900 log::info!(
7901 "[KAPPA-PHASE] phase=efs call={} design_revision={:?} theta_norm={:.4e} log_kappa_norm={:.4e} elapsed_s={:.4}",
7902 kphase_efs_calls.get(),
7903 design_revision,
7904 theta_norm,
7905 log_kappa_norm,
7906 elapsed_s,
7907 );
7908 let ExactJointEfsEvaluation { evaluation, mode } =
7909 eval_result.map_err(|reason| EstimationError::TrialPointRefused {
7910 reason,
7911 })?;
7912 ctx.cache.invalidate_objective_memo();
7919 ctx.cache.store_cost_only(theta, evaluation.cost);
7920 ctx.install_terminal_mode(theta, evaluation.cost, mode);
7921 Ok(evaluation)
7922 },
7923 ),
7924 );
7925 let mut obj = obj
7926 .with_seed_inner_state(
7927 move |_: &mut &mut NBlockExactJointState<'_, Mode>, beta: &Array1<f64>| {
7928 (seed_inner_beta_fn)(beta)
7929 },
7930 )
7931 .with_exact_polish(|ctx: &mut &mut NBlockExactJointState<'_, Mode>| {
7932 if !ctx.staged_pilot_active {
7933 return false;
7934 }
7935 ctx.cache.invalidate_objective_memo();
7940 ctx.terminal_mode = None;
7941 ctx.row_set = gam_problem::outer_subsample::RowSet::All;
7942 ctx.staged_pilot_active = false;
7943 true
7944 })
7945 .with_terminal_eval_order(if analytic_outer_hessian_available {
7976 OuterEvalOrder::ValueGradientHessian
7977 } else {
7978 OuterEvalOrder::ValueAndGradient
7979 });
7980
7981 problem
7982 .run_certified(&mut obj, "n-block exact-joint spatial")
7983 .map_err(|error| error.to_string())?
7984 }; let kphase_total_s = kphase_optim_start.elapsed().as_secs_f64();
7994 log::info!(
7995 "[KAPPA-PHASE-SUMMARY] log_kappa_dim={} n_cost={} cost_total_s={:.4} n_eval={} eval_total_s={:.4} n_efs={} efs_total_s={:.4} optim_total_s={:.4}",
7996 kphase_log_kappa_dim,
7997 kphase_cost_calls.get(),
7998 kphase_cost_total_s.get(),
7999 kphase_eval_calls.get(),
8000 kphase_eval_total_s.get(),
8001 kphase_efs_calls.get(),
8002 kphase_efs_total_s.get(),
8003 kphase_total_s,
8004 );
8005 let timing = SpatialLengthScaleOptimizationTiming {
8006 log_kappa_dim: kphase_log_kappa_dim,
8007 cost_calls: kphase_cost_calls.get(),
8008 cost_total_s: kphase_cost_total_s.get(),
8009 eval_calls: kphase_eval_calls.get(),
8010 eval_total_s: kphase_eval_total_s.get(),
8011 efs_calls: kphase_efs_calls.get(),
8012 efs_total_s: kphase_efs_total_s.get(),
8013 slow_path_resets: 0,
8014 design_revision_delta: 0,
8015 nfree_skip_row_touches: 0,
8016 nfree_miss_shape: 0,
8017 nfree_miss_value: 0,
8018 nfree_miss_gradient: 0,
8019 nfree_miss_penalty: 0,
8020 nfree_miss_revision: 0,
8021 nfree_miss_second_order: 0,
8022 nfree_miss_other: 0,
8023 exact_polish_ran: false,
8028 polish_slow_path_resets: 0,
8029 polish_nfree_skip_row_touches: 0,
8030 optim_total_s: kphase_total_s,
8031 };
8032
8033 if !matches!(state.row_set, gam_problem::outer_subsample::RowSet::All) {
8034 return Err(
8035 "n-block exact-joint spatial optimization returned before its exact full-data transition"
8036 .to_string(),
8037 );
8038 }
8039 let certified_outer = result;
8040 let theta_star = certified_outer.rho().clone();
8041
8042 state.ensure_theta(&theta_star)?;
8047 let (mode_theta, mode_objective, mode) = state.terminal_mode.take().ok_or_else(|| {
8048 "n-block exact-joint spatial optimization produced a certificate without retaining the owned terminal coefficient mode"
8049 .to_string()
8050 })?;
8051 if !theta_values_match(&mode_theta, &theta_star) {
8052 return Err(
8053 "n-block exact-joint spatial terminal coefficient mode does not bitwise match the certified hyperparameter vector"
8054 .to_string(),
8055 );
8056 }
8057 if mode_objective.to_bits() != certified_outer.final_value().to_bits() {
8058 return Err(format!(
8059 "n-block exact-joint spatial terminal coefficient mode objective does not bitwise match the certified objective: mode={mode_objective:.17e}, certified={:.17e}",
8060 certified_outer.final_value(),
8061 ));
8062 }
8063
8064 let resolved_specs: Vec<TermCollectionSpec> = collect_specs(&state.cache);
8065 let designs: Vec<TermCollectionDesign> = collect_designs(&state.cache);
8066
8067 let fit = fit_fn(
8068 &theta_star,
8069 &resolved_specs,
8070 &designs,
8071 SpatialFitProvenance::Certified {
8072 outer: &certified_outer,
8073 mode,
8074 },
8075 )?;
8076
8077 for spec in &resolved_specs {
8078 log_spatial_aniso_scales(spec);
8079 }
8080
8081 Ok(SpatialLengthScaleOptimizationResult {
8082 resolved_specs,
8083 designs,
8084 fit,
8085 certified_outer: Some(certified_outer),
8086 timing: Some(timing),
8087 })
8088}
8089
8090fn try_exact_joint_latent_coord_optimization(
8091 data: ArrayView2<'_, f64>,
8092 y: ArrayView1<'_, f64>,
8093 weights: ArrayView1<'_, f64>,
8094 offset: ArrayView1<'_, f64>,
8095 resolvedspec: &TermCollectionSpec,
8096 best: &FittedTermCollection,
8097 family: LikelihoodSpec,
8098 options: &FitOptions,
8099 latent: &StandardLatentCoordConfig,
8100) -> Result<FittedTermCollectionWithSpec, EstimationError> {
8101 use gam_problem::{DeclaredHessianForm, Derivative, OuterEval};
8102 use gam_solve::rho_optimizer::OuterEvalOrder;
8103
8104 let rho_dim = best.fit.lambdas.len();
8105 let latent_flat_dim = latent.values.len();
8106 if latent_flat_dim == 0 {
8107 crate::bail_invalid_estim!(
8108 "latent-coordinate optimization requires a non-empty latent block"
8109 );
8110 }
8111 let direct_hypers =
8112 latent_coord_initial_direct_hypers(latent.values.id_mode(), latent.values.latent_dim())?;
8113 let analytic_rho_count = latent
8114 .analytic_penalties
8115 .as_ref()
8116 .map_or(0, |registry| registry.total_rho_count());
8117 let latent_coord_ext_dim = latent_flat_dim + analytic_rho_count + direct_hypers.len();
8118
8119 let mut theta0 = Array1::<f64>::zeros(rho_dim + latent_coord_ext_dim);
8120 theta0
8121 .slice_mut(s![..rho_dim])
8122 .assign(&best.fit.lambdas.mapv(f64::ln));
8123 theta0
8124 .slice_mut(s![rho_dim..rho_dim + latent_flat_dim])
8125 .assign(latent.values.as_flat());
8126 if !direct_hypers.is_empty() {
8127 let direct_start = rho_dim + latent_flat_dim + analytic_rho_count;
8128 theta0
8129 .slice_mut(s![direct_start..direct_start + direct_hypers.len()])
8130 .assign(&direct_hypers);
8131 }
8132
8133 let mut lower = Array1::<f64>::from_elem(theta0.len(), -12.0);
8134 let mut upper = Array1::<f64>::from_elem(theta0.len(), 12.0);
8135 let latent_bound = latent
8136 .values
8137 .as_flat()
8138 .iter()
8139 .fold(1.0_f64, |acc, &v| acc.max(v.abs()))
8140 + 10.0;
8141 for axis in rho_dim..rho_dim + latent_flat_dim {
8142 lower[axis] = -latent_bound;
8143 upper[axis] = latent_bound;
8144 }
8145 if let Some(registry) = latent.analytic_penalties.as_ref() {
8146 let (domain_lower, domain_upper) = registry
8147 .rho_domain_bounds()
8148 .map_err(EstimationError::InvalidInput)?;
8149 let start = rho_dim + latent_flat_dim;
8150 for local in 0..analytic_rho_count {
8151 lower[start + local] = lower[start + local].max(domain_lower[local]);
8152 upper[start + local] = upper[start + local].min(domain_upper[local]);
8153 if lower[start + local] >= upper[start + local] {
8154 return Err(EstimationError::InvalidInput(format!(
8155 "analytic-penalty rho domain has no searchable interval at coordinate {local}: lower={}, upper={}",
8156 lower[start + local],
8157 upper[start + local]
8158 )));
8159 }
8160 }
8161 }
8162
8163 struct LatentJointContext<'d> {
8164 rho_dim: usize,
8165 cache: SingleBlockLatentCoordDesignCache,
8166 evaluator: gam_solve::estimate::ExternalJointHyperEvaluator<'d>,
8167 }
8168
8169 impl<'d> LatentJointContext<'d> {
8170 fn eval_full(
8171 &mut self,
8172 theta: &Array1<f64>,
8173 order: OuterEvalOrder,
8174 ) -> Result<(f64, Array1<f64>, gam_problem::HessianValue), EstimationError> {
8175 if let Some(eval) = self.cache.memoized_eval(theta) {
8176 return Ok(eval);
8177 }
8178 self.cache
8179 .ensure_theta(theta)
8180 .map_err(EstimationError::InvalidInput)?;
8181 let hyper_dirs = self
8182 .cache
8183 .hyper_dirs()
8184 .map_err(EstimationError::InvalidInput)?;
8185 let design_revision = Some(self.cache.design_revision());
8186 let registry_for_key = self.cache.analytic_penalties();
8187 self.evaluator
8188 .set_analytic_penalty_registry(registry_for_key.as_deref());
8189 let mut eval = evaluate_joint_reml_outer_eval_at_theta(
8190 &mut self.evaluator,
8191 self.cache.design(),
8192 theta,
8193 self.rho_dim,
8194 hyper_dirs,
8195 None,
8196 order,
8197 design_revision,
8198 )?;
8199 let latent = self.cache.latent().map_err(EstimationError::InvalidInput)?;
8200 if let Some(registry) = registry_for_key {
8201 add_analytic_penalty_objective_to_eval(
8202 theta,
8203 self.rho_dim,
8204 latent.as_ref(),
8205 registry.as_ref(),
8206 &mut eval,
8207 )?;
8208 }
8209 add_latent_id_objective_to_eval(
8210 theta,
8211 self.rho_dim,
8212 self.cache.analytic_penalty_rho_count(),
8213 latent.as_ref(),
8214 &mut eval,
8215 )?;
8216 self.cache.store_eval(eval.clone());
8217 Ok(eval)
8218 }
8219
8220 fn eval_efs(
8221 &mut self,
8222 theta: &Array1<f64>,
8223 ) -> Result<gam_problem::EfsEval, EstimationError> {
8224 self.cache
8225 .ensure_theta(theta)
8226 .map_err(EstimationError::InvalidInput)?;
8227 let hyper_dirs = self
8228 .cache
8229 .hyper_dirs()
8230 .map_err(EstimationError::InvalidInput)?;
8231 let registry_for_key = self.cache.analytic_penalties();
8232 self.evaluator
8233 .set_analytic_penalty_registry(registry_for_key.as_deref());
8234 let mut efs = evaluate_joint_reml_efs_at_theta(
8235 &mut self.evaluator,
8236 self.cache.design(),
8237 theta,
8238 self.rho_dim,
8239 hyper_dirs,
8240 None,
8241 Some(self.cache.design_revision()),
8242 )?;
8243 if let Some(registry) = registry_for_key {
8244 let latent = self.cache.latent().map_err(EstimationError::InvalidInput)?;
8245 let contribution = analytic_penalty_objective_contribution(
8246 theta,
8247 self.rho_dim,
8248 latent.as_ref(),
8249 registry.as_ref(),
8250 )?;
8251 efs.cost += contribution.cost;
8252 if let (Some(psi_gradient), Some(psi_indices)) =
8253 (efs.psi_gradient.as_mut(), efs.psi_indices.as_ref())
8254 {
8255 if psi_gradient.len() != psi_indices.len() {
8256 crate::bail_invalid_estim!(
8257 "latent-coordinate analytic penalty EFS psi gradient length mismatch: gradient={}, indices={}",
8258 psi_gradient.len(),
8259 psi_indices.len()
8260 );
8261 }
8262 for (local_idx, &theta_idx) in psi_indices.iter().enumerate() {
8263 psi_gradient[local_idx] += contribution.gradient[theta_idx];
8264 }
8265 }
8266 }
8267 Ok(efs)
8268 }
8269
8270 fn eval_cost(&mut self, theta: &Array1<f64>) -> f64 {
8271 if let Some(cost) = self.cache.memoized_cost(theta) {
8272 return cost;
8273 }
8274 if self.cache.ensure_theta(theta).is_err() {
8275 return f64::INFINITY;
8276 }
8277 let design_revision = Some(self.cache.design_revision());
8278 let registry_for_key = self.cache.analytic_penalties();
8279 self.evaluator
8280 .set_analytic_penalty_registry(registry_for_key.as_deref());
8281 let result = {
8282 let design = self.cache.design();
8283 self.evaluator.evaluate_cost_only(
8284 &design.design,
8285 &design.penalties,
8286 &design.nullspace_dims,
8287 design.linear_constraints.clone(),
8288 theta,
8289 self.rho_dim,
8290 None,
8291 "latent-coordinate-joint cost-only",
8292 design_revision,
8293 )
8294 };
8295 match result {
8296 Ok(cost) => {
8297 let latent = match self.cache.latent() {
8298 Ok(latent) => latent,
8299 Err(_) => return f64::INFINITY,
8300 };
8301 let contribution = match latent_id_objective_contribution(
8302 theta,
8303 self.rho_dim,
8304 self.cache.analytic_penalty_rho_count(),
8305 latent.as_ref(),
8306 ) {
8307 Ok(contribution) => contribution,
8308 Err(_) => return f64::INFINITY,
8309 };
8310 let cost = cost + contribution.cost;
8311 let cost = if let Some(registry) = registry_for_key {
8312 match analytic_penalty_objective_contribution(
8313 theta,
8314 self.rho_dim,
8315 latent.as_ref(),
8316 registry.as_ref(),
8317 ) {
8318 Ok(contribution) => cost + contribution.cost,
8319 Err(_) => return f64::INFINITY,
8320 }
8321 } else {
8322 cost
8323 };
8324 self.cache.store_cost(cost);
8325 cost
8326 }
8327 Err(_) => f64::INFINITY,
8328 }
8329 }
8330 }
8331
8332 let effective_offset = best
8333 .design
8334 .compose_offset(offset, "latent-coordinate joint fit")
8335 .map_err(EstimationError::BasisError)?;
8336 let mut ctx = LatentJointContext {
8337 rho_dim,
8338 cache: SingleBlockLatentCoordDesignCache::new(
8339 data.to_owned(),
8340 resolvedspec.clone(),
8341 best.design.clone(),
8342 latent,
8343 rho_dim,
8344 )
8345 .map_err(EstimationError::InvalidInput)?,
8346 evaluator: gam_solve::estimate::ExternalJointHyperEvaluator::new(
8347 y,
8348 weights,
8349 &best.design.design,
8350 effective_offset.view(),
8351 &best.design.penalties,
8352 &external_opts_for_design(&family, &best.design, options),
8353 "latent-coordinate-joint",
8354 )?,
8355 };
8356 let registry_for_key = ctx.cache.analytic_penalties();
8357 ctx.evaluator
8358 .set_analytic_penalty_registry(registry_for_key.as_deref());
8359 ctx.evaluator
8360 .set_persistent_latent_values_fingerprint(latent.values.id_mode());
8361 if let Some(cached_t) = ctx
8362 .evaluator
8363 .load_persistent_latent_values(latent.values.n_obs(), latent.values.latent_dim())
8364 {
8365 let cached_t: Array2<f64> = cached_t;
8366 for (dst, src) in theta0
8367 .slice_mut(s![rho_dim..rho_dim + latent_flat_dim])
8368 .iter_mut()
8369 .zip(cached_t.iter())
8370 {
8371 *dst = *src;
8372 }
8373 }
8374
8375 let problem = exact_joint_multistart_outer_problem(
8376 &theta0,
8377 &lower,
8378 &upper,
8379 rho_dim,
8380 latent_coord_ext_dim,
8381 theta0.len(),
8382 Derivative::Analytic,
8383 DeclaredHessianForm::Unavailable,
8384 false,
8385 seed_risk_profile_for_likelihood_family(&family),
8386 options.tol,
8387 options.max_iter.max(1),
8388 Some(5.0),
8389 Some(0.5),
8390 None,
8391 Some((data.nrows(), best.design.design.ncols().max(1))),
8394 !constant_curvature_term_indices(resolvedspec).is_empty(),
8397 false,
8399 )?;
8400
8401 let eval_outer = |ctx: &mut &mut LatentJointContext<'_>,
8402 theta: &Array1<f64>,
8403 order: OuterEvalOrder|
8404 -> Result<OuterEval, EstimationError> {
8405 let (cost, gradient, hessian) = ctx.eval_full(theta, order)?;
8406 Ok(OuterEval {
8407 cost,
8408 gradient,
8409 hessian,
8410 inner_beta_hint: None,
8411 })
8412 };
8413
8414 let result = {
8415 let obj = problem.build_objective_with_eval_order(
8416 &mut ctx,
8417 |ctx: &mut &mut LatentJointContext<'_>, theta: &Array1<f64>| Ok(ctx.eval_cost(theta)),
8418 |ctx: &mut &mut LatentJointContext<'_>, theta: &Array1<f64>| {
8419 eval_outer(ctx, theta, OuterEvalOrder::ValueAndGradient)
8420 },
8421 |ctx: &mut &mut LatentJointContext<'_>, theta: &Array1<f64>, order: OuterEvalOrder| {
8422 eval_outer(ctx, theta, order)
8423 },
8424 Some(|ctx: &mut &mut LatentJointContext<'_>| {
8425 ctx.cache.reset();
8426 }),
8427 Some(|ctx: &mut &mut LatentJointContext<'_>, theta: &Array1<f64>| ctx.eval_efs(theta)),
8428 );
8429 let mut obj = obj.with_criterion_invariance(
8434 |ctx: &mut &mut LatentJointContext<'_>, rho: &Array1<f64>| {
8435 ctx.evaluator.criterion_invariant_directions(rho)
8436 },
8437 );
8438
8439 problem
8440 .run(&mut obj, "latent-coordinate joint REML")
8441 .map_err(|e| {
8442 EstimationError::InvalidInput(format!(
8443 "latent-coordinate joint optimization failed after exhausting strategy fallbacks: {e}"
8444 ))
8445 })?
8446 };
8447 if !result.converged() {
8448 crate::bail_invalid_estim!(
8449 "latent-coordinate joint optimization did not converge after {} iterations (final_objective={:.6e}, final_grad_norm={})",
8450 result.iterations,
8451 result.final_value,
8452 result.final_grad_norm_report(),
8453 );
8454 }
8455
8456 let theta_star = result.rho;
8457 let selected_lambdas = Array1::from_vec(
8458 gam_problem::checked_exp_log_strengths(
8459 theta_star.slice(s![..rho_dim]).iter().copied(),
8460 )
8461 .map_err(|error| {
8462 EstimationError::InvalidInput(format!(
8463 "selected latent-coordinate smoothing coordinate is outside the canonical log-strength domain: {error}"
8464 ))
8465 })?,
8466 );
8467 let mut final_data = data.to_owned();
8468 let flat_t = theta_star
8469 .slice(s![rho_dim..rho_dim + latent_flat_dim])
8470 .to_owned();
8471 let mut fitted_latent_values =
8472 Array2::<f64>::zeros((latent.values.n_obs(), latent.values.latent_dim()));
8473 for n in 0..latent.values.n_obs() {
8474 for axis in 0..latent.values.latent_dim() {
8475 let value = flat_t[n * latent.values.latent_dim() + axis];
8476 fitted_latent_values[[n, axis]] = value;
8477 final_data[[n, latent.feature_cols[axis]]] = value;
8478 }
8479 }
8480 let optimized = fit_term_collection_forspecwith_heuristic_lambdas(
8481 final_data.view(),
8482 y,
8483 weights,
8484 offset,
8485 resolvedspec,
8486 selected_lambdas.as_slice(),
8487 family,
8488 options,
8489 )?;
8490 ctx.evaluator
8491 .store_persistent_latent_values(&fitted_latent_values);
8492 let mut fit = optimized.fit;
8493 fit.set_criterion(Some(result.final_value));
8494 Ok(FittedTermCollectionWithSpec {
8495 fit,
8496 design: optimized.design,
8497 resolvedspec: resolvedspec.clone(),
8498 adaptive_diagnostics: optimized.adaptive_diagnostics,
8499 kappa_timing: None,
8500 })
8501}
8502
8503pub fn fit_term_collectionwith_latent_coord_optimization(
8504 data: ArrayView2<'_, f64>,
8505 y: Array1<f64>,
8506 weights: Array1<f64>,
8507 offset: Array1<f64>,
8508 spec: &TermCollectionSpec,
8509 latent: &StandardLatentCoordConfig,
8510 family: LikelihoodSpec,
8511 options: &FitOptions,
8512) -> Result<FittedTermCollectionWithSpec, EstimationError> {
8513 let n = data.nrows();
8514 if !(y.len() == n && weights.len() == n && offset.len() == n) {
8515 crate::bail_invalid_estim!(
8516 "fit_term_collectionwith_latent_coord_optimization row mismatch: n={}, y={}, weights={}, offset={}",
8517 n,
8518 y.len(),
8519 weights.len(),
8520 offset.len()
8521 );
8522 }
8523 let best = fit_term_collection_forspec(
8524 data,
8525 y.view(),
8526 weights.view(),
8527 offset.view(),
8528 spec,
8529 family.clone(),
8530 options,
8531 )?;
8532 let resolvedspec = freeze_term_collection_from_design(spec, &best.design)?;
8533 try_exact_joint_latent_coord_optimization(
8534 data,
8535 y.view(),
8536 weights.view(),
8537 offset.view(),
8538 &resolvedspec,
8539 &best,
8540 family,
8541 options,
8542 latent,
8543 )
8544}
8545
8546fn select_isotropic_matern_range_basin(
8563 data: ArrayView2<'_, f64>,
8564 y: ArrayView1<'_, f64>,
8565 weights: ArrayView1<'_, f64>,
8566 offset: ArrayView1<'_, f64>,
8567 mut resolvedspec: TermCollectionSpec,
8568 mut best: FittedTermCollection,
8569 family: &LikelihoodSpec,
8570 options: &FitOptions,
8571 kappa_options: &SpatialLengthScaleOptimizationOptions,
8572 spatial_terms: &[usize],
8573) -> Result<(TermCollectionSpec, FittedTermCollection), EstimationError> {
8574 if has_aniso_terms(&resolvedspec, spatial_terms)
8578 || !constant_curvature_term_indices(&resolvedspec).is_empty()
8579 {
8580 return Ok((resolvedspec, best));
8581 }
8582
8583 let mut best_score = fit_score(&best.fit);
8584 if !best_score.is_finite() {
8585 crate::bail_invalid_estim!(
8586 "isotropic Matérn basin selection received a non-finite incumbent profile"
8587 );
8588 }
8589
8590 for &term_idx in spatial_terms {
8591 let Some(SmoothBasisSpec::Matern {
8592 feature_cols,
8593 spec: matern,
8594 ..
8595 }) = resolvedspec
8596 .smooth_terms
8597 .get(term_idx)
8598 .map(|term| &term.basis)
8599 else {
8600 continue;
8601 };
8602 let num_centers = gam_terms::basis::center_strategy_num_centers(&matern.center_strategy)
8603 .ok_or_else(|| {
8604 EstimationError::InvalidInput(format!(
8605 "resolved isotropic Matérn term {term_idx} has no finite center count"
8606 ))
8607 })?;
8608 let companion_length_scale = matern_low_rank_center_resolution_length_scale(
8609 data,
8610 feature_cols,
8611 num_centers,
8612 )
8613 .ok_or_else(|| {
8614 EstimationError::InvalidInput(format!(
8615 "resolved isotropic Matérn term {term_idx} has no finite center-resolution range"
8616 ))
8617 })?;
8618 let (psi_long_bound, psi_short_bound) =
8619 spatial_term_psi_bounds(data, &resolvedspec, term_idx, kappa_options)
8620 .map_err(EstimationError::BasisError)?;
8621 let psi_long = (-companion_length_scale.ln()).clamp(psi_long_bound, psi_short_bound);
8622 let long_length_scale = (-psi_long).exp();
8623 if !(long_length_scale.is_finite() && long_length_scale > 0.0) {
8624 crate::bail_invalid_estim!(
8625 "isotropic Matérn term {term_idx} produced an invalid long-range endpoint from psi={psi_long}"
8626 );
8627 }
8628 if get_spatial_length_scale(&resolvedspec, term_idx)
8629 .is_some_and(|current| current == long_length_scale)
8630 {
8631 continue;
8632 }
8633
8634 let mut endpoint_spec = resolvedspec.clone();
8635 set_spatial_length_scale(&mut endpoint_spec, term_idx, long_length_scale)?;
8636 let endpoint = fit_term_collection_forspecwith_heuristic_lambdas(
8646 data,
8647 y,
8648 weights,
8649 offset,
8650 &endpoint_spec,
8651 best.fit.lambdas.as_slice(),
8652 family.clone(),
8653 options,
8654 )?;
8655 let endpoint_score = fit_score(&endpoint.fit);
8656 if !endpoint_score.is_finite() {
8657 crate::bail_invalid_estim!(
8658 "isotropic Matérn term {term_idx} long-range endpoint returned a non-finite profiled REML score"
8659 );
8660 }
8661
8662 if endpoint_score < best_score {
8663 log::info!(
8664 "[spatial-kappa] term {term_idx} selected certified long-range basin: \
8665 length_scale={long_length_scale:.6}, profiled REML {endpoint_score:.6} \
8666 < short-basin {best_score:.6}"
8667 );
8668 resolvedspec = freeze_term_collection_from_design(&endpoint_spec, &endpoint.design)?;
8669 best = endpoint;
8670 best_score = endpoint_score;
8671 } else {
8672 log::info!(
8673 "[spatial-kappa] term {term_idx} retained certified short-range basin: \
8674 profiled REML {best_score:.6} <= long-endpoint {endpoint_score:.6} \
8675 at length_scale={long_length_scale:.6}"
8676 );
8677 }
8678 }
8679
8680 Ok((resolvedspec, best))
8681}
8682
8683pub fn fit_term_collectionwith_spatial_length_scale_optimization(
8684 data: ArrayView2<'_, f64>,
8685 y: Array1<f64>,
8686 weights: Array1<f64>,
8687 offset: Array1<f64>,
8688 spec: &TermCollectionSpec,
8689 family: LikelihoodSpec,
8690 options: &FitOptions,
8691 kappa_options: &SpatialLengthScaleOptimizationOptions,
8692) -> Result<FittedTermCollectionWithSpec, EstimationError> {
8693 let mut resolvedspec = spec.clone();
8709 let n = data.nrows();
8710 if !(y.len() == n && weights.len() == n && offset.len() == n) {
8711 crate::bail_invalid_estim!(
8712 "fit_term_collectionwith_spatial_length_scale_optimization row mismatch: n={}, y={}, weights={}, offset={}",
8713 n,
8714 y.len(),
8715 weights.len(),
8716 offset.len()
8717 );
8718 }
8719 seed_measure_jet_auto_ranges(data, y.view(), weights.view(), &mut resolvedspec);
8728 let spatial_terms = spatial_length_scale_term_indices(&resolvedspec);
8729 if !kappa_options.enabled || spatial_terms.is_empty() {
8730 let out = fit_term_collection_forspec(
8731 data,
8732 y.view(),
8733 weights.view(),
8734 offset.view(),
8735 &resolvedspec,
8736 family,
8737 options,
8738 )?;
8739 let resolvedspec = freeze_term_collection_from_design(&resolvedspec, &out.design)?;
8740 return Ok(FittedTermCollectionWithSpec {
8741 fit: out.fit,
8742 design: out.design,
8743 resolvedspec,
8744 adaptive_diagnostics: out.adaptive_diagnostics,
8745 kappa_timing: None,
8746 });
8747 }
8748 if kappa_options.max_outer_iter == 0 {
8749 crate::bail_invalid_estim!("spatial kappa optimization requires max_outer_iter >= 1");
8750 }
8751 if !(kappa_options.log_step.is_finite() && kappa_options.log_step > 0.0) {
8752 crate::bail_invalid_estim!("spatial kappa optimization requires log_step > 0");
8753 }
8754 if !(kappa_options.min_length_scale.is_finite()
8755 && kappa_options.max_length_scale.is_finite()
8756 && kappa_options.min_length_scale > 0.0
8757 && kappa_options.max_length_scale >= kappa_options.min_length_scale)
8758 {
8759 crate::bail_invalid_estim!(
8760 "spatial kappa optimization requires valid positive length_scale bounds"
8761 );
8762 }
8763
8764 let projected_scales =
8781 project_spatial_length_scales_in_spec(&mut resolvedspec, &spatial_terms, kappa_options)?;
8782 for &(term_idx, raw, projected) in &projected_scales {
8783 log::info!(
8784 "[spatial-kappa] term {term_idx}: length_scale projected onto the caller's window \
8785 before the baseline fit: {raw:.6e} -> {projected:.6e} \
8786 (window=[{:.6e}, {:.6e}])",
8787 kappa_options.min_length_scale,
8788 kappa_options.max_length_scale,
8789 );
8790 }
8791
8792 let pilot_threshold = kappa_options.pilot_subsample_threshold;
8793 if pilot_threshold > 0 && n > pilot_threshold * 2 {
8794 log::info!(
8795 "[spatial-kappa] n={n} exceeds pilot threshold {}; using pilot geometry only for deterministic anisotropy initialization",
8796 pilot_threshold * 2,
8797 );
8798 apply_spatial_anisotropy_pilot_initializer(
8799 data,
8800 &mut resolvedspec,
8801 &spatial_terms,
8802 pilot_threshold,
8803 kappa_options,
8804 )?;
8805 }
8806
8807 apply_response_aware_anisotropy_seed(data, y.view(), &mut resolvedspec, &spatial_terms);
8816
8817 let free_curvature_terms: Vec<usize> = constant_curvature_term_indices(&resolvedspec)
8852 .into_iter()
8853 .filter(|&term_idx| !constant_curvature_kappa_is_fixed(&resolvedspec, term_idx))
8854 .collect();
8855 let pinned_kappa_free_range_terms: Vec<usize> =
8856 constant_curvature_term_indices(&resolvedspec)
8857 .into_iter()
8858 .filter(|&term_idx| {
8859 constant_curvature_kappa_is_fixed(&resolvedspec, term_idx)
8860 && !constant_curvature_length_scale_is_fixed(&resolvedspec, term_idx)
8861 })
8862 .collect();
8863 if !free_curvature_terms.is_empty() {
8864 validate_constant_curvature_profile_inputs(weights.view(), offset.view(), &family)?;
8865 }
8866 if !pinned_kappa_free_range_terms.is_empty()
8867 && validate_constant_curvature_profile_inputs(weights.view(), offset.view(), &family)
8868 .is_ok()
8869 {
8870 for term_idx in pinned_kappa_free_range_terms {
8871 let length_scale_hat =
8872 constant_curvature_range_only_optimum(data, y.view(), &resolvedspec, term_idx)?;
8873 if let Some(SmoothBasisSpec::ConstantCurvature { spec: cc, .. }) = resolvedspec
8874 .smooth_terms
8875 .get_mut(term_idx)
8876 .map(|term| &mut term.basis)
8877 {
8878 cc.length_scale = length_scale_hat;
8882 }
8883 }
8884 }
8885 for term_idx in free_curvature_terms {
8886 let psi_hat = constant_curvature_kappa_profile_optimum(
8887 data,
8888 y.view(),
8889 &resolvedspec,
8890 term_idx,
8891 options,
8892 )?;
8893 if let Some(SmoothBasisSpec::ConstantCurvature { spec: cc, .. }) = resolvedspec
8894 .smooth_terms
8895 .get_mut(term_idx)
8896 .map(|term| &mut term.basis)
8897 {
8898 cc.kappa = psi_hat.kappa;
8899 cc.length_scale = psi_hat.length_scale;
8904 }
8905 }
8906
8907 let baseline_options = superseded_fit_options(options);
8908 let best = fit_term_collection_forspec(
8909 data,
8910 y.view(),
8911 weights.view(),
8912 offset.view(),
8913 &resolvedspec,
8914 family.clone(),
8915 &baseline_options,
8916 )?;
8917 resolvedspec = freeze_term_collection_from_design(&resolvedspec, &best.design)?;
8918 let spatial_terms: Vec<usize> = spatial_length_scale_term_indices(&resolvedspec)
8936 .into_iter()
8937 .filter(|&term_idx| constant_curvature_term_spec(&resolvedspec, term_idx).is_none())
8938 .collect();
8939 let (next_spec, best) = select_isotropic_matern_range_basin(
8940 data,
8941 y.view(),
8942 weights.view(),
8943 offset.view(),
8944 resolvedspec,
8945 best,
8946 &family,
8947 &baseline_options,
8948 kappa_options,
8949 &spatial_terms,
8950 )?;
8951 resolvedspec = next_spec;
8952 sync_aniso_contrasts_from_metadata(&mut resolvedspec, &best.design.smooth);
8956 if spatial_terms.is_empty() {
8957 let fitted = fit_term_collection_forspecwith_heuristic_lambdas(
8958 data,
8959 y.view(),
8960 weights.view(),
8961 offset.view(),
8962 &resolvedspec,
8963 best.fit.lambdas.as_slice(),
8964 family,
8965 options,
8966 )?;
8967 return Ok(FittedTermCollectionWithSpec {
8968 fit: fitted.fit,
8969 design: fitted.design,
8970 resolvedspec,
8971 adaptive_diagnostics: fitted.adaptive_diagnostics,
8972 kappa_timing: None,
8973 });
8974 }
8975 let initial_score = fit_score(&best.fit);
8976 if !initial_score.is_finite() {
8977 crate::bail_invalid_estim!(
8978 "spatial kappa optimization received a non-finite initial profiled score"
8979 );
8980 }
8981 let exact_joint = match try_exact_joint_spatial_length_scale_optimization(
8982 data,
8983 y.view(),
8984 weights.view(),
8985 offset.view(),
8986 &resolvedspec,
8987 &best,
8988 family.clone(),
8989 options,
8990 kappa_options,
8991 &spatial_terms,
8992 )? {
8993 JointSpatialKappaOutcome::Optimized(optimized) => *optimized,
8994 JointSpatialKappaOutcome::DeclinedKeepIncumbent {
8995 baseline_score,
8996 optimized_score,
8997 } => {
8998 log::info!(
9007 "[spatial-kappa] joint kappa optimization DECLINED its own candidate (incumbent={baseline_score:.12e}, candidate={optimized_score:.12e}, regression={:.3e}); shipping the incumbent scalar-route fit at the incumbent κ, which is what the decline means. Not an unavailability.",
9008 optimized_score - baseline_score,
9009 );
9010 let fitted = fit_term_collection_forspecwith_heuristic_lambdas(
9011 data,
9012 y.view(),
9013 weights.view(),
9014 offset.view(),
9015 &resolvedspec,
9016 best.fit.lambdas.as_slice(),
9017 family,
9018 options,
9019 )?;
9020 return Ok(FittedTermCollectionWithSpec {
9021 fit: fitted.fit,
9022 design: fitted.design,
9023 resolvedspec,
9024 adaptive_diagnostics: fitted.adaptive_diagnostics,
9025 kappa_timing: None,
9026 });
9027 }
9028 JointSpatialKappaOutcome::Unavailable => {
9029 return Err(EstimationError::RemlOptimizationFailed(
9030 "spatial kappa optimization is unavailable for one or more eligible spatial terms"
9031 .to_string(),
9032 ));
9033 }
9034 };
9035 let exact_joint = require_available_spatial_optimization_result(Ok(Some(exact_joint)))?;
9036 let exact_score = fit_score(&exact_joint.fit);
9037
9038 if exact_score.is_finite() && exact_score <= initial_score {
9054 log_spatial_aniso_scales(&exact_joint.resolvedspec);
9055 return Ok(exact_joint);
9056 }
9057 log::info!(
9058 "[spatial-kappa] the optimized-κ fit scores {exact_score:.12e} against the incumbent's \
9059 {initial_score:.12e} (regression {:.3e}); shipping the INCUMBENT, which is the better \
9060 of the two fits this call has in hand. A refinement that does not improve on the fit \
9061 it refines is not a reason to have no fit (#2748).",
9062 exact_score - initial_score,
9063 );
9064 let fitted = fit_term_collection_forspecwith_heuristic_lambdas(
9065 data,
9066 y.view(),
9067 weights.view(),
9068 offset.view(),
9069 &resolvedspec,
9070 best.fit.lambdas.as_slice(),
9071 family,
9072 options,
9073 )?;
9074 Ok(FittedTermCollectionWithSpec {
9075 fit: fitted.fit,
9076 design: fitted.design,
9077 resolvedspec,
9078 adaptive_diagnostics: fitted.adaptive_diagnostics,
9079 kappa_timing: None,
9080 })
9081}
9082
9083#[derive(Clone, Debug)]
9089pub struct CurvatureInference {
9090 pub term_idx: usize,
9092 pub kappa_hat: f64,
9095 pub ci: gam_geometry::curvature_estimand::KappaProfileCi,
9097 pub flatness: gam_geometry::curvature_estimand::FlatnessTest,
9101 pub length_scale_hat: f64,
9111 pub length_scale_estimated: bool,
9113 pub length_scale_support: gam_geometry::curvature_estimand::RangeEstimateSupport,
9121}
9122
9123fn curvature_profile_lr_endpoint<F>(
9135 profile: &mut F,
9136 kappa_hat: f64,
9137 value_hat: f64,
9138 bound: f64,
9139 half_threshold: f64,
9140 x_tolerance: f64,
9141 score_tolerance: f64,
9142) -> Result<(f64, bool), String>
9143where
9144 F: FnMut(f64) -> Result<(f64, f64), String>,
9145{
9146 let direction = (bound - kappa_hat).signum();
9147 let span = (bound - kappa_hat).abs();
9148 if direction == 0.0 || span <= x_tolerance {
9149 return Ok((bound, true));
9150 }
9151
9152 let (bound_value, bound_score) = profile(bound)?;
9153 let outward_score = direction * bound_score;
9154 if outward_score < -score_tolerance {
9155 return Err(format!(
9156 "curvature profile is not outward-monotone at chart bound {bound}: \
9157 outward score {outward_score:.6e} is below tolerance {score_tolerance:.6e}"
9158 ));
9159 }
9160 let value_tolerance = score_tolerance * span;
9161 if bound_value < value_hat - value_tolerance {
9162 return Err(format!(
9163 "fitted curvature is not the minimum of its inference profile: \
9164 V(bound={bound})={bound_value:.6e} < V(kappa_hat)={value_hat:.6e}"
9165 ));
9166 }
9167 let bound_residual = bound_value - value_hat - half_threshold;
9168 if bound_residual < 0.0 {
9169 return Ok((bound, true));
9170 }
9171 if bound_residual == 0.0 {
9172 return Ok((bound, false));
9173 }
9174
9175 let mut inside_x = kappa_hat;
9180 let mut outside_x = bound;
9181 let mut outside_residual = bound_residual;
9182 let mut outside_score = bound_score;
9183 while (outside_x - inside_x).abs() > x_tolerance {
9184 let lo = inside_x.min(outside_x);
9185 let hi = inside_x.max(outside_x);
9186 let width = hi - lo;
9187 let central_lo = lo + 0.25 * width;
9188 let central_hi = hi - 0.25 * width;
9189 let newton = outside_x - outside_residual / outside_score;
9190 let probe = if newton.is_finite() && newton > central_lo && newton < central_hi {
9191 newton
9192 } else {
9193 lo + 0.5 * width
9194 };
9195 if !(probe > lo && probe < hi) {
9196 break;
9197 }
9198 let (value, score) = profile(probe)?;
9199 let outward_score = direction * score;
9200 if outward_score < -score_tolerance {
9201 return Err(format!(
9202 "curvature profile changed direction before its likelihood crossing at \
9203 kappa={probe}: outward score {outward_score:.6e} is below tolerance \
9204 {score_tolerance:.6e}"
9205 ));
9206 }
9207 let residual = value - value_hat - half_threshold;
9208 if residual >= 0.0 {
9209 outside_x = probe;
9210 outside_residual = residual;
9211 outside_score = score;
9212 } else {
9213 inside_x = probe;
9214 }
9215 }
9216 let midpoint = inside_x + 0.5 * (outside_x - inside_x);
9226 let refined = outside_x - outside_residual / outside_score;
9227 let lo = inside_x.min(outside_x);
9228 let hi = inside_x.max(outside_x);
9229 let endpoint = if refined.is_finite() && refined >= lo && refined <= hi {
9230 refined
9231 } else {
9232 midpoint
9233 };
9234 Ok((endpoint, false))
9235}
9236
9237fn curvature_profile_ci_from_analytic_score<F>(
9238 profile: &mut F,
9239 kappa_hat: f64,
9240 kappa_min: f64,
9241 kappa_max: f64,
9242 level: f64,
9243 relative_tolerance: f64,
9244) -> Result<gam_geometry::curvature_estimand::KappaProfileCi, String>
9245where
9246 F: FnMut(f64) -> Result<(f64, f64), String>,
9247{
9248 if !(kappa_min < kappa_max && kappa_hat >= kappa_min && kappa_hat <= kappa_max) {
9249 return Err("curvature profile requires kappa_hat inside valid chart bounds".to_string());
9250 }
9251 if !(level > 0.0 && level < 1.0) {
9252 return Err("curvature profile level must lie in (0, 1)".to_string());
9253 }
9254 let z = gam_geometry::curvature_estimand::wald_half_width(1.0, level)
9255 .ok_or_else(|| "curvature profile threshold is not finite".to_string())?;
9256 let half_threshold = 0.5 * z * z;
9257 let (value_hat, score_hat) = profile(kappa_hat)?;
9258 let relative_tolerance = relative_tolerance.max(f64::EPSILON.sqrt());
9259 let x_tolerance = relative_tolerance * (1.0 + kappa_min.abs().max(kappa_max.abs()));
9260 let score_tolerance = relative_tolerance * (1.0 + value_hat.abs());
9261 let at_lower = (kappa_hat - kappa_min).abs() <= x_tolerance;
9267 let at_upper = (kappa_hat - kappa_max).abs() <= x_tolerance;
9268 let kappa_hat_support = if at_lower {
9269 gam_geometry::curvature_estimand::KappaEstimateSupport::RailedAtLowerBound
9270 } else if at_upper {
9271 gam_geometry::curvature_estimand::KappaEstimateSupport::RailedAtUpperBound
9272 } else {
9273 gam_geometry::curvature_estimand::KappaEstimateSupport::Interior
9274 };
9275 let stationary = if at_lower {
9276 score_hat >= -score_tolerance
9277 } else if at_upper {
9278 score_hat <= score_tolerance
9279 } else {
9280 score_hat.abs() <= score_tolerance
9281 };
9282 if !stationary {
9283 return Err(format!(
9289 "curvature inference rejected a non-stationary point estimate: \
9290 kappa_hat={kappa_hat}, score={score_hat:.6e}, \
9291 stationarity_bound={score_tolerance:.6e}; \
9292 box=[{kappa_min}, {kappa_max}], gap_to_lower={:.6e}, gap_to_upper={:.6e}, \
9293 rail_tolerance={x_tolerance:.6e}, classified={}",
9294 kappa_hat - kappa_min,
9295 kappa_max - kappa_hat,
9296 kappa_hat_support.label()
9297 ));
9298 }
9299
9300 let (ci_lo, lo_at_bound) = curvature_profile_lr_endpoint(
9301 profile,
9302 kappa_hat,
9303 value_hat,
9304 kappa_min,
9305 half_threshold,
9306 x_tolerance,
9307 score_tolerance,
9308 )?;
9309 let (ci_hi, hi_at_bound) = curvature_profile_lr_endpoint(
9310 profile,
9311 kappa_hat,
9312 value_hat,
9313 kappa_max,
9314 half_threshold,
9315 x_tolerance,
9316 score_tolerance,
9317 )?;
9318 let verdict = if ci_lo > 0.0 {
9319 gam_geometry::curvature_estimand::CurvatureVerdict::Spherical
9320 } else if ci_hi < 0.0 {
9321 gam_geometry::curvature_estimand::CurvatureVerdict::Hyperbolic
9322 } else {
9323 gam_geometry::curvature_estimand::CurvatureVerdict::Flat
9324 };
9325 Ok(gam_geometry::curvature_estimand::KappaProfileCi {
9326 kappa_hat,
9327 ci_lo,
9328 ci_hi,
9329 lo_at_bound,
9330 hi_at_bound,
9331 kappa_hat_support,
9332 verdict,
9333 })
9334}
9335
9336pub fn curvature_inference_forspec(
9337 data: ArrayView2<'_, f64>,
9338 y: ArrayView1<'_, f64>,
9339 weights: ArrayView1<'_, f64>,
9340 offset: ArrayView1<'_, f64>,
9341 resolvedspec: &TermCollectionSpec,
9342 term_idx: usize,
9343 family: LikelihoodSpec,
9344 options: &FitOptions,
9345 level: f64,
9346) -> Result<CurvatureInference, EstimationError> {
9347 let kappa_hat = get_constant_curvature_kappa(resolvedspec, term_idx).ok_or_else(|| {
9348 EstimationError::InvalidInput(format!(
9349 "curvature_inference_forspec: term {term_idx} is not a constant-curvature smooth"
9350 ))
9351 })?;
9352 if constant_curvature_kappa_is_fixed(resolvedspec, term_idx) {
9353 crate::bail_invalid_estim!(
9354 "curvature inference requires an estimated curvature; term {term_idx} has user-pinned kappa={kappa_hat}"
9355 );
9356 }
9357 if y.len() != data.nrows() || weights.len() != data.nrows() || offset.len() != data.nrows() {
9358 crate::bail_invalid_estim!(
9359 "curvature inference row mismatch: data={}, y={}, weights={}, offset={}",
9360 data.nrows(),
9361 y.len(),
9362 weights.len(),
9363 offset.len(),
9364 );
9365 }
9366 validate_constant_curvature_profile_inputs(weights, offset, &family)?;
9367 let (kappa_min, kappa_max) = constant_curvature_kappa_bounds(data, resolvedspec, term_idx);
9368 let (feature_cols, base_spec) = match resolvedspec
9369 .smooth_terms
9370 .get(term_idx)
9371 .map(|term| &term.basis)
9372 {
9373 Some(SmoothBasisSpec::ConstantCurvature {
9374 feature_cols, spec, ..
9375 }) => (feature_cols, spec.clone()),
9376 _ => {
9377 return Err(EstimationError::InvalidInput(format!(
9378 "constant-curvature κ profile: smooth term {term_idx} is not a \
9379 constant-curvature basis"
9380 )));
9381 }
9382 };
9383 let x_term = select_columns(data, feature_cols).map_err(EstimationError::from)?;
9384 let profile = ConstantCurvatureProfile::new(x_term.view(), y, base_spec)?;
9385
9386 let mut v_p = |kappa: f64| -> Result<(f64, f64), String> {
9389 if !kappa.is_finite() {
9390 return Err(format!("V_p probed a non-finite κ = {kappa}"));
9391 }
9392 let (value, score, _curvature) = profile.evaluate(kappa).map_err(|error| {
9393 format!("analytic curvature profile at kappa={kappa} failed: {error}")
9394 })?;
9395 Ok((value, score))
9396 };
9397 let ci = curvature_profile_ci_from_analytic_score(
9398 &mut v_p,
9399 kappa_hat,
9400 kappa_min,
9401 kappa_max,
9402 level,
9403 options.tol,
9404 )
9405 .map_err(EstimationError::RemlOptimizationFailed)?;
9406 let flatness = gam_geometry::curvature_estimand::flatness_lr_test(
9407 |kappa| v_p(kappa).map(|(value, _)| value),
9408 kappa_hat,
9409 )
9410 .map_err(EstimationError::RemlOptimizationFailed)?;
9411
9412 let (eta_hat, _, range_outcome) = profile.minimize_over_eta(kappa_hat)?;
9413 Ok(CurvatureInference {
9414 term_idx,
9415 kappa_hat,
9416 ci,
9417 flatness,
9418 length_scale_hat: eta_hat.exp(),
9419 length_scale_estimated: profile.eta_bounds.is_some(),
9420 length_scale_support: range_outcome.support(),
9421 })
9422}
9423
9424#[cfg(test)]
9425mod curvature_profile_score_tests {
9426 use super::*;
9427
9428 #[test]
9429 fn analytic_profile_score_finds_exact_quadratic_lr_crossings() {
9430 let kappa_hat = -0.37;
9431 let curvature = 16.0;
9432 let level = 0.95;
9433 let mut profile = |kappa: f64| -> Result<(f64, f64), String> {
9434 let displacement = kappa - kappa_hat;
9435 Ok((
9436 7.0 + 0.5 * curvature * displacement * displacement,
9437 curvature * displacement,
9438 ))
9439 };
9440 let ci = curvature_profile_ci_from_analytic_score(
9441 &mut profile,
9442 kappa_hat,
9443 -3.0,
9444 3.0,
9445 level,
9446 1.0e-10,
9447 )
9448 .expect("analytic quadratic profile CI");
9449 let z = gam_geometry::curvature_estimand::wald_half_width(1.0, level)
9450 .expect("valid normal quantile");
9451 let expected_half_width = z / curvature.sqrt();
9452 assert!((ci.ci_lo - (kappa_hat - expected_half_width)).abs() <= 1.0e-8);
9453 assert!((ci.ci_hi - (kappa_hat + expected_half_width)).abs() <= 1.0e-8);
9454 assert!(!ci.lo_at_bound && !ci.hi_at_bound);
9455 }
9456
9457 #[test]
9458 fn analytic_profile_marks_chart_bound_when_wilks_set_never_crosses() {
9459 let mut profile =
9460 |kappa: f64| -> Result<(f64, f64), String> { Ok((0.5 * kappa * kappa, kappa)) };
9461 let ci =
9462 curvature_profile_ci_from_analytic_score(&mut profile, 0.0, -0.1, 0.1, 0.95, 1.0e-10)
9463 .expect("open bounded profile CI");
9464 assert_eq!(ci.ci_lo, -0.1);
9465 assert_eq!(ci.ci_hi, 0.1);
9466 assert!(ci.lo_at_bound && ci.hi_at_bound);
9467 assert_eq!(
9470 ci.kappa_hat_support,
9471 gam_geometry::curvature_estimand::KappaEstimateSupport::Interior
9472 );
9473 }
9474
9475 #[test]
9481 fn a_railed_point_estimate_is_accepted_and_declared_by_the_analytic_route_2687() {
9482 let kappa_max = 1.388_888_888_888_888_9_f64;
9485 let mut monotone = |kappa: f64| -> Result<(f64, f64), String> { Ok((-kappa, -1.0)) };
9486 let ci = curvature_profile_ci_from_analytic_score(
9487 &mut monotone,
9488 kappa_max,
9489 -kappa_max,
9490 kappa_max,
9491 0.95,
9492 1.0e-10,
9493 )
9494 .expect("a boundary optimum with the score pointing out of the box is stationary");
9495 assert_eq!(
9496 ci.kappa_hat_support,
9497 gam_geometry::curvature_estimand::KappaEstimateSupport::RailedAtUpperBound,
9498 "κ̂ = {kappa_max} is the box's own upper end"
9499 );
9500 let mut increasing = |kappa: f64| -> Result<(f64, f64), String> { Ok((kappa, 1.0)) };
9503 let ci_lo = curvature_profile_ci_from_analytic_score(
9504 &mut increasing,
9505 -kappa_max,
9506 -kappa_max,
9507 kappa_max,
9508 0.95,
9509 1.0e-10,
9510 )
9511 .expect("the mirrored boundary optimum");
9512 assert_eq!(
9513 ci_lo.kappa_hat_support,
9514 gam_geometry::curvature_estimand::KappaEstimateSupport::RailedAtLowerBound
9515 );
9516 let mut interior_slope = |kappa: f64| -> Result<(f64, f64), String> { Ok((-kappa, -1.0)) };
9519 assert!(
9520 curvature_profile_ci_from_analytic_score(
9521 &mut interior_slope,
9522 0.0,
9523 -kappa_max,
9524 kappa_max,
9525 0.95,
9526 1.0e-10,
9527 )
9528 .is_err(),
9529 "a non-stationary INTERIOR point is not an optimum and must still be refused"
9530 );
9531 }
9532}
9533
9534#[cfg(test)]
9535mod nfree_gate_tests {
9536 use super::nfree_skip_gate_status_from_parts;
9537
9538 #[test]
9539 fn value_only_nfree_gate_does_not_require_basis_skip_witness() {
9540 let gate = nfree_skip_gate_status_from_parts(
9541 true, true, false, false, true, true, false, false, );
9550 assert!(
9551 gate.would_skip(false),
9552 "value-only κ cost probes must stay n-free when the Gram value is certified; \
9553 the reduced-basis skip witness is required only for beta/gradient probes"
9554 );
9555 }
9556
9557 #[test]
9558 fn gradient_nfree_gate_still_requires_basis_skip_witness() {
9559 let gate =
9560 nfree_skip_gate_status_from_parts(true, true, false, true, true, true, false, true);
9561 assert!(
9562 !gate.would_skip(true),
9563 "gradient probes return beta/gradient objects in a reduced basis and must not \
9564 skip the row lane without the reduced-basis witness"
9565 );
9566 }
9567}