Skip to main content

gam_terms/smooth/
scale_contract.rs

1//! Exhaustive coordinate-scale contracts for smooth basis construction.
2//!
3//! A scale contract is not a human-readable tag.  It is a typed description of
4//! the pullback that must be applied to every coordinate-bearing part of a
5//! basis: design, penalty, coordinate derivatives, null geometry, and
6//! dimensionful hyperparameters.  [`SmoothBasisSpec::scale_contract`] matches
7//! every enum arm without a wildcard, so adding a new basis variant is a compile
8//! error until its law is declared here.
9
10use super::*;
11use crate::basis::{OneDimensionalBoundary, SphereMethod};
12
13/// Builder-level basis family.  Variants that share a `SmoothBasisSpec` arm but
14/// use different mathematics (cyclic/open/natural-cubic, Wahba/harmonic sphere,
15/// pure/hybrid Duchon, and factor-wrapper flavours) remain distinct here.
16#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
17pub enum BasisScaleFamily {
18    ByVariableNumeric,
19    ByVariableFactor,
20    FactorSumToZero,
21    OpenBSpline,
22    CyclicBSpline,
23    NaturalCubic,
24    BySmoothNumeric,
25    BySmoothFactor,
26    FactorSmoothFs,
27    FactorSmoothSz,
28    FactorSmoothRe,
29    ThinPlate,
30    SphereWahba,
31    SphereHarmonic,
32    ConstantCurvature,
33    Matern,
34    MeasureJet,
35    PureDuchon,
36    HybridDuchon,
37    Pca,
38    TensorBSpline,
39}
40
41impl BasisScaleFamily {
42    /// Canonical registry order.  The registry-completeness test walks this
43    /// array and compares it with a concrete `SmoothBasisSpec` zoo.
44    pub const ALL: [Self; 21] = [
45        Self::ByVariableNumeric,
46        Self::ByVariableFactor,
47        Self::FactorSumToZero,
48        Self::OpenBSpline,
49        Self::CyclicBSpline,
50        Self::NaturalCubic,
51        Self::BySmoothNumeric,
52        Self::BySmoothFactor,
53        Self::FactorSmoothFs,
54        Self::FactorSmoothSz,
55        Self::FactorSmoothRe,
56        Self::ThinPlate,
57        Self::SphereWahba,
58        Self::SphereHarmonic,
59        Self::ConstantCurvature,
60        Self::Matern,
61        Self::MeasureJet,
62        Self::PureDuchon,
63        Self::HybridDuchon,
64        Self::Pca,
65        Self::TensorBSpline,
66    ];
67}
68
69/// Coordinate action under which the declared law is exact.
70#[derive(Debug, Clone, Copy, PartialEq, Eq)]
71pub enum BasisCoordinateScaleAction {
72    /// `x' = a*x + b`, `a > 0`; all knots/endpoints move by the same map.
73    PositiveAffineAbscissa,
74    /// Each tensor margin has its own positive affine map.
75    IndependentPositiveAffineAxes,
76    /// `x' = a*x + b`, with one `a > 0` shared by every Euclidean axis;
77    /// centers, periodic lengths, and kernel range move with the same map.
78    /// (A diagonal anisotropy cannot represent a general rotated metric, so
79    /// rotation is deliberately not claimed by this scale-only contract.)
80    UniformEuclideanScale,
81    /// Degrees/radians are two coordinate encodings of the same point on S².
82    IntrinsicAngularUnitConversion,
83    /// `x'=a*x`, `kappa'=kappa/a^2`, `ell'=a*ell` in the stereographic chart.
84    ConstantCurvatureChartSimilarity,
85    /// `x'=a*x`, `mean'=a*mean`, `loadings'=loadings/a`, preserving PCA scores.
86    PcaScoreGauge,
87    /// Multiply an inner smooth by a numeric `by` coordinate.
88    NumericModulation,
89    /// Replicate/gate an inner smooth using dimensionless categorical labels.
90    DiscreteReplication,
91}
92
93/// Transformation of the emitted design matrix.
94#[derive(Debug, Clone, Copy, PartialEq, Eq)]
95pub enum BasisDesignScaleLaw {
96    /// The declared joint coordinate/parameter action leaves every design entry
97    /// invariant (up to floating-point roundoff).
98    Invariant,
99    /// Scaling only the numeric multiplier by `a` sends `X -> a*X`; scaling an
100    /// inner coordinate follows the child's law.
101    NumericMultiplierDegreeOne,
102    /// Wrapper design is a row gate/replication of its invariant inner design.
103    ReplicatedInner,
104    /// Tensor design is the row-wise product of invariant marginal designs.
105    TensorProductOfMarginals,
106    /// A random-intercept/slope block has alternating degree-zero and
107    /// degree-one columns in each factor level: intercept columns are invariant
108    /// and slope columns gain one power of the abscissa scale.
109    RandomInterceptSlopeDegreesZeroAndOne,
110    /// Every design entry gains exactly one power of the chart scale `a`.
111    ///
112    /// The constant-curvature smooth's kernel is `ℓ·(e^{−d_κ/ℓ} − 1)`, which is
113    /// a LENGTH: the `ℓ` factor is what stops the realized design collapsing
114    /// like `1/ℓ` and dragging `ρ̂` with it, and what makes the `ℓ → ∞` face the
115    /// finite geodesic-distance kernel rather than nothing (gam#2747,
116    /// `constant_curvature_kernel_matrix`). Under `x' = a·x`, `κ' = κ/a²`,
117    /// `ℓ' = a·ℓ` the exponent `d_κ/ℓ` is invariant and the prefactor is not, so
118    /// `X' = a·X` exactly.
119    ///
120    /// This is equivariance, not a broken invariance: the design gains a power
121    /// the penalty gains too, and a common positive factor on `(X, S)` is
122    /// absorbed exactly by `λ`, so the fitted function is unchanged. A basis
123    /// whose columns carry units is ordinary — a parametric `x` column does —
124    /// and the alternative (a dimensionless kernel) is what confounded the range
125    /// coordinate with the smoothing parameter.
126    ChartSimilarityDegreeOne,
127}
128
129/// Transformation of active penalty matrices.
130#[derive(Debug, Clone, Copy, PartialEq, Eq)]
131pub enum BasisPenaltyScaleLaw {
132    /// The raw functional has the stated coordinate homogeneity, while the
133    /// emitted unit-Frobenius matrix is invariant and its recorded
134    /// `normalization_scale` carries that power.
135    FrobeniusNormalizedRawPower(i32),
136    /// The raw derivative-energy functional has exact homogeneity
137    /// `a^(1 - 2*order)`; keeping `order` typed avoids truncating a structural
138    /// basis parameter into a fixed-width exponent.
139    FrobeniusNormalizedDerivativeOrder { order: usize },
140    /// Each tensor marginal carries its own `1 - 2*marginal_order` raw power;
141    /// every emitted normalized Kronecker penalty is invariant.
142    FrobeniusNormalizedPerMarginalDerivativeOrder,
143    /// Kernel/energy construction and its normalization are invariant under the
144    /// complete declared parameter pullback.
145    FrobeniusNormalizedInvariant,
146    /// The physical RKHS Gram is emitted without arbitrary normalization but
147    /// is itself invariant under the complete coordinate/parameter pullback.
148    InvariantPhysicalRkhsGram,
149    /// Penalties are copied/congruence-transformed into wrapper blocks without
150    /// introducing a coordinate scale of their own.
151    ReplicatedInner,
152    /// PCA uses the empirical function-mass Gram `X_score'X_score/n`; invariant
153    /// scores therefore give an exactly invariant penalty and null geometry.
154    InvariantFunctionMass,
155    /// The physical RKHS Gram is emitted without arbitrary normalization and
156    /// gains exactly one power of the chart scale, for the same reason its
157    /// design does (see
158    /// [`BasisDesignScaleLaw::ChartSimilarityDegreeOne`]): design and penalty
159    /// are the two blocks of ONE kernel, so they cannot scale differently
160    /// without the penalty ceasing to be the roughness of the design.
161    PhysicalRkhsGramDegreeOne,
162}
163
164/// Transformation of analytic input/hyperparameter derivatives.
165#[derive(Debug, Clone, Copy, PartialEq, Eq)]
166pub enum BasisDerivativeScaleLaw {
167    /// An order-q coordinate derivative transforms as `D'^q = a^-q D^q`.
168    InverseCoordinatePower { maximum_order: usize },
169    /// The same coordinate law, while derivatives in `log(kappa)=-log(ell)`
170    /// remain invariant under `ell -> a*ell`.
171    InverseCoordinatePowerAndInvariantLogRange { maximum_order: usize },
172    /// Under `kappa'=kappa/a^2`, first and second derivatives with respect to
173    /// the numeric curvature coordinate gain factors `a^3` and `a^5` — the
174    /// `a^2` and `a^4` from differentiating in a coordinate that itself scales
175    /// like `a^-2`, times the one power the kernel VALUE carries (see
176    /// [`BasisDesignScaleLaw::ChartSimilarityDegreeOne`]).
177    ConstantCurvatureParameterPowers,
178    /// Derivatives live in the intrinsic angular chart; changing degrees to
179    /// radians applies the ordinary inverse unit-conversion chain rule.
180    IntrinsicAngularChainRule,
181    /// Wrapper derivatives obey the product rule with the numeric multiplier.
182    NumericModulationProductRule,
183    /// Discrete gates/replications introduce no derivative coordinate.
184    DelegatedToInner,
185    /// Tensor partial derivatives take the inverse power on the differentiated
186    /// margin and multiply the invariant values of all other margins.
187    TensorMarginalProductRule { maximum_order: usize },
188}
189
190/// Transformation of the structural penalty null space / identifiability
191/// section.  This is explicit because equal penalty ranks alone do not prove
192/// that the same unpenalized functions survived a coordinate change.
193#[derive(Debug, Clone, Copy, PartialEq, Eq)]
194pub enum BasisNullGeometryScaleLaw {
195    PolynomialPullback,
196    CyclicConstantMode,
197    TensorProductPullback,
198    EuclideanPolynomialPullback,
199    CenterConstraintPullback,
200    MeasureJetAffineHeadPullback,
201    IntrinsicHarmonicSubspace,
202    ConstantCurvatureCenterConstraint,
203    PcaScoreCongruence,
204    FullRankRandomEffect,
205    ReplicatedInner,
206}
207
208/// Every dimensionful object that must move with a coordinate-unit change.
209#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
210pub enum DimensionfulBasisParameter {
211    Knots,
212    DomainEndpoints,
213    Periods,
214    Centers,
215    IsotropicInputScale,
216    LengthScale,
217    MeasureJetScaleBand,
218    MeasureJetCoordinateNoise,
219    Curvature,
220    PcaCenterMean,
221    PcaLoadings,
222    NumericByMultiplier,
223}
224
225/// Integer homogeneity of a dimensionful parameter under `x -> a*x`.
226/// `power=1` means multiply by `a`, `-1` divide by `a`, and `-2` divide by
227/// `a^2`.  Angular units use `power=1` with the degree/radian conversion in
228/// place of `a`.
229#[derive(Debug, Clone, Copy, PartialEq, Eq)]
230pub struct DimensionfulParameterScale {
231    pub parameter: DimensionfulBasisParameter,
232    pub power: i8,
233}
234
235impl DimensionfulParameterScale {
236    const fn new(parameter: DimensionfulBasisParameter, power: i8) -> Self {
237        Self { parameter, power }
238    }
239}
240
241/// How construction realizes the input frame.  This is deliberately private:
242/// callers consume it through [`BasisScaleContract::normalize_euclidean_frame`]
243/// rather than branching on a second public policy enum.
244#[derive(Debug, Clone, Copy, PartialEq, Eq)]
245enum InputFrameNormalization {
246    Parameterized,
247    AutoStandardizedOriginalUnits,
248    AutoStandardizedFreshOriginalReplayRealized,
249    Intrinsic,
250    PcaGauge,
251    Delegated,
252}
253
254/// Complete scale declaration for one node of a smooth-basis tree.
255#[derive(Debug, Clone, PartialEq, Eq)]
256pub struct BasisScaleContract {
257    pub family: BasisScaleFamily,
258    pub coordinate_action: BasisCoordinateScaleAction,
259    pub design: BasisDesignScaleLaw,
260    pub penalty: BasisPenaltyScaleLaw,
261    pub derivatives: BasisDerivativeScaleLaw,
262    pub null_geometry: BasisNullGeometryScaleLaw,
263    pub dimensionful_parameters: Vec<DimensionfulParameterScale>,
264    /// Inner/marginal contracts in construction order.  Wrappers have one
265    /// child; tensors have one child per marginal; leaf bases have none.
266    pub children: Vec<BasisScaleContract>,
267    input_frame: InputFrameNormalization,
268}
269
270/// Realized Euclidean input frame returned to the spatial constructors.
271pub(super) struct NormalizedEuclideanFrame {
272    pub coordinates: Array2<f64>,
273    pub input_scale: crate::IsotropicScale,
274    /// The realized range, in the same frame as `coordinates`.  The tag is
275    /// the contract's output: the private [`InputFrameNormalization`] decides
276    /// whether a conversion was applied, and the type carries that decision
277    /// forward instead of leaving readers to re-derive it by variant name.
278    pub length_scale: Option<crate::StandardizedUnits>,
279}
280
281fn standardize_resolved_center_strategy(
282    strategy: &mut crate::basis::CenterStrategy,
283    input_scale: crate::IsotropicScale,
284) {
285    match strategy {
286        crate::basis::CenterStrategy::Auto(inner) => {
287            standardize_resolved_center_strategy(inner, input_scale);
288        }
289        crate::basis::CenterStrategy::DuchonSpectral { knots, .. } => {
290            standardize_resolved_center_strategy(knots, input_scale);
291        }
292        crate::basis::CenterStrategy::UserProvided(centers) => {
293            input_scale.standardize(centers);
294        }
295        crate::basis::CenterStrategy::EqualMass { .. }
296        | crate::basis::CenterStrategy::EqualMassCovarRepresentative { .. }
297        | crate::basis::CenterStrategy::FarthestPoint { .. }
298        | crate::basis::CenterStrategy::KMeans { .. }
299        | crate::basis::CenterStrategy::UniformGrid { .. } => {}
300    }
301}
302
303impl BasisScaleContract {
304    fn leaf(
305        family: BasisScaleFamily,
306        coordinate_action: BasisCoordinateScaleAction,
307        design: BasisDesignScaleLaw,
308        penalty: BasisPenaltyScaleLaw,
309        derivatives: BasisDerivativeScaleLaw,
310        null_geometry: BasisNullGeometryScaleLaw,
311        dimensionful_parameters: Vec<DimensionfulParameterScale>,
312        input_frame: InputFrameNormalization,
313    ) -> Self {
314        Self {
315            family,
316            coordinate_action,
317            design,
318            penalty,
319            derivatives,
320            null_geometry,
321            dimensionful_parameters,
322            children: Vec::new(),
323            input_frame,
324        }
325    }
326
327    fn wrapper(
328        family: BasisScaleFamily,
329        coordinate_action: BasisCoordinateScaleAction,
330        design: BasisDesignScaleLaw,
331        derivatives: BasisDerivativeScaleLaw,
332        child: BasisScaleContract,
333        dimensionful_parameters: Vec<DimensionfulParameterScale>,
334    ) -> Self {
335        Self {
336            family,
337            coordinate_action,
338            design,
339            penalty: BasisPenaltyScaleLaw::ReplicatedInner,
340            derivatives,
341            null_geometry: BasisNullGeometryScaleLaw::ReplicatedInner,
342            dimensionful_parameters,
343            children: vec![child],
344            input_frame: InputFrameNormalization::Delegated,
345        }
346    }
347
348    /// Standardize a Euclidean spatial input using this family's declared
349    /// original-unit/replay law.  This is the sole construction path for the
350    /// ThinPlate, Matérn, Duchon, and MeasureJet `input_scale` fields.
351    ///
352    /// `spec_length_scale` is deliberately UNTAGGED: whether the spec declares
353    /// its range in original or realized units is exactly the question this
354    /// contract answers, and it answers it from `input_frame` plus whether
355    /// this is a fresh build or a replay.  The returned range is tagged,
356    /// because from there on the frame is settled (#2636).
357    pub(super) fn normalize_euclidean_frame(
358        &self,
359        mut coordinates: Array2<f64>,
360        stored_scale: Option<crate::IsotropicScale>,
361        spec_length_scale: Option<f64>,
362        center_strategy: &mut crate::basis::CenterStrategy,
363    ) -> Result<NormalizedEuclideanFrame, BasisError> {
364        let replay = stored_scale.is_some();
365        let replay_range_is_already_realized = match self.input_frame {
366            InputFrameNormalization::AutoStandardizedOriginalUnits => false,
367            InputFrameNormalization::AutoStandardizedFreshOriginalReplayRealized => true,
368            InputFrameNormalization::Parameterized
369            | InputFrameNormalization::Intrinsic
370            | InputFrameNormalization::PcaGauge
371            | InputFrameNormalization::Delegated => {
372                return Err(BasisError::InvalidInput(format!(
373                    "basis {:?} does not declare Euclidean auto-standardization",
374                    self.family
375                )));
376            }
377        };
378        let input_scale = match stored_scale {
379            Some(scale) => scale,
380            None => estimate_isotropic_scale(coordinates.view())?,
381        };
382
383        input_scale.standardize(&mut coordinates);
384        // A center matrix is a coordinate-bearing part of the kernel, not a
385        // dimensionless tuning artifact. Fresh user-provided centers and
386        // centers resolved by model materialization are expressed in the same
387        // original input frame as `coordinates`, so the frame transition must
388        // pull back both together. Leaving the centers untouched evaluates
389        // K(x / sigma, c) instead of K(x / sigma, c / sigma): translation and
390        // scale invariance are then genuinely broken, even though the error can
391        // look like harmless roundoff when sigma is close to one (#2623).
392        //
393        // Frozen replay is deliberately different: its center matrix was
394        // emitted by the builder in the realized standardized frame. The saved
395        // `stored_scale` marks that state, so scaling it again would double
396        // divide fit-time geometry at prediction.
397        if !replay {
398            standardize_resolved_center_strategy(center_strategy, input_scale);
399        }
400        let transformed_length = if replay_range_is_already_realized && replay {
401            // The frozen spec already holds the realized range; re-applying
402            // the conversion here is the double-divide the tag now forbids
403            // everywhere downstream.
404            spec_length_scale.map(crate::StandardizedUnits::new)
405        } else if replay_range_is_already_realized {
406            spec_length_scale.map(|ell| {
407                if ell > 0.0 {
408                    input_scale.to_standardized_units(crate::OriginalUnits::new(ell))
409                } else {
410                    // A non-positive sentinel is a mode marker, not a length;
411                    // it carries no units and must survive verbatim.
412                    crate::StandardizedUnits::new(ell)
413                }
414            })
415        } else {
416            spec_length_scale
417                .map(|ell| input_scale.to_standardized_units(crate::OriginalUnits::new(ell)))
418        };
419
420        Ok(NormalizedEuclideanFrame {
421            coordinates,
422            input_scale,
423            length_scale: transformed_length,
424        })
425    }
426}
427
428fn scale(power: i8, parameter: DimensionfulBasisParameter) -> DimensionfulParameterScale {
429    DimensionfulParameterScale::new(parameter, power)
430}
431
432fn bspline_contract(spec: &BSplineBasisSpec) -> BasisScaleContract {
433    let (family, null_geometry, mut parameters) = match (&spec.knotspec, &spec.boundary) {
434        (BSplineKnotSpec::NaturalCubicRegression { .. }, _) => (
435            BasisScaleFamily::NaturalCubic,
436            BasisNullGeometryScaleLaw::PolynomialPullback,
437            vec![scale(1, DimensionfulBasisParameter::Knots)],
438        ),
439        (BSplineKnotSpec::PeriodicUniform { .. }, _)
440        | (_, OneDimensionalBoundary::Cyclic { .. }) => (
441            BasisScaleFamily::CyclicBSpline,
442            BasisNullGeometryScaleLaw::CyclicConstantMode,
443            vec![
444                scale(1, DimensionfulBasisParameter::Knots),
445                scale(1, DimensionfulBasisParameter::DomainEndpoints),
446                scale(1, DimensionfulBasisParameter::Periods),
447            ],
448        ),
449        (BSplineKnotSpec::Generate { .. }, OneDimensionalBoundary::Open)
450        | (BSplineKnotSpec::Automatic { .. }, OneDimensionalBoundary::Open)
451        | (BSplineKnotSpec::Provided(_), OneDimensionalBoundary::Open) => (
452            BasisScaleFamily::OpenBSpline,
453            BasisNullGeometryScaleLaw::PolynomialPullback,
454            vec![scale(1, DimensionfulBasisParameter::Knots)],
455        ),
456    };
457    if matches!(spec.knotspec, BSplineKnotSpec::Generate { .. }) {
458        parameters.push(scale(1, DimensionfulBasisParameter::DomainEndpoints));
459    }
460    let penalty = if family == BasisScaleFamily::NaturalCubic {
461        BasisPenaltyScaleLaw::FrobeniusNormalizedRawPower(-3)
462    } else {
463        BasisPenaltyScaleLaw::FrobeniusNormalizedDerivativeOrder {
464            order: spec.penalty_order,
465        }
466    };
467    BasisScaleContract::leaf(
468        family,
469        BasisCoordinateScaleAction::PositiveAffineAbscissa,
470        BasisDesignScaleLaw::Invariant,
471        penalty,
472        BasisDerivativeScaleLaw::InverseCoordinatePower { maximum_order: 2 },
473        null_geometry,
474        parameters,
475        InputFrameNormalization::Parameterized,
476    )
477}
478
479fn spatial_parameters(
480    include_length: bool,
481    include_periods: bool,
482) -> Vec<DimensionfulParameterScale> {
483    let mut parameters = vec![
484        scale(1, DimensionfulBasisParameter::Centers),
485        scale(1, DimensionfulBasisParameter::IsotropicInputScale),
486    ];
487    if include_length {
488        parameters.push(scale(1, DimensionfulBasisParameter::LengthScale));
489    }
490    if include_periods {
491        parameters.push(scale(1, DimensionfulBasisParameter::Periods));
492    }
493    parameters
494}
495
496impl SmoothBasisSpec {
497    /// Return the complete scale contract for this concrete basis tree.
498    ///
499    /// There is intentionally no wildcard arm: a new `SmoothBasisSpec` variant
500    /// cannot compile until its design, penalty, derivative, null-space, and
501    /// dimensionful-parameter laws are all selected here.
502    pub fn scale_contract(&self) -> BasisScaleContract {
503        match self {
504            SmoothBasisSpec::ByVariable { inner, by, .. } => match by {
505                ByVariableSpec::Numeric => BasisScaleContract::wrapper(
506                    BasisScaleFamily::ByVariableNumeric,
507                    BasisCoordinateScaleAction::NumericModulation,
508                    BasisDesignScaleLaw::NumericMultiplierDegreeOne,
509                    BasisDerivativeScaleLaw::NumericModulationProductRule,
510                    inner.scale_contract(),
511                    vec![scale(1, DimensionfulBasisParameter::NumericByMultiplier)],
512                ),
513                ByVariableSpec::Level { .. } => BasisScaleContract::wrapper(
514                    BasisScaleFamily::ByVariableFactor,
515                    BasisCoordinateScaleAction::DiscreteReplication,
516                    BasisDesignScaleLaw::ReplicatedInner,
517                    BasisDerivativeScaleLaw::DelegatedToInner,
518                    inner.scale_contract(),
519                    Vec::new(),
520                ),
521            },
522            SmoothBasisSpec::FactorSumToZero { inner, .. } => BasisScaleContract::wrapper(
523                BasisScaleFamily::FactorSumToZero,
524                BasisCoordinateScaleAction::DiscreteReplication,
525                BasisDesignScaleLaw::ReplicatedInner,
526                BasisDerivativeScaleLaw::DelegatedToInner,
527                inner.scale_contract(),
528                Vec::new(),
529            ),
530            SmoothBasisSpec::BSpline1D { spec, .. } => bspline_contract(spec),
531            SmoothBasisSpec::BySmooth { smooth, by_kind } => match by_kind {
532                ByVarKind::Numeric { .. } => BasisScaleContract::wrapper(
533                    BasisScaleFamily::BySmoothNumeric,
534                    BasisCoordinateScaleAction::NumericModulation,
535                    BasisDesignScaleLaw::NumericMultiplierDegreeOne,
536                    BasisDerivativeScaleLaw::NumericModulationProductRule,
537                    smooth.scale_contract(),
538                    vec![scale(1, DimensionfulBasisParameter::NumericByMultiplier)],
539                ),
540                ByVarKind::Factor { .. } => BasisScaleContract::wrapper(
541                    BasisScaleFamily::BySmoothFactor,
542                    BasisCoordinateScaleAction::DiscreteReplication,
543                    BasisDesignScaleLaw::ReplicatedInner,
544                    BasisDerivativeScaleLaw::DelegatedToInner,
545                    smooth.scale_contract(),
546                    Vec::new(),
547                ),
548            },
549            SmoothBasisSpec::FactorSmooth { spec } => match &spec.flavour {
550                FactorSmoothFlavour::Fs { .. } => BasisScaleContract::wrapper(
551                    BasisScaleFamily::FactorSmoothFs,
552                    BasisCoordinateScaleAction::DiscreteReplication,
553                    BasisDesignScaleLaw::ReplicatedInner,
554                    BasisDerivativeScaleLaw::DelegatedToInner,
555                    bspline_contract(&spec.marginal),
556                    Vec::new(),
557                ),
558                FactorSmoothFlavour::Sz => BasisScaleContract::wrapper(
559                    BasisScaleFamily::FactorSmoothSz,
560                    BasisCoordinateScaleAction::DiscreteReplication,
561                    BasisDesignScaleLaw::ReplicatedInner,
562                    BasisDerivativeScaleLaw::DelegatedToInner,
563                    bspline_contract(&spec.marginal),
564                    Vec::new(),
565                ),
566                FactorSmoothFlavour::Re => BasisScaleContract {
567                    family: BasisScaleFamily::FactorSmoothRe,
568                    coordinate_action: BasisCoordinateScaleAction::PositiveAffineAbscissa,
569                    design: BasisDesignScaleLaw::RandomInterceptSlopeDegreesZeroAndOne,
570                    penalty: BasisPenaltyScaleLaw::FrobeniusNormalizedInvariant,
571                    derivatives: BasisDerivativeScaleLaw::InverseCoordinatePower {
572                        maximum_order: 1,
573                    },
574                    null_geometry: BasisNullGeometryScaleLaw::FullRankRandomEffect,
575                    dimensionful_parameters: vec![scale(1, DimensionfulBasisParameter::Knots)],
576                    children: vec![bspline_contract(&spec.marginal)],
577                    input_frame: InputFrameNormalization::Delegated,
578                },
579            },
580            SmoothBasisSpec::ThinPlate { .. } => BasisScaleContract::leaf(
581                BasisScaleFamily::ThinPlate,
582                BasisCoordinateScaleAction::UniformEuclideanScale,
583                BasisDesignScaleLaw::Invariant,
584                BasisPenaltyScaleLaw::FrobeniusNormalizedInvariant,
585                BasisDerivativeScaleLaw::InverseCoordinatePowerAndInvariantLogRange {
586                    maximum_order: 2,
587                },
588                BasisNullGeometryScaleLaw::EuclideanPolynomialPullback,
589                spatial_parameters(true, true),
590                InputFrameNormalization::AutoStandardizedOriginalUnits,
591            ),
592            SmoothBasisSpec::Sphere { spec, .. } => {
593                let (family, parameters) = match spec.method {
594                    SphereMethod::Wahba => (
595                        BasisScaleFamily::SphereWahba,
596                        vec![scale(1, DimensionfulBasisParameter::Centers)],
597                    ),
598                    SphereMethod::Harmonic => (BasisScaleFamily::SphereHarmonic, Vec::new()),
599                };
600                BasisScaleContract::leaf(
601                    family,
602                    BasisCoordinateScaleAction::IntrinsicAngularUnitConversion,
603                    BasisDesignScaleLaw::Invariant,
604                    BasisPenaltyScaleLaw::FrobeniusNormalizedInvariant,
605                    BasisDerivativeScaleLaw::IntrinsicAngularChainRule,
606                    BasisNullGeometryScaleLaw::IntrinsicHarmonicSubspace,
607                    parameters,
608                    InputFrameNormalization::Intrinsic,
609                )
610            }
611            SmoothBasisSpec::ConstantCurvature { .. } => BasisScaleContract::leaf(
612                BasisScaleFamily::ConstantCurvature,
613                BasisCoordinateScaleAction::ConstantCurvatureChartSimilarity,
614                BasisDesignScaleLaw::ChartSimilarityDegreeOne,
615                BasisPenaltyScaleLaw::PhysicalRkhsGramDegreeOne,
616                BasisDerivativeScaleLaw::ConstantCurvatureParameterPowers,
617                BasisNullGeometryScaleLaw::ConstantCurvatureCenterConstraint,
618                vec![
619                    scale(1, DimensionfulBasisParameter::Centers),
620                    scale(1, DimensionfulBasisParameter::LengthScale),
621                    scale(-2, DimensionfulBasisParameter::Curvature),
622                ],
623                InputFrameNormalization::Intrinsic,
624            ),
625            SmoothBasisSpec::Matern { .. } => BasisScaleContract::leaf(
626                BasisScaleFamily::Matern,
627                BasisCoordinateScaleAction::UniformEuclideanScale,
628                BasisDesignScaleLaw::Invariant,
629                BasisPenaltyScaleLaw::FrobeniusNormalizedInvariant,
630                BasisDerivativeScaleLaw::InverseCoordinatePowerAndInvariantLogRange {
631                    maximum_order: 2,
632                },
633                BasisNullGeometryScaleLaw::CenterConstraintPullback,
634                spatial_parameters(true, true),
635                InputFrameNormalization::AutoStandardizedOriginalUnits,
636            ),
637            SmoothBasisSpec::MeasureJet { .. } => BasisScaleContract::leaf(
638                BasisScaleFamily::MeasureJet,
639                BasisCoordinateScaleAction::UniformEuclideanScale,
640                BasisDesignScaleLaw::Invariant,
641                BasisPenaltyScaleLaw::FrobeniusNormalizedInvariant,
642                BasisDerivativeScaleLaw::InverseCoordinatePowerAndInvariantLogRange {
643                    maximum_order: 2,
644                },
645                BasisNullGeometryScaleLaw::MeasureJetAffineHeadPullback,
646                {
647                    let mut parameters = spatial_parameters(true, false);
648                    parameters.push(scale(1, DimensionfulBasisParameter::MeasureJetScaleBand));
649                    parameters.push(scale(
650                        1,
651                        DimensionfulBasisParameter::MeasureJetCoordinateNoise,
652                    ));
653                    parameters
654                },
655                InputFrameNormalization::AutoStandardizedFreshOriginalReplayRealized,
656            ),
657            SmoothBasisSpec::Duchon { spec, .. } => {
658                let family = if spec.length_scale.is_some() {
659                    BasisScaleFamily::HybridDuchon
660                } else {
661                    BasisScaleFamily::PureDuchon
662                };
663                BasisScaleContract::leaf(
664                    family,
665                    BasisCoordinateScaleAction::UniformEuclideanScale,
666                    BasisDesignScaleLaw::Invariant,
667                    BasisPenaltyScaleLaw::FrobeniusNormalizedInvariant,
668                    if spec.length_scale.is_some() {
669                        BasisDerivativeScaleLaw::InverseCoordinatePowerAndInvariantLogRange {
670                            maximum_order: 2,
671                        }
672                    } else {
673                        BasisDerivativeScaleLaw::InverseCoordinatePower { maximum_order: 2 }
674                    },
675                    BasisNullGeometryScaleLaw::EuclideanPolynomialPullback,
676                    spatial_parameters(spec.length_scale.is_some(), true),
677                    InputFrameNormalization::AutoStandardizedOriginalUnits,
678                )
679            }
680            SmoothBasisSpec::Pca { .. } => BasisScaleContract::leaf(
681                BasisScaleFamily::Pca,
682                BasisCoordinateScaleAction::PcaScoreGauge,
683                BasisDesignScaleLaw::Invariant,
684                BasisPenaltyScaleLaw::InvariantFunctionMass,
685                BasisDerivativeScaleLaw::InverseCoordinatePower { maximum_order: 1 },
686                BasisNullGeometryScaleLaw::PcaScoreCongruence,
687                vec![
688                    scale(1, DimensionfulBasisParameter::PcaCenterMean),
689                    scale(-1, DimensionfulBasisParameter::PcaLoadings),
690                ],
691                InputFrameNormalization::PcaGauge,
692            ),
693            SmoothBasisSpec::TensorBSpline { spec, .. } => BasisScaleContract {
694                family: BasisScaleFamily::TensorBSpline,
695                coordinate_action: BasisCoordinateScaleAction::IndependentPositiveAffineAxes,
696                design: BasisDesignScaleLaw::TensorProductOfMarginals,
697                penalty: BasisPenaltyScaleLaw::FrobeniusNormalizedPerMarginalDerivativeOrder,
698                derivatives: BasisDerivativeScaleLaw::TensorMarginalProductRule {
699                    maximum_order: 2,
700                },
701                null_geometry: BasisNullGeometryScaleLaw::TensorProductPullback,
702                dimensionful_parameters: Vec::new(),
703                children: spec.marginalspecs.iter().map(bspline_contract).collect(),
704                input_frame: InputFrameNormalization::Delegated,
705            },
706        }
707    }
708
709    /// Validate every scale-bearing field before construction or frozen replay.
710    /// Wrapper recursion is exhaustive; no basis can bypass this check.
711    pub fn validate_scale_configuration(&self) -> Result<(), BasisError> {
712        let contract = self.scale_contract();
713        match self {
714            SmoothBasisSpec::ByVariable { inner, .. }
715            | SmoothBasisSpec::FactorSumToZero { inner, .. } => {
716                inner.validate_scale_configuration()
717            }
718            SmoothBasisSpec::BSpline1D { .. } => Ok(()),
719            SmoothBasisSpec::BySmooth { smooth, .. } => smooth.validate_scale_configuration(),
720            SmoothBasisSpec::FactorSmooth { .. } => Ok(()),
721            SmoothBasisSpec::ThinPlate {
722                feature_cols,
723                input_scale: _,
724                ..
725            }
726            | SmoothBasisSpec::Matern {
727                feature_cols,
728                input_scale: _,
729                ..
730            }
731            | SmoothBasisSpec::MeasureJet {
732                feature_cols,
733                input_scale: _,
734                ..
735            }
736            | SmoothBasisSpec::Duchon {
737                feature_cols,
738                input_scale: _,
739                ..
740            } => {
741                if feature_cols.is_empty() {
742                    return Err(BasisError::InvalidInput(format!(
743                        "basis {:?} requires at least one coordinate axis",
744                        contract.family
745                    )));
746                }
747                Ok(())
748            }
749            SmoothBasisSpec::Sphere { feature_cols, .. }
750            | SmoothBasisSpec::ConstantCurvature { feature_cols, .. } => {
751                if feature_cols.is_empty() {
752                    return Err(BasisError::InvalidInput(format!(
753                        "basis {:?} requires at least one coordinate axis",
754                        contract.family
755                    )));
756                }
757                Ok(())
758            }
759            SmoothBasisSpec::Pca { .. } => Ok(()),
760            SmoothBasisSpec::TensorBSpline { feature_cols, spec } => {
761                if feature_cols.len() != spec.marginalspecs.len() {
762                    return Err(BasisError::DimensionMismatch(format!(
763                        "TensorBSpline has {} feature axes but {} marginal scale contracts",
764                        feature_cols.len(),
765                        spec.marginalspecs.len()
766                    )));
767                }
768                Ok(())
769            }
770        }
771    }
772}
773
774#[cfg(test)]
775mod tests {
776    use super::*;
777    use crate::basis::{
778        BSplineBoundaryConditions, BasisWorkspace, ConstantCurvatureIdentifiability,
779        DuchonNullspaceOrder, DuchonOperatorPenaltySpec, MaternIdentifiability, MaternNu,
780        MeasureJetIdentifiability, SphereWahbaKernel, SphericalSplineIdentifiability,
781        build_constant_curvature_basis, build_spherical_spline_basis,
782        constant_curvature_kernel_kappa_jets,
783    };
784    use ndarray::{Array1, Array2, array};
785    use std::collections::HashSet;
786
787    fn open_marginal(scale: f64) -> BSplineBasisSpec {
788        BSplineBasisSpec {
789            degree: 3,
790            penalty_order: 2,
791            knotspec: BSplineKnotSpec::Provided(Array1::from(vec![
792                0.0,
793                0.0,
794                0.0,
795                0.0,
796                0.25 * scale,
797                0.6 * scale,
798                scale,
799                scale,
800                scale,
801                scale,
802            ])),
803            double_penalty: false,
804            identifiability: BSplineIdentifiability::None,
805            boundary: OneDimensionalBoundary::Open,
806            boundary_conditions: BSplineBoundaryConditions::default(),
807        }
808    }
809
810    fn basis(feature_col: usize, scale: f64) -> SmoothBasisSpec {
811        SmoothBasisSpec::BSpline1D {
812            feature_col,
813            spec: open_marginal(scale),
814        }
815    }
816
817    fn factor_spec(flavour: FactorSmoothFlavour) -> SmoothBasisSpec {
818        SmoothBasisSpec::FactorSmooth {
819            spec: FactorSmoothSpec {
820                continuous_cols: vec![0],
821                group_col: 1,
822                marginal: open_marginal(1.0),
823                flavour,
824                group_frozen_levels: Some(vec![0.0_f64.to_bits(), 1.0_f64.to_bits()]),
825                frozen_global_orthogonality: None,
826            },
827        }
828    }
829
830    fn scale_contract_zoo() -> Vec<SmoothBasisSpec> {
831        let by_level_bits = 1.0_f64.to_bits();
832        let mut sphere_harmonic = SphericalSplineBasisSpec::default();
833        sphere_harmonic.method = SphereMethod::Harmonic;
834        sphere_harmonic.max_degree = Some(3);
835        vec![
836            SmoothBasisSpec::ByVariable {
837                inner: Box::new(basis(0, 1.0)),
838                by_col: 1,
839                kind: BySmoothKind::Numeric,
840                by: ByVariableSpec::Numeric,
841            },
842            SmoothBasisSpec::ByVariable {
843                inner: Box::new(basis(0, 1.0)),
844                by_col: 1,
845                kind: BySmoothKind::Level {
846                    level_bits: by_level_bits,
847                },
848                by: ByVariableSpec::Level {
849                    value_bits: by_level_bits,
850                    label: "one".to_string(),
851                },
852            },
853            SmoothBasisSpec::FactorSumToZero {
854                inner: Box::new(basis(0, 1.0)),
855                by_col: 1,
856                levels: vec![0.0_f64.to_bits(), by_level_bits],
857                frozen_global_orthogonality: None,
858            },
859            basis(0, 1.0),
860            SmoothBasisSpec::BSpline1D {
861                feature_col: 0,
862                spec: BSplineBasisSpec {
863                    knotspec: BSplineKnotSpec::PeriodicUniform {
864                        data_range: (0.0, 1.0),
865                        num_basis: 8,
866                    },
867                    boundary: OneDimensionalBoundary::Cyclic {
868                        start: 0.0,
869                        end: 1.0,
870                    },
871                    ..open_marginal(1.0)
872                },
873            },
874            SmoothBasisSpec::BSpline1D {
875                feature_col: 0,
876                spec: BSplineBasisSpec {
877                    knotspec: BSplineKnotSpec::NaturalCubicRegression {
878                        knots: array![0.0, 0.2, 0.5, 0.8, 1.0],
879                    },
880                    ..open_marginal(1.0)
881                },
882            },
883            SmoothBasisSpec::BySmooth {
884                smooth: Box::new(basis(0, 1.0)),
885                by_kind: ByVarKind::Numeric { feature_col: 1 },
886            },
887            SmoothBasisSpec::BySmooth {
888                smooth: Box::new(basis(0, 1.0)),
889                by_kind: ByVarKind::Factor {
890                    feature_col: 1,
891                    ordered: false,
892                    frozen_levels: Some(vec![0.0_f64.to_bits(), by_level_bits]),
893                },
894            },
895            factor_spec(FactorSmoothFlavour::Fs {
896                m_null_penalty_orders: vec![1],
897            }),
898            factor_spec(FactorSmoothFlavour::Sz),
899            factor_spec(FactorSmoothFlavour::Re),
900            SmoothBasisSpec::ThinPlate {
901                feature_cols: vec![0, 1],
902                spec: ThinPlateBasisSpec {
903                    center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
904                    periodic: None,
905                    length_scale: 0.7,
906                    double_penalty: false,
907                    identifiability: SpatialIdentifiability::None,
908                    radial_reparam: None,
909                },
910                input_scale: None,
911            },
912            SmoothBasisSpec::Sphere {
913                feature_cols: vec![0, 1],
914                spec: SphericalSplineBasisSpec::default(),
915            },
916            SmoothBasisSpec::Sphere {
917                feature_cols: vec![0, 1],
918                spec: sphere_harmonic,
919            },
920            SmoothBasisSpec::ConstantCurvature {
921                feature_cols: vec![0, 1],
922                spec: ConstantCurvatureBasisSpec::default(),
923            },
924            SmoothBasisSpec::Matern {
925                feature_cols: vec![0, 1],
926                spec: MaternBasisSpec {
927                    center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
928                    periodic: None,
929                    length_scale: crate::basis::MaternLengthScale::fixed(0.7),
930                    nu: MaternNu::ThreeHalves,
931                    include_intercept: false,
932                    double_penalty: false,
933                    identifiability: MaternIdentifiability::None,
934                    aniso_log_scales: None,
935                },
936                input_scale: None,
937            },
938            SmoothBasisSpec::MeasureJet {
939                feature_cols: vec![0, 1],
940                spec: MeasureJetBasisSpec::default(),
941                input_scale: None,
942            },
943            SmoothBasisSpec::Duchon {
944                feature_cols: vec![0, 1],
945                spec: DuchonBasisSpec {
946                    center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
947                    periodic: None,
948                    length_scale: None,
949                    power: 0.0,
950                    nullspace_order: DuchonNullspaceOrder::Linear,
951                    identifiability: SpatialIdentifiability::None,
952                    aniso_log_scales: None,
953                    operator_penalties: DuchonOperatorPenaltySpec::default(),
954                    boundary: OneDimensionalBoundary::Open,
955                    radial_reparam: None,
956                },
957                input_scale: None,
958            },
959            SmoothBasisSpec::Duchon {
960                feature_cols: vec![0, 1],
961                spec: DuchonBasisSpec {
962                    center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
963                    periodic: None,
964                    length_scale: Some(0.7),
965                    power: 1.0,
966                    nullspace_order: DuchonNullspaceOrder::Linear,
967                    identifiability: SpatialIdentifiability::None,
968                    aniso_log_scales: None,
969                    operator_penalties: DuchonOperatorPenaltySpec::default(),
970                    boundary: OneDimensionalBoundary::Open,
971                    radial_reparam: None,
972                },
973                input_scale: None,
974            },
975            SmoothBasisSpec::Pca {
976                feature_cols: vec![0, 1],
977                basis_matrix: array![[1.0], [0.0]],
978                centered: true,
979                smooth_penalty: 1.0,
980                center_mean: Some(array![0.0, 0.0]),
981                pca_basis_path: None,
982                chunk_size: 32,
983            },
984            SmoothBasisSpec::TensorBSpline {
985                feature_cols: vec![0, 1],
986                spec: TensorBSplineSpec {
987                    marginalspecs: vec![open_marginal(1.0), open_marginal(1.0)],
988                    periods: vec![None, None],
989                    double_penalty: false,
990                    identifiability: TensorBSplineIdentifiability::None,
991                    penalty_decomposition: TensorBSplinePenaltyDecomposition::MarginalKroneckerSum,
992                },
993            },
994        ]
995    }
996
997    fn zoo_basis(family: BasisScaleFamily) -> SmoothBasisSpec {
998        scale_contract_zoo()
999            .into_iter()
1000            .find(|basis| basis.scale_contract().family == family)
1001            .unwrap_or_else(|| panic!("scale-contract zoo is missing {family:?}"))
1002    }
1003
1004    #[test]
1005    fn scale_contract_registry_is_exhaustive_unique_and_typed_2315() {
1006        let zoo = scale_contract_zoo();
1007        let observed: HashSet<_> = zoo
1008            .iter()
1009            .map(|basis| basis.scale_contract().family)
1010            .collect();
1011        let expected: HashSet<_> = BasisScaleFamily::ALL.into_iter().collect();
1012        assert_eq!(zoo.len(), BasisScaleFamily::ALL.len());
1013        assert_eq!(observed, expected);
1014
1015        for basis in &zoo {
1016            let contract = basis.scale_contract();
1017            match contract.family {
1018                BasisScaleFamily::ByVariableNumeric
1019                | BasisScaleFamily::ByVariableFactor
1020                | BasisScaleFamily::FactorSumToZero
1021                | BasisScaleFamily::BySmoothNumeric
1022                | BasisScaleFamily::BySmoothFactor
1023                | BasisScaleFamily::FactorSmoothFs
1024                | BasisScaleFamily::FactorSmoothSz
1025                | BasisScaleFamily::FactorSmoothRe => {
1026                    assert_eq!(contract.children.len(), 1, "{:?}", contract.family);
1027                }
1028                BasisScaleFamily::TensorBSpline => {
1029                    assert_eq!(contract.children.len(), 2);
1030                }
1031                BasisScaleFamily::OpenBSpline
1032                | BasisScaleFamily::CyclicBSpline
1033                | BasisScaleFamily::NaturalCubic
1034                | BasisScaleFamily::ThinPlate
1035                | BasisScaleFamily::SphereWahba
1036                | BasisScaleFamily::SphereHarmonic
1037                | BasisScaleFamily::ConstantCurvature
1038                | BasisScaleFamily::Matern
1039                | BasisScaleFamily::MeasureJet
1040                | BasisScaleFamily::PureDuchon
1041                | BasisScaleFamily::HybridDuchon
1042                | BasisScaleFamily::Pca => assert!(contract.children.is_empty()),
1043            }
1044            let unique_parameters: HashSet<_> = contract
1045                .dimensionful_parameters
1046                .iter()
1047                .map(|parameter| parameter.parameter)
1048                .collect();
1049            assert_eq!(
1050                unique_parameters.len(),
1051                contract.dimensionful_parameters.len(),
1052                "duplicate dimensionful parameter in {:?}",
1053                contract.family
1054            );
1055        }
1056    }
1057
1058    fn assert_matrix_close(actual: &Array2<f64>, expected: &Array2<f64>, tolerance: f64) {
1059        assert_matrix_close_with_context(actual, expected, tolerance, "matrix");
1060    }
1061
1062    fn assert_matrix_close_with_context(
1063        actual: &Array2<f64>,
1064        expected: &Array2<f64>,
1065        tolerance: f64,
1066        context: &str,
1067    ) {
1068        assert_eq!(actual.dim(), expected.dim());
1069        for ((row, col), &target) in expected.indexed_iter() {
1070            let observed = actual[[row, col]];
1071            assert!(
1072                (observed - target).abs() <= tolerance * (1.0 + target.abs()),
1073                "{context}: matrix[{row},{col}] differs: observed={observed:.16e}, target={target:.16e}; actual={actual:?}; expected={expected:?}"
1074            );
1075        }
1076    }
1077
1078    fn assert_build_geometry_close(
1079        actual: &BasisBuildResult,
1080        expected: &BasisBuildResult,
1081        tolerance: f64,
1082    ) {
1083        assert_matrix_close(
1084            &actual.design.to_dense(),
1085            &expected.design.to_dense(),
1086            tolerance,
1087        );
1088        assert_eq!(
1089            actual.active_penalties.len(),
1090            expected.active_penalties.len()
1091        );
1092        for (observed, target) in actual
1093            .active_penalties
1094            .iter()
1095            .zip(expected.active_penalties.iter())
1096        {
1097            assert_eq!(observed.info.source, target.info.source);
1098            assert_eq!(observed.info.effective_rank, target.info.effective_rank);
1099            assert_eq!(observed.nullity, target.nullity);
1100            assert_matrix_close(&observed.matrix, &target.matrix, tolerance);
1101        }
1102    }
1103
1104    /// [`assert_build_geometry_close`] for a basis that is EQUIVARIANT with
1105    /// weight one rather than invariant: one common `scale` must restore both
1106    /// the design and every active penalty. Taking the same factor out of both
1107    /// is what distinguishes a units change from a model change.
1108    fn assert_build_geometry_scaled_close(
1109        actual: &BasisBuildResult,
1110        expected: &BasisBuildResult,
1111        scale: f64,
1112        tolerance: f64,
1113    ) {
1114        assert_matrix_close(
1115            &actual.design.to_dense().mapv(|value| value / scale),
1116            &expected.design.to_dense(),
1117            tolerance,
1118        );
1119        assert_eq!(
1120            actual.active_penalties.len(),
1121            expected.active_penalties.len()
1122        );
1123        for (observed, target) in actual
1124            .active_penalties
1125            .iter()
1126            .zip(expected.active_penalties.iter())
1127        {
1128            assert_eq!(observed.info.source, target.info.source);
1129            assert_eq!(observed.info.effective_rank, target.info.effective_rank);
1130            assert_eq!(observed.nullity, target.nullity);
1131            assert_matrix_close(
1132                &observed.matrix.mapv(|value| value / scale),
1133                &target.matrix,
1134                tolerance,
1135            );
1136        }
1137    }
1138
1139    fn assert_local_geometry_close(
1140        actual: &LocalSmoothTermBuild,
1141        expected: &LocalSmoothTermBuild,
1142        tolerance: f64,
1143        context: &str,
1144    ) {
1145        assert_matrix_close_with_context(
1146            &actual.design.to_dense(),
1147            &expected.design.to_dense(),
1148            tolerance,
1149            &format!("{context} design"),
1150        );
1151        assert_eq!(
1152            actual.active_penalties.len(),
1153            expected.active_penalties.len()
1154        );
1155        for (observed, target) in actual
1156            .active_penalties
1157            .iter()
1158            .zip(expected.active_penalties.iter())
1159        {
1160            assert_eq!(observed.info.source, target.info.source);
1161            assert_eq!(observed.info.effective_rank, target.info.effective_rank);
1162            assert_eq!(observed.nullity, target.nullity);
1163            assert_matrix_close_with_context(
1164                &observed.matrix,
1165                &target.matrix,
1166                tolerance,
1167                &format!("{context} penalty source {:?}", target.info.source),
1168            );
1169        }
1170    }
1171
1172    fn scaled_cyclic_marginal(scale: f64) -> BSplineBasisSpec {
1173        BSplineBasisSpec {
1174            degree: 3,
1175            penalty_order: 2,
1176            knotspec: BSplineKnotSpec::PeriodicUniform {
1177                data_range: (-0.4 * scale, 1.6 * scale),
1178                num_basis: 9,
1179            },
1180            double_penalty: false,
1181            identifiability: BSplineIdentifiability::None,
1182            boundary: OneDimensionalBoundary::Cyclic {
1183                start: -0.4 * scale,
1184                end: 1.6 * scale,
1185            },
1186            boundary_conditions: BSplineBoundaryConditions::default(),
1187        }
1188    }
1189
1190    fn scaled_natural_marginal(scale: f64) -> BSplineBasisSpec {
1191        BSplineBasisSpec {
1192            degree: 3,
1193            penalty_order: 2,
1194            knotspec: BSplineKnotSpec::NaturalCubicRegression {
1195                knots: array![0.0, 0.17, 0.43, 0.71, 1.0].mapv(|value| value * scale),
1196            },
1197            double_penalty: false,
1198            identifiability: BSplineIdentifiability::None,
1199            boundary: OneDimensionalBoundary::Open,
1200            boundary_conditions: BSplineBoundaryConditions::default(),
1201        }
1202    }
1203
1204    #[test]
1205    fn declared_spline_builders_obey_design_penalty_and_null_scale_laws_2315() {
1206        let open_points = array![0.0, 0.05, 0.3, 0.55, 0.72, 0.95, 1.0];
1207        let cyclic_points = array![-2.4, -0.4, -0.13, 0.2, 1.1, 1.6, 3.6];
1208        let natural_points = array![-0.35, 0.0, 0.09, 0.43, 0.86, 1.0, 1.28];
1209
1210        for (family, points, build_spec) in [
1211            (
1212                BasisScaleFamily::OpenBSpline,
1213                open_points,
1214                open_marginal as fn(f64) -> BSplineBasisSpec,
1215            ),
1216            (
1217                BasisScaleFamily::CyclicBSpline,
1218                cyclic_points,
1219                scaled_cyclic_marginal as fn(f64) -> BSplineBasisSpec,
1220            ),
1221            (
1222                BasisScaleFamily::NaturalCubic,
1223                natural_points,
1224                scaled_natural_marginal as fn(f64) -> BSplineBasisSpec,
1225            ),
1226        ] {
1227            let reference_spec = build_spec(1.0);
1228            assert_eq!(bspline_contract(&reference_spec).family, family);
1229            let reference = build_bspline_basis_1d(points.view(), &reference_spec)
1230                .expect("reference scalar basis");
1231            for factor in [1e-9_f64, 1.0, 1e9] {
1232                let actual = build_bspline_basis_1d(
1233                    points.mapv(|value| value * factor).view(),
1234                    &build_spec(factor),
1235                )
1236                .expect("rescaled scalar basis");
1237                assert_build_geometry_close(&actual, &reference, 8e-10);
1238                for (observed, target) in actual
1239                    .active_penalties
1240                    .iter()
1241                    .zip(reference.active_penalties.iter())
1242                {
1243                    let rescaled = observed.info.normalization_scale * factor.powi(3);
1244                    assert!(
1245                        (rescaled - target.info.normalization_scale).abs()
1246                            <= 2e-8 * (1.0 + target.info.normalization_scale.abs()),
1247                        "{family:?} raw penalty normalizer violated a^-3 at factor {factor}"
1248                    );
1249                }
1250            }
1251        }
1252
1253        let tensor_data = array![
1254            [0.00, 0.13],
1255            [0.08, 0.91],
1256            [0.22, 0.37],
1257            [0.41, 1.00],
1258            [0.58, 0.02],
1259            [0.73, 0.66],
1260            [0.89, 0.48],
1261            [1.00, 0.00]
1262        ];
1263        let build_tensor = |x_scale: f64, y_scale: f64| {
1264            let mut scaled = tensor_data.clone();
1265            scaled.column_mut(0).mapv_inplace(|value| value * x_scale);
1266            scaled.column_mut(1).mapv_inplace(|value| value * y_scale);
1267            let spec = TensorBSplineSpec {
1268                marginalspecs: vec![open_marginal(x_scale), open_marginal(y_scale)],
1269                periods: vec![None, None],
1270                double_penalty: false,
1271                identifiability: TensorBSplineIdentifiability::None,
1272                penalty_decomposition: TensorBSplinePenaltyDecomposition::MarginalKroneckerSum,
1273            };
1274            build_tensor_bspline_basis(scaled.view(), &[0, 1], &spec)
1275                .expect("rescaled tensor basis")
1276        };
1277        let reference = build_tensor(1.0, 1.0);
1278        for x_scale in [1e-9_f64, 1.0, 1e9] {
1279            for y_scale in [1e-9_f64, 1.0, 1e9] {
1280                let actual = build_tensor(x_scale, y_scale);
1281                assert_build_geometry_close(&actual, &reference, 1e-9);
1282            }
1283        }
1284    }
1285
1286    fn wrapper_basis(family: BasisScaleFamily, abscissa_scale: f64) -> SmoothBasisSpec {
1287        let inner = || Box::new(basis(0, abscissa_scale));
1288        let levels = vec![0.0_f64.to_bits(), 1.0_f64.to_bits(), 2.0_f64.to_bits()];
1289        let factor_smooth = |flavour| SmoothBasisSpec::FactorSmooth {
1290            spec: FactorSmoothSpec {
1291                continuous_cols: vec![0],
1292                group_col: 1,
1293                marginal: open_marginal(abscissa_scale),
1294                flavour,
1295                group_frozen_levels: Some(levels.clone()),
1296                frozen_global_orthogonality: None,
1297            },
1298        };
1299        match family {
1300            BasisScaleFamily::ByVariableNumeric => SmoothBasisSpec::ByVariable {
1301                inner: inner(),
1302                by_col: 1,
1303                kind: BySmoothKind::Numeric,
1304                by: ByVariableSpec::Numeric,
1305            },
1306            BasisScaleFamily::ByVariableFactor => SmoothBasisSpec::ByVariable {
1307                inner: inner(),
1308                by_col: 1,
1309                kind: BySmoothKind::Level {
1310                    level_bits: 1.0_f64.to_bits(),
1311                },
1312                by: ByVariableSpec::Level {
1313                    value_bits: 1.0_f64.to_bits(),
1314                    label: "one".to_string(),
1315                },
1316            },
1317            BasisScaleFamily::FactorSumToZero => SmoothBasisSpec::FactorSumToZero {
1318                inner: inner(),
1319                by_col: 1,
1320                levels,
1321                frozen_global_orthogonality: None,
1322            },
1323            BasisScaleFamily::BySmoothNumeric => SmoothBasisSpec::BySmooth {
1324                smooth: inner(),
1325                by_kind: ByVarKind::Numeric { feature_col: 1 },
1326            },
1327            BasisScaleFamily::BySmoothFactor => SmoothBasisSpec::BySmooth {
1328                smooth: inner(),
1329                by_kind: ByVarKind::Factor {
1330                    feature_col: 1,
1331                    ordered: false,
1332                    frozen_levels: Some(levels),
1333                },
1334            },
1335            BasisScaleFamily::FactorSmoothFs => factor_smooth(FactorSmoothFlavour::Fs {
1336                m_null_penalty_orders: vec![1],
1337            }),
1338            BasisScaleFamily::FactorSmoothSz => factor_smooth(FactorSmoothFlavour::Sz),
1339            BasisScaleFamily::FactorSmoothRe => factor_smooth(FactorSmoothFlavour::Re),
1340            BasisScaleFamily::OpenBSpline
1341            | BasisScaleFamily::CyclicBSpline
1342            | BasisScaleFamily::NaturalCubic
1343            | BasisScaleFamily::ThinPlate
1344            | BasisScaleFamily::SphereWahba
1345            | BasisScaleFamily::SphereHarmonic
1346            | BasisScaleFamily::ConstantCurvature
1347            | BasisScaleFamily::Matern
1348            | BasisScaleFamily::MeasureJet
1349            | BasisScaleFamily::PureDuchon
1350            | BasisScaleFamily::HybridDuchon
1351            | BasisScaleFamily::Pca
1352            | BasisScaleFamily::TensorBSpline => {
1353                panic!("{family:?} is not a wrapper fixture")
1354            }
1355        }
1356    }
1357
1358    fn build_local(data: &Array2<f64>, basis: SmoothBasisSpec) -> LocalSmoothTermBuild {
1359        build_single_local_smooth_term(
1360            data.view(),
1361            &SmoothTermSpec {
1362                frozen_parametric_residualization: None,
1363                name: "scale-contract-wrapper".to_string(),
1364                basis,
1365                shape: ShapeConstraint::None,
1366                joint_null_rotation: None,
1367            },
1368            &mut BasisWorkspace::new(),
1369        )
1370        .expect("wrapper scale fixture must build")
1371    }
1372
1373    #[test]
1374    fn every_wrapper_preserves_its_declared_inner_abscissa_pullback_2315() {
1375        let data = array![
1376            [0.00, 0.0],
1377            [0.08, 1.0],
1378            [0.17, 2.0],
1379            [0.26, 0.0],
1380            [0.35, 1.0],
1381            [0.44, 2.0],
1382            [0.56, 0.0],
1383            [0.65, 1.0],
1384            [0.74, 2.0],
1385            [0.83, 0.0],
1386            [0.92, 1.0],
1387            [1.00, 2.0]
1388        ];
1389        let invariant_families = [
1390            BasisScaleFamily::ByVariableNumeric,
1391            BasisScaleFamily::ByVariableFactor,
1392            BasisScaleFamily::FactorSumToZero,
1393            BasisScaleFamily::BySmoothNumeric,
1394            BasisScaleFamily::BySmoothFactor,
1395            BasisScaleFamily::FactorSmoothFs,
1396            BasisScaleFamily::FactorSmoothSz,
1397        ];
1398        for family in invariant_families {
1399            let reference = build_local(&data, wrapper_basis(family, 1.0));
1400            for factor in [1e-9_f64, 1.0, 1e9] {
1401                let mut scaled = data.clone();
1402                scaled.column_mut(0).mapv_inplace(|value| factor * value);
1403                let actual = build_local(&scaled, wrapper_basis(family, factor));
1404                assert_local_geometry_close(
1405                    &actual,
1406                    &reference,
1407                    2e-8,
1408                    &format!("{family:?} factor={factor:e}"),
1409                );
1410            }
1411        }
1412
1413        // `bs="re"` is a random intercept+slope, so each level carries one
1414        // invariant intercept column and one degree-one slope column.
1415        let family = BasisScaleFamily::FactorSmoothRe;
1416        let reference = build_local(&data, wrapper_basis(family, 1.0));
1417        for factor in [1e-9_f64, 1.0, 1e9] {
1418            let mut scaled = data.clone();
1419            scaled.column_mut(0).mapv_inplace(|value| factor * value);
1420            let actual = build_local(&scaled, wrapper_basis(family, factor));
1421            let mut pulled_back = actual.design.to_dense();
1422            for slope_col in (1..pulled_back.ncols()).step_by(2) {
1423                pulled_back
1424                    .column_mut(slope_col)
1425                    .mapv_inplace(|value| value / factor);
1426            }
1427            assert_matrix_close(&pulled_back, &reference.design.to_dense(), 2e-9);
1428            assert_eq!(
1429                actual.active_penalties.len(),
1430                reference.active_penalties.len()
1431            );
1432            for (observed, target) in actual
1433                .active_penalties
1434                .iter()
1435                .zip(reference.active_penalties.iter())
1436            {
1437                assert_matrix_close(&observed.matrix, &target.matrix, 2e-10);
1438                assert_eq!(observed.nullity, target.nullity);
1439            }
1440            assert_eq!(
1441                joint_unpenalized_dim(actual.dim, &actual.active_penalties),
1442                0,
1443                "the combined random-intercept/slope penalty must be full rank"
1444            );
1445        }
1446    }
1447
1448    #[test]
1449    fn numeric_modulator_has_exact_degree_one_design_and_invariant_penalty_2315() {
1450        let data = array![
1451            [0.00, 0.4],
1452            [0.08, 0.7],
1453            [0.17, 1.1],
1454            [0.26, 0.8],
1455            [0.35, 1.4],
1456            [0.44, 0.6],
1457            [0.56, 1.2],
1458            [0.65, 0.9],
1459            [0.74, 1.5],
1460            [0.83, 0.5],
1461            [0.92, 1.3],
1462            [1.00, 1.0]
1463        ];
1464        for family in [
1465            BasisScaleFamily::ByVariableNumeric,
1466            BasisScaleFamily::BySmoothNumeric,
1467        ] {
1468            let reference = build_local(&data, wrapper_basis(family, 1.0));
1469            for factor in [1e-9_f64, 1.0, 1e9] {
1470                let mut scaled = data.clone();
1471                scaled.column_mut(1).mapv_inplace(|value| factor * value);
1472                let actual = build_local(&scaled, wrapper_basis(family, 1.0));
1473                let pulled_back = actual.design.to_dense().mapv(|value| value / factor);
1474                assert_matrix_close(&pulled_back, &reference.design.to_dense(), 2e-10);
1475                assert_eq!(
1476                    actual.active_penalties.len(),
1477                    reference.active_penalties.len()
1478                );
1479                for (observed, target) in actual
1480                    .active_penalties
1481                    .iter()
1482                    .zip(reference.active_penalties.iter())
1483                {
1484                    assert_matrix_close(&observed.matrix, &target.matrix, 2e-10);
1485                }
1486            }
1487        }
1488    }
1489
1490    fn spherical_spec(method: SphereMethod, radians: bool) -> SphericalSplineBasisSpec {
1491        SphericalSplineBasisSpec {
1492            center_strategy: CenterStrategy::FarthestPoint { num_centers: 6 },
1493            penalty_order: 2,
1494            double_penalty: false,
1495            radians,
1496            method,
1497            max_degree: Some(3),
1498            wahba_kernel: SphereWahbaKernel::Sobolev,
1499            identifiability: SphericalSplineIdentifiability::CenterSumToZero,
1500        }
1501    }
1502
1503    #[test]
1504    fn sphere_constant_curvature_and_pca_obey_their_non_euclidean_gauges_2315() {
1505        let degrees = array![
1506            [-62.0, -150.0],
1507            [-41.0, -77.0],
1508            [-18.0, -12.0],
1509            [4.0, 39.0],
1510            [23.0, 101.0],
1511            [47.0, 166.0],
1512            [66.0, -115.0],
1513            [11.0, -171.0]
1514        ];
1515        let radians = degrees.mapv(f64::to_radians);
1516        for method in [SphereMethod::Wahba, SphereMethod::Harmonic] {
1517            let in_degrees =
1518                build_spherical_spline_basis(degrees.view(), &spherical_spec(method, false))
1519                    .expect("degree-encoded sphere basis");
1520            let in_radians =
1521                build_spherical_spline_basis(radians.view(), &spherical_spec(method, true))
1522                    .expect("radian-encoded sphere basis");
1523            assert_build_geometry_close(&in_radians, &in_degrees, 2e-9);
1524        }
1525
1526        let chart_data = array![
1527            [-0.42, -0.18],
1528            [-0.31, 0.22],
1529            [-0.08, -0.34],
1530            [0.13, 0.29],
1531            [0.27, -0.11],
1532            [0.38, 0.17]
1533        ];
1534        let centers = array![[-0.36, -0.04], [-0.12, 0.25], [0.16, -0.21], [0.34, 0.13]];
1535        let kappa = -0.7_f64;
1536        let length_scale = 0.55_f64;
1537        let reference_spec = ConstantCurvatureBasisSpec {
1538            center_strategy: CenterStrategy::UserProvided(centers.clone()),
1539            kappa,
1540            kappa_fixed: true,
1541            length_scale,
1542            length_scale_fixed: true,
1543            double_penalty: false,
1544            identifiability: ConstantCurvatureIdentifiability::CenterSumToZero,
1545        };
1546        let reference = build_constant_curvature_basis(chart_data.view(), &reference_spec)
1547            .expect("reference constant-curvature basis");
1548        let (_, dk_reference, dkk_reference) = constant_curvature_kernel_kappa_jets(
1549            chart_data.view(),
1550            centers.view(),
1551            kappa,
1552            length_scale,
1553        )
1554        .expect("reference curvature jets");
1555        for factor in [1e-9_f64, 1.0, 1e9] {
1556            let scaled_data = chart_data.mapv(|value| value * factor);
1557            let scaled_centers = centers.mapv(|value| value * factor);
1558            let scaled_kappa = kappa / factor.powi(2);
1559            let scaled_length = length_scale * factor;
1560            let actual = build_constant_curvature_basis(
1561                scaled_data.view(),
1562                &ConstantCurvatureBasisSpec {
1563                    center_strategy: CenterStrategy::UserProvided(scaled_centers.clone()),
1564                    kappa: scaled_kappa,
1565                    kappa_fixed: true,
1566                    length_scale: scaled_length,
1567                    length_scale_fixed: true,
1568                    double_penalty: false,
1569                    identifiability: ConstantCurvatureIdentifiability::CenterSumToZero,
1570                },
1571            )
1572            .expect("rescaled constant-curvature basis");
1573            // Degree ONE, not invariant (gam#2747): the kernel `ℓ·(e^{−d/ℓ}−1)`
1574            // is a length, so the chart similarity multiplies design AND
1575            // penalty by `a`. Dividing it back out is the whole content of the
1576            // declared `ChartSimilarityDegreeOne` / `PhysicalRkhsGramDegreeOne`
1577            // laws — and the fact that ONE factor restores BOTH blocks is the
1578            // check that matters, because a common positive factor on `(X, S)`
1579            // is exactly what `λ` absorbs.
1580            assert_build_geometry_scaled_close(&actual, &reference, factor, 2e-8);
1581            let (_, dk, dkk) = constant_curvature_kernel_kappa_jets(
1582                scaled_data.view(),
1583                scaled_centers.view(),
1584                scaled_kappa,
1585                scaled_length,
1586            )
1587            .expect("rescaled curvature jets");
1588            assert_matrix_close(
1589                &dk.mapv(|value| value / factor.powi(3)),
1590                &dk_reference,
1591                2e-8,
1592            );
1593            assert_matrix_close(
1594                &dkk.mapv(|value| value / factor.powi(5)),
1595                &dkk_reference,
1596                3e-7,
1597            );
1598        }
1599
1600        let pca_data = array![
1601            [-1.0, 0.3],
1602            [-0.4, 1.1],
1603            [0.2, -0.7],
1604            [0.8, 0.5],
1605            [1.3, -0.2],
1606            [1.7, 0.9]
1607        ];
1608        let center_mean = array![0.35, 0.15];
1609        let loadings = array![[0.8, -0.3], [0.6, 0.9]];
1610        let reference = build_pca_smooth_basis(
1611            pca_data.view(),
1612            &[0, 1],
1613            &loadings,
1614            true,
1615            1.7,
1616            Some(&center_mean),
1617            None,
1618            32,
1619        )
1620        .expect("reference PCA score gauge");
1621        for factor in [1e-9_f64, 1.0, 1e9] {
1622            let actual = build_pca_smooth_basis(
1623                pca_data.mapv(|value| factor * value).view(),
1624                &[0, 1],
1625                &loadings.mapv(|value| value / factor),
1626                true,
1627                1.7,
1628                Some(&center_mean.mapv(|value| factor * value)),
1629                None,
1630                32,
1631            )
1632            .expect("rescaled PCA score gauge");
1633            assert_build_geometry_close(&actual, &reference, 2e-10);
1634        }
1635    }
1636
1637    #[test]
1638    fn euclidean_registry_frames_are_similarity_invariant_for_every_declared_family_2315() {
1639        let coordinates = array![
1640            [-0.9, 0.2],
1641            [-0.4, 0.8],
1642            [0.1, -0.6],
1643            [0.7, 0.4],
1644            [1.2, -0.1]
1645        ];
1646        let cases = [
1647            (BasisScaleFamily::ThinPlate, Some(0.7)),
1648            (BasisScaleFamily::Matern, Some(0.7)),
1649            (BasisScaleFamily::MeasureJet, Some(0.7)),
1650            (BasisScaleFamily::PureDuchon, None),
1651            (BasisScaleFamily::HybridDuchon, Some(0.7)),
1652        ];
1653        for (family, length_scale) in cases {
1654            let basis = zoo_basis(family);
1655            let contract = basis.scale_contract();
1656            assert_eq!(contract.family, family);
1657            let mut reference_centers = CenterStrategy::FarthestPoint { num_centers: 3 };
1658            let reference = contract
1659                .normalize_euclidean_frame(
1660                    coordinates.clone(),
1661                    None,
1662                    length_scale,
1663                    &mut reference_centers,
1664                )
1665                .expect("reference scale frame");
1666            for factor in [1e-9_f64, 1.0, 1e9] {
1667                let mut actual_centers = CenterStrategy::FarthestPoint { num_centers: 3 };
1668                let actual = contract
1669                    .normalize_euclidean_frame(
1670                        coordinates.mapv(|value| factor * value),
1671                        None,
1672                        length_scale.map(|ell| factor * ell),
1673                        &mut actual_centers,
1674                    )
1675                    .expect("rescaled frame");
1676                assert_matrix_close(&actual.coordinates, &reference.coordinates, 3e-12);
1677                match (actual.length_scale, reference.length_scale) {
1678                    (Some(observed), Some(target)) => {
1679                        let (observed, target) =
1680                            (observed.standardized_value(), target.standardized_value());
1681                        assert!(
1682                            (observed - target).abs() <= 3e-12 * (1.0 + target.abs()),
1683                            "{family:?} effective range changed at factor {factor}: {observed} vs {target}"
1684                        );
1685                    }
1686                    (None, None) => {}
1687                    pair => panic!("{family:?} changed optional range shape: {pair:?}"),
1688                }
1689            }
1690        }
1691    }
1692
1693    #[test]
1694    fn euclidean_frame_scales_nested_resolved_centers_once_2623() {
1695        let coordinates = array![[-1.2, 0.4], [-0.3, 1.1], [0.6, -0.8], [1.7, 0.2]];
1696        let original_centers = array![[-1.2, 0.4], [0.6, -0.8], [1.7, 0.2]];
1697        let mut strategy = CenterStrategy::DuchonSpectral {
1698            knots: Box::new(CenterStrategy::UserProvided(original_centers.clone())),
1699            basis: crate::basis::DuchonSpectralBasis::Fresh { rank: 3 },
1700        };
1701        let contract = zoo_basis(BasisScaleFamily::PureDuchon).scale_contract();
1702        let fresh = contract
1703            .normalize_euclidean_frame(coordinates.clone(), None, None, &mut strategy)
1704            .expect("fresh Euclidean frame");
1705        let expected_centers = original_centers.mapv(|value| value / fresh.input_scale.get());
1706        let CenterStrategy::DuchonSpectral { knots, .. } = &strategy else {
1707            panic!("Duchon spectral strategy changed shape");
1708        };
1709        let CenterStrategy::UserProvided(fresh_centers) = knots.as_ref() else {
1710            panic!("Duchon spectral knot strategy changed shape");
1711        };
1712        assert_matrix_close(fresh_centers, &expected_centers, 0.0);
1713
1714        let frozen_centers = fresh_centers.clone();
1715        let replay = contract
1716            .normalize_euclidean_frame(coordinates, Some(fresh.input_scale), None, &mut strategy)
1717            .expect("frozen Euclidean replay frame");
1718        assert_matrix_close(&replay.coordinates, &fresh.coordinates, 0.0);
1719        let CenterStrategy::DuchonSpectral { knots, .. } = &strategy else {
1720            panic!("Duchon spectral strategy changed shape on replay");
1721        };
1722        let CenterStrategy::UserProvided(replay_centers) = knots.as_ref() else {
1723            panic!("Duchon spectral knot strategy changed shape on replay");
1724        };
1725        assert_matrix_close(replay_centers, &frozen_centers, 0.0);
1726    }
1727
1728    fn euclidean_basis(family: BasisScaleFamily, factor: f64) -> SmoothBasisSpec {
1729        let centers = CenterStrategy::FarthestPoint { num_centers: 8 };
1730        match family {
1731            BasisScaleFamily::ThinPlate => SmoothBasisSpec::ThinPlate {
1732                feature_cols: vec![0, 1],
1733                spec: ThinPlateBasisSpec {
1734                    center_strategy: centers,
1735                    periodic: None,
1736                    length_scale: 0.55 * factor,
1737                    double_penalty: false,
1738                    identifiability: SpatialIdentifiability::None,
1739                    radial_reparam: None,
1740                },
1741                input_scale: None,
1742            },
1743            BasisScaleFamily::Matern => SmoothBasisSpec::Matern {
1744                feature_cols: vec![0, 1],
1745                spec: MaternBasisSpec {
1746                    center_strategy: centers,
1747                    periodic: None,
1748                    length_scale: crate::basis::MaternLengthScale::fixed(0.55 * factor),
1749                    nu: MaternNu::ThreeHalves,
1750                    include_intercept: false,
1751                    double_penalty: false,
1752                    identifiability: MaternIdentifiability::None,
1753                    aniso_log_scales: None,
1754                },
1755                input_scale: None,
1756            },
1757            BasisScaleFamily::MeasureJet => SmoothBasisSpec::MeasureJet {
1758                feature_cols: vec![0, 1],
1759                spec: MeasureJetBasisSpec {
1760                    center_strategy: centers,
1761                    order_s: 1.5,
1762                    alpha: 1.0,
1763                    tau0: 1e-3,
1764                    num_scales: 3,
1765                    length_scale: 0.55 * factor,
1766                    double_penalty: false,
1767                    learn_length_scale: false,
1768                    multiscale: false,
1769                    identifiability: MeasureJetIdentifiability::CenterSumToZero,
1770                    frozen_quadrature: None,
1771                },
1772                input_scale: None,
1773            },
1774            BasisScaleFamily::PureDuchon | BasisScaleFamily::HybridDuchon => {
1775                SmoothBasisSpec::Duchon {
1776                    feature_cols: vec![0, 1],
1777                    spec: DuchonBasisSpec {
1778                        center_strategy: centers,
1779                        periodic: None,
1780                        length_scale: (family == BasisScaleFamily::HybridDuchon)
1781                            .then_some(0.55 * factor),
1782                        power: if family == BasisScaleFamily::HybridDuchon {
1783                            1.0
1784                        } else {
1785                            0.0
1786                        },
1787                        nullspace_order: DuchonNullspaceOrder::Linear,
1788                        identifiability: SpatialIdentifiability::None,
1789                        aniso_log_scales: None,
1790                        operator_penalties: DuchonOperatorPenaltySpec::all_disabled(),
1791                        boundary: OneDimensionalBoundary::Open,
1792                        radial_reparam: None,
1793                    },
1794                    input_scale: None,
1795                }
1796            }
1797            BasisScaleFamily::ByVariableNumeric
1798            | BasisScaleFamily::ByVariableFactor
1799            | BasisScaleFamily::FactorSumToZero
1800            | BasisScaleFamily::OpenBSpline
1801            | BasisScaleFamily::CyclicBSpline
1802            | BasisScaleFamily::NaturalCubic
1803            | BasisScaleFamily::BySmoothNumeric
1804            | BasisScaleFamily::BySmoothFactor
1805            | BasisScaleFamily::FactorSmoothFs
1806            | BasisScaleFamily::FactorSmoothSz
1807            | BasisScaleFamily::FactorSmoothRe
1808            | BasisScaleFamily::SphereWahba
1809            | BasisScaleFamily::SphereHarmonic
1810            | BasisScaleFamily::ConstantCurvature
1811            | BasisScaleFamily::Pca
1812            | BasisScaleFamily::TensorBSpline => {
1813                panic!("{family:?} is not a Euclidean spatial fixture")
1814            }
1815        }
1816    }
1817
1818    #[test]
1819    fn every_euclidean_spatial_builder_obeys_its_registry_pullback_2315() {
1820        let data = array![
1821            [-0.95, -0.22],
1822            [-0.82, 0.31],
1823            [-0.63, -0.47],
1824            [-0.48, 0.62],
1825            [-0.27, -0.08],
1826            [-0.11, 0.41],
1827            [0.06, -0.55],
1828            [0.21, 0.17],
1829            [0.38, 0.73],
1830            [0.54, -0.31],
1831            [0.69, 0.49],
1832            [0.83, -0.66],
1833            [0.97, 0.04],
1834            [1.08, 0.58],
1835            [1.19, -0.39],
1836            [1.31, 0.26]
1837        ];
1838        for family in [
1839            BasisScaleFamily::ThinPlate,
1840            BasisScaleFamily::Matern,
1841            BasisScaleFamily::MeasureJet,
1842            BasisScaleFamily::PureDuchon,
1843            BasisScaleFamily::HybridDuchon,
1844        ] {
1845            let reference = build_local(&data, euclidean_basis(family, 1.0));
1846            for factor in [1e-9_f64, 1.0, 1e9] {
1847                let scaled = data.mapv(|value| factor * value);
1848                let actual = build_local(&scaled, euclidean_basis(family, factor));
1849                assert_local_geometry_close(
1850                    &actual,
1851                    &reference,
1852                    3e-7,
1853                    &format!("{family:?} factor={factor:e}"),
1854                );
1855            }
1856        }
1857    }
1858
1859    #[test]
1860    fn legacy_or_vector_isotropic_scale_wire_states_are_refused_2319() {
1861        let basis = zoo_basis(BasisScaleFamily::Matern);
1862        let mut legacy = serde_json::to_value(&basis).unwrap();
1863        let fields = legacy
1864            .get_mut("Matern")
1865            .and_then(serde_json::Value::as_object_mut)
1866            .unwrap();
1867        fields.remove("input_scale");
1868        fields.insert("input_scales".to_string(), serde_json::json!([1.0, 2.0]));
1869        assert!(serde_json::from_value::<SmoothBasisSpec>(legacy).is_err());
1870
1871        let mut vector = serde_json::to_value(&basis).unwrap();
1872        vector
1873            .get_mut("Matern")
1874            .and_then(serde_json::Value::as_object_mut)
1875            .unwrap()
1876            .insert("input_scale".to_string(), serde_json::json!([1.0, 1.0]));
1877        assert!(serde_json::from_value::<SmoothBasisSpec>(vector).is_err());
1878    }
1879
1880    #[test]
1881    fn frozen_spatial_spec_requires_the_realized_scalar_frame_2319() {
1882        let frozen = TermCollectionSpec {
1883            linear_terms: Vec::new(),
1884            random_effect_terms: Vec::new(),
1885            smooth_terms: vec![SmoothTermSpec {
1886                frozen_parametric_residualization: None,
1887                name: "matern".to_string(),
1888                basis: zoo_basis(BasisScaleFamily::Matern),
1889                shape: ShapeConstraint::None,
1890                joint_null_rotation: None,
1891            }],
1892        };
1893
1894        let error = frozen
1895            .validate_frozen("model payload")
1896            .expect_err("an unresolved isotropic frame is not frozen");
1897        assert!(error.contains("Matern input_scale is missing"), "{error}");
1898    }
1899}