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::{BSplineBoundaryConditions, BasisWorkspace, DuchonNullspaceOrder, DuchonOperatorPenaltySpec, MaternIdentifiability, MaternNu, MeasureJetIdentifiability};
778    use ndarray::{Array1, Array2, array};
779    use std::collections::HashSet;
780
781    fn open_marginal(scale: f64) -> BSplineBasisSpec {
782        BSplineBasisSpec {
783            degree: 3,
784            penalty_order: 2,
785            knotspec: BSplineKnotSpec::Provided(Array1::from(vec![
786                0.0,
787                0.0,
788                0.0,
789                0.0,
790                0.25 * scale,
791                0.6 * scale,
792                scale,
793                scale,
794                scale,
795                scale,
796            ])),
797            double_penalty: false,
798            identifiability: BSplineIdentifiability::None,
799            boundary: OneDimensionalBoundary::Open,
800            boundary_conditions: BSplineBoundaryConditions::default(),
801        }
802    }
803
804    fn basis(feature_col: usize, scale: f64) -> SmoothBasisSpec {
805        SmoothBasisSpec::BSpline1D {
806            feature_col,
807            spec: open_marginal(scale),
808        }
809    }
810
811    fn factor_spec(flavour: FactorSmoothFlavour) -> SmoothBasisSpec {
812        SmoothBasisSpec::FactorSmooth {
813            spec: FactorSmoothSpec {
814                continuous_cols: vec![0],
815                group_col: 1,
816                marginal: open_marginal(1.0),
817                flavour,
818                group_frozen_levels: Some(vec![0.0_f64.to_bits(), 1.0_f64.to_bits()]),
819                frozen_global_orthogonality: None,
820            },
821        }
822    }
823
824    fn scale_contract_zoo() -> Vec<SmoothBasisSpec> {
825        let by_level_bits = 1.0_f64.to_bits();
826        let mut sphere_harmonic = SphericalSplineBasisSpec::default();
827        sphere_harmonic.method = SphereMethod::Harmonic;
828        sphere_harmonic.max_degree = Some(3);
829        vec![
830            SmoothBasisSpec::ByVariable {
831                inner: Box::new(basis(0, 1.0)),
832                by_col: 1,
833                kind: BySmoothKind::Numeric,
834                by: ByVariableSpec::Numeric,
835            },
836            SmoothBasisSpec::ByVariable {
837                inner: Box::new(basis(0, 1.0)),
838                by_col: 1,
839                kind: BySmoothKind::Level {
840                    level_bits: by_level_bits,
841                },
842                by: ByVariableSpec::Level {
843                    value_bits: by_level_bits,
844                    label: "one".to_string(),
845                },
846            },
847            SmoothBasisSpec::FactorSumToZero {
848                inner: Box::new(basis(0, 1.0)),
849                by_col: 1,
850                levels: vec![0.0_f64.to_bits(), by_level_bits],
851                frozen_global_orthogonality: None,
852            },
853            basis(0, 1.0),
854            SmoothBasisSpec::BSpline1D {
855                feature_col: 0,
856                spec: BSplineBasisSpec {
857                    knotspec: BSplineKnotSpec::PeriodicUniform {
858                        data_range: (0.0, 1.0),
859                        num_basis: 8,
860                    },
861                    boundary: OneDimensionalBoundary::Cyclic {
862                        start: 0.0,
863                        end: 1.0,
864                    },
865                    ..open_marginal(1.0)
866                },
867            },
868            SmoothBasisSpec::BSpline1D {
869                feature_col: 0,
870                spec: BSplineBasisSpec {
871                    knotspec: BSplineKnotSpec::NaturalCubicRegression {
872                        knots: array![0.0, 0.2, 0.5, 0.8, 1.0],
873                    },
874                    ..open_marginal(1.0)
875                },
876            },
877            SmoothBasisSpec::BySmooth {
878                smooth: Box::new(basis(0, 1.0)),
879                by_kind: ByVarKind::Numeric { feature_col: 1 },
880            },
881            SmoothBasisSpec::BySmooth {
882                smooth: Box::new(basis(0, 1.0)),
883                by_kind: ByVarKind::Factor {
884                    feature_col: 1,
885                    ordered: false,
886                    frozen_levels: Some(vec![0.0_f64.to_bits(), by_level_bits]),
887                },
888            },
889            factor_spec(FactorSmoothFlavour::Fs {}),
890            factor_spec(FactorSmoothFlavour::Sz),
891            factor_spec(FactorSmoothFlavour::Re),
892            SmoothBasisSpec::ThinPlate {
893                feature_cols: vec![0, 1],
894                spec: ThinPlateBasisSpec {
895                    center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
896                    periodic: None,
897                    length_scale: 0.7,
898                    double_penalty: false,
899                    identifiability: SpatialIdentifiability::None,
900                    radial_reparam: None,
901                },
902                input_scale: None,
903            },
904            SmoothBasisSpec::Sphere {
905                feature_cols: vec![0, 1],
906                spec: SphericalSplineBasisSpec::default(),
907            },
908            SmoothBasisSpec::Sphere {
909                feature_cols: vec![0, 1],
910                spec: sphere_harmonic,
911            },
912            SmoothBasisSpec::ConstantCurvature {
913                feature_cols: vec![0, 1],
914                spec: ConstantCurvatureBasisSpec::default(),
915            },
916            SmoothBasisSpec::Matern {
917                feature_cols: vec![0, 1],
918                spec: MaternBasisSpec {
919                    center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
920                    periodic: None,
921                    length_scale: crate::basis::MaternLengthScale::fixed(0.7),
922                    nu: MaternNu::ThreeHalves,
923                    include_intercept: false,
924                    double_penalty: false,
925                    identifiability: MaternIdentifiability::None,
926                    aniso_log_scales: None,
927                },
928                input_scale: None,
929            },
930            SmoothBasisSpec::MeasureJet {
931                feature_cols: vec![0, 1],
932                spec: MeasureJetBasisSpec::default(),
933                input_scale: None,
934            },
935            SmoothBasisSpec::Duchon {
936                feature_cols: vec![0, 1],
937                spec: DuchonBasisSpec {
938                    center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
939                    periodic: None,
940                    length_scale: None,
941                    power: 0.0,
942                    nullspace_order: DuchonNullspaceOrder::Linear,
943                    identifiability: SpatialIdentifiability::None,
944                    aniso_log_scales: None,
945                    operator_penalties: DuchonOperatorPenaltySpec::default(),
946                    boundary: OneDimensionalBoundary::Open,
947                    radial_reparam: None,
948                },
949                input_scale: None,
950            },
951            SmoothBasisSpec::Duchon {
952                feature_cols: vec![0, 1],
953                spec: DuchonBasisSpec {
954                    center_strategy: CenterStrategy::FarthestPoint { num_centers: 8 },
955                    periodic: None,
956                    length_scale: Some(0.7),
957                    power: 1.0,
958                    nullspace_order: DuchonNullspaceOrder::Linear,
959                    identifiability: SpatialIdentifiability::None,
960                    aniso_log_scales: None,
961                    operator_penalties: DuchonOperatorPenaltySpec::default(),
962                    boundary: OneDimensionalBoundary::Open,
963                    radial_reparam: None,
964                },
965                input_scale: None,
966            },
967            SmoothBasisSpec::Pca {
968                feature_cols: vec![0, 1],
969                basis_matrix: array![[1.0], [0.0]],
970                centered: true,
971                smooth_penalty: 1.0,
972                center_mean: Some(array![0.0, 0.0]),
973                pca_basis_path: None,
974                chunk_size: 32,
975            },
976            SmoothBasisSpec::TensorBSpline {
977                feature_cols: vec![0, 1],
978                spec: TensorBSplineSpec {
979                    marginalspecs: vec![open_marginal(1.0), open_marginal(1.0)],
980                    periods: vec![None, None],
981                    double_penalty: false,
982                    identifiability: TensorBSplineIdentifiability::None,
983                    penalty_decomposition: TensorBSplinePenaltyDecomposition::MarginalKroneckerSum,
984                },
985            },
986        ]
987    }
988
989    fn zoo_basis(family: BasisScaleFamily) -> SmoothBasisSpec {
990        scale_contract_zoo()
991            .into_iter()
992            .find(|basis| basis.scale_contract().family == family)
993            .unwrap_or_else(|| panic!("scale-contract zoo is missing {family:?}"))
994    }
995
996    #[test]
997    fn scale_contract_registry_is_exhaustive_unique_and_typed_2315() {
998        let zoo = scale_contract_zoo();
999        let observed: HashSet<_> = zoo
1000            .iter()
1001            .map(|basis| basis.scale_contract().family)
1002            .collect();
1003        let expected: HashSet<_> = BasisScaleFamily::ALL.into_iter().collect();
1004        assert_eq!(zoo.len(), BasisScaleFamily::ALL.len());
1005        assert_eq!(observed, expected);
1006
1007        for basis in &zoo {
1008            let contract = basis.scale_contract();
1009            match contract.family {
1010                BasisScaleFamily::ByVariableNumeric
1011                | BasisScaleFamily::ByVariableFactor
1012                | BasisScaleFamily::FactorSumToZero
1013                | BasisScaleFamily::BySmoothNumeric
1014                | BasisScaleFamily::BySmoothFactor
1015                | BasisScaleFamily::FactorSmoothFs
1016                | BasisScaleFamily::FactorSmoothSz
1017                | BasisScaleFamily::FactorSmoothRe => {
1018                    assert_eq!(contract.children.len(), 1, "{:?}", contract.family);
1019                }
1020                BasisScaleFamily::TensorBSpline => {
1021                    assert_eq!(contract.children.len(), 2);
1022                }
1023                BasisScaleFamily::OpenBSpline
1024                | BasisScaleFamily::CyclicBSpline
1025                | BasisScaleFamily::NaturalCubic
1026                | BasisScaleFamily::ThinPlate
1027                | BasisScaleFamily::SphereWahba
1028                | BasisScaleFamily::SphereHarmonic
1029                | BasisScaleFamily::ConstantCurvature
1030                | BasisScaleFamily::Matern
1031                | BasisScaleFamily::MeasureJet
1032                | BasisScaleFamily::PureDuchon
1033                | BasisScaleFamily::HybridDuchon
1034                | BasisScaleFamily::Pca => assert!(contract.children.is_empty()),
1035            }
1036            let unique_parameters: HashSet<_> = contract
1037                .dimensionful_parameters
1038                .iter()
1039                .map(|parameter| parameter.parameter)
1040                .collect();
1041            assert_eq!(
1042                unique_parameters.len(),
1043                contract.dimensionful_parameters.len(),
1044                "duplicate dimensionful parameter in {:?}",
1045                contract.family
1046            );
1047        }
1048    }
1049
1050    fn assert_matrix_close(actual: &Array2<f64>, expected: &Array2<f64>, tolerance: f64) {
1051        assert_matrix_close_with_context(actual, expected, tolerance, "matrix");
1052    }
1053
1054    fn assert_matrix_close_with_context(
1055        actual: &Array2<f64>,
1056        expected: &Array2<f64>,
1057        tolerance: f64,
1058        context: &str,
1059    ) {
1060        assert_eq!(actual.dim(), expected.dim());
1061        for ((row, col), &target) in expected.indexed_iter() {
1062            let observed = actual[[row, col]];
1063            assert!(
1064                (observed - target).abs() <= tolerance * (1.0 + target.abs()),
1065                "{context}: matrix[{row},{col}] differs: observed={observed:.16e}, target={target:.16e}; actual={actual:?}; expected={expected:?}"
1066            );
1067        }
1068    }
1069
1070    fn assert_build_geometry_close(
1071        actual: &BasisBuildResult,
1072        expected: &BasisBuildResult,
1073        tolerance: f64,
1074    ) {
1075        assert_matrix_close(
1076            &actual.design.to_dense(),
1077            &expected.design.to_dense(),
1078            tolerance,
1079        );
1080        assert_eq!(
1081            actual.active_penalties.len(),
1082            expected.active_penalties.len()
1083        );
1084        for (observed, target) in actual
1085            .active_penalties
1086            .iter()
1087            .zip(expected.active_penalties.iter())
1088        {
1089            assert_eq!(observed.info.source, target.info.source);
1090            assert_eq!(observed.info.effective_rank, target.info.effective_rank);
1091            assert_eq!(observed.nullity, target.nullity);
1092            assert_matrix_close(&observed.matrix, &target.matrix, tolerance);
1093        }
1094    }
1095
1096    fn assert_local_geometry_close(
1097        actual: &LocalSmoothTermBuild,
1098        expected: &LocalSmoothTermBuild,
1099        tolerance: f64,
1100        context: &str,
1101    ) {
1102        assert_matrix_close_with_context(
1103            &actual.design.to_dense(),
1104            &expected.design.to_dense(),
1105            tolerance,
1106            &format!("{context} design"),
1107        );
1108        assert_eq!(
1109            actual.active_penalties.len(),
1110            expected.active_penalties.len()
1111        );
1112        for (observed, target) in actual
1113            .active_penalties
1114            .iter()
1115            .zip(expected.active_penalties.iter())
1116        {
1117            assert_eq!(observed.info.source, target.info.source);
1118            assert_eq!(observed.info.effective_rank, target.info.effective_rank);
1119            assert_eq!(observed.nullity, target.nullity);
1120            assert_matrix_close_with_context(
1121                &observed.matrix,
1122                &target.matrix,
1123                tolerance,
1124                &format!("{context} penalty source {:?}", target.info.source),
1125            );
1126        }
1127    }
1128
1129    fn scaled_cyclic_marginal(scale: f64) -> BSplineBasisSpec {
1130        BSplineBasisSpec {
1131            degree: 3,
1132            penalty_order: 2,
1133            knotspec: BSplineKnotSpec::PeriodicUniform {
1134                data_range: (-0.4 * scale, 1.6 * scale),
1135                num_basis: 9,
1136            },
1137            double_penalty: false,
1138            identifiability: BSplineIdentifiability::None,
1139            boundary: OneDimensionalBoundary::Cyclic {
1140                start: -0.4 * scale,
1141                end: 1.6 * scale,
1142            },
1143            boundary_conditions: BSplineBoundaryConditions::default(),
1144        }
1145    }
1146
1147    fn scaled_natural_marginal(scale: f64) -> BSplineBasisSpec {
1148        BSplineBasisSpec {
1149            degree: 3,
1150            penalty_order: 2,
1151            knotspec: BSplineKnotSpec::NaturalCubicRegression {
1152                knots: array![0.0, 0.17, 0.43, 0.71, 1.0].mapv(|value| value * scale),
1153            },
1154            double_penalty: false,
1155            identifiability: BSplineIdentifiability::None,
1156            boundary: OneDimensionalBoundary::Open,
1157            boundary_conditions: BSplineBoundaryConditions::default(),
1158        }
1159    }
1160
1161    #[test]
1162    fn declared_spline_builders_obey_design_penalty_and_null_scale_laws_2315() {
1163        let open_points = array![0.0, 0.05, 0.3, 0.55, 0.72, 0.95, 1.0];
1164        let cyclic_points = array![-2.4, -0.4, -0.13, 0.2, 1.1, 1.6, 3.6];
1165        let natural_points = array![-0.35, 0.0, 0.09, 0.43, 0.86, 1.0, 1.28];
1166
1167        for (family, points, build_spec) in [
1168            (
1169                BasisScaleFamily::OpenBSpline,
1170                open_points,
1171                open_marginal as fn(f64) -> BSplineBasisSpec,
1172            ),
1173            (
1174                BasisScaleFamily::CyclicBSpline,
1175                cyclic_points,
1176                scaled_cyclic_marginal as fn(f64) -> BSplineBasisSpec,
1177            ),
1178            (
1179                BasisScaleFamily::NaturalCubic,
1180                natural_points,
1181                scaled_natural_marginal as fn(f64) -> BSplineBasisSpec,
1182            ),
1183        ] {
1184            let reference_spec = build_spec(1.0);
1185            assert_eq!(bspline_contract(&reference_spec).family, family);
1186            let reference = build_bspline_basis_1d(points.view(), &reference_spec)
1187                .expect("reference scalar basis");
1188            for factor in [1e-9_f64, 1.0, 1e9] {
1189                let actual = build_bspline_basis_1d(
1190                    points.mapv(|value| value * factor).view(),
1191                    &build_spec(factor),
1192                )
1193                .expect("rescaled scalar basis");
1194                assert_build_geometry_close(&actual, &reference, 8e-10);
1195                for (observed, target) in actual
1196                    .active_penalties
1197                    .iter()
1198                    .zip(reference.active_penalties.iter())
1199                {
1200                    let rescaled = observed.info.normalization_scale * factor.powi(3);
1201                    assert!(
1202                        (rescaled - target.info.normalization_scale).abs()
1203                            <= 2e-8 * (1.0 + target.info.normalization_scale.abs()),
1204                        "{family:?} raw penalty normalizer violated a^-3 at factor {factor}"
1205                    );
1206                }
1207            }
1208        }
1209
1210        let tensor_data = array![
1211            [0.00, 0.13],
1212            [0.08, 0.91],
1213            [0.22, 0.37],
1214            [0.41, 1.00],
1215            [0.58, 0.02],
1216            [0.73, 0.66],
1217            [0.89, 0.48],
1218            [1.00, 0.00]
1219        ];
1220        let build_tensor = |x_scale: f64, y_scale: f64| {
1221            let mut scaled = tensor_data.clone();
1222            scaled.column_mut(0).mapv_inplace(|value| value * x_scale);
1223            scaled.column_mut(1).mapv_inplace(|value| value * y_scale);
1224            let spec = TensorBSplineSpec {
1225                marginalspecs: vec![open_marginal(x_scale), open_marginal(y_scale)],
1226                periods: vec![None, None],
1227                double_penalty: false,
1228                identifiability: TensorBSplineIdentifiability::None,
1229                penalty_decomposition: TensorBSplinePenaltyDecomposition::MarginalKroneckerSum,
1230            };
1231            build_tensor_bspline_basis(scaled.view(), &[0, 1], &spec)
1232                .expect("rescaled tensor basis")
1233        };
1234        let reference = build_tensor(1.0, 1.0);
1235        for x_scale in [1e-9_f64, 1.0, 1e9] {
1236            for y_scale in [1e-9_f64, 1.0, 1e9] {
1237                let actual = build_tensor(x_scale, y_scale);
1238                assert_build_geometry_close(&actual, &reference, 1e-9);
1239            }
1240        }
1241    }
1242
1243    fn wrapper_basis(family: BasisScaleFamily, abscissa_scale: f64) -> SmoothBasisSpec {
1244        let inner = || Box::new(basis(0, abscissa_scale));
1245        let levels = vec![0.0_f64.to_bits(), 1.0_f64.to_bits(), 2.0_f64.to_bits()];
1246        let factor_smooth = |flavour| SmoothBasisSpec::FactorSmooth {
1247            spec: FactorSmoothSpec {
1248                continuous_cols: vec![0],
1249                group_col: 1,
1250                marginal: open_marginal(abscissa_scale),
1251                flavour,
1252                group_frozen_levels: Some(levels.clone()),
1253                frozen_global_orthogonality: None,
1254            },
1255        };
1256        match family {
1257            BasisScaleFamily::ByVariableNumeric => SmoothBasisSpec::ByVariable {
1258                inner: inner(),
1259                by_col: 1,
1260                kind: BySmoothKind::Numeric,
1261                by: ByVariableSpec::Numeric,
1262            },
1263            BasisScaleFamily::ByVariableFactor => SmoothBasisSpec::ByVariable {
1264                inner: inner(),
1265                by_col: 1,
1266                kind: BySmoothKind::Level {
1267                    level_bits: 1.0_f64.to_bits(),
1268                },
1269                by: ByVariableSpec::Level {
1270                    value_bits: 1.0_f64.to_bits(),
1271                    label: "one".to_string(),
1272                },
1273            },
1274            BasisScaleFamily::FactorSumToZero => SmoothBasisSpec::FactorSumToZero {
1275                inner: inner(),
1276                by_col: 1,
1277                levels,
1278                frozen_global_orthogonality: None,
1279            },
1280            BasisScaleFamily::BySmoothNumeric => SmoothBasisSpec::BySmooth {
1281                smooth: inner(),
1282                by_kind: ByVarKind::Numeric { feature_col: 1 },
1283            },
1284            BasisScaleFamily::BySmoothFactor => SmoothBasisSpec::BySmooth {
1285                smooth: inner(),
1286                by_kind: ByVarKind::Factor {
1287                    feature_col: 1,
1288                    ordered: false,
1289                    frozen_levels: Some(levels),
1290                },
1291            },
1292            BasisScaleFamily::FactorSmoothFs => factor_smooth(FactorSmoothFlavour::Fs {}),
1293            BasisScaleFamily::FactorSmoothSz => factor_smooth(FactorSmoothFlavour::Sz),
1294            BasisScaleFamily::FactorSmoothRe => factor_smooth(FactorSmoothFlavour::Re),
1295            BasisScaleFamily::OpenBSpline
1296            | BasisScaleFamily::CyclicBSpline
1297            | BasisScaleFamily::NaturalCubic
1298            | BasisScaleFamily::ThinPlate
1299            | BasisScaleFamily::SphereWahba
1300            | BasisScaleFamily::SphereHarmonic
1301            | BasisScaleFamily::ConstantCurvature
1302            | BasisScaleFamily::Matern
1303            | BasisScaleFamily::MeasureJet
1304            | BasisScaleFamily::PureDuchon
1305            | BasisScaleFamily::HybridDuchon
1306            | BasisScaleFamily::Pca
1307            | BasisScaleFamily::TensorBSpline => {
1308                panic!("{family:?} is not a wrapper fixture")
1309            }
1310        }
1311    }
1312
1313    fn build_local(data: &Array2<f64>, basis: SmoothBasisSpec) -> LocalSmoothTermBuild {
1314        build_single_local_smooth_term(
1315            data.view(),
1316            &SmoothTermSpec {
1317                frozen_parametric_residualization: None,
1318                name: "scale-contract-wrapper".to_string(),
1319                basis,
1320                shape: ShapeConstraint::None,
1321                joint_null_rotation: None,
1322            },
1323            &mut BasisWorkspace::new(),
1324        )
1325        .expect("wrapper scale fixture must build")
1326    }
1327
1328    #[test]
1329    fn every_wrapper_preserves_its_declared_inner_abscissa_pullback_2315() {
1330        let data = array![
1331            [0.00, 0.0],
1332            [0.08, 1.0],
1333            [0.17, 2.0],
1334            [0.26, 0.0],
1335            [0.35, 1.0],
1336            [0.44, 2.0],
1337            [0.56, 0.0],
1338            [0.65, 1.0],
1339            [0.74, 2.0],
1340            [0.83, 0.0],
1341            [0.92, 1.0],
1342            [1.00, 2.0]
1343        ];
1344        let invariant_families = [
1345            BasisScaleFamily::ByVariableNumeric,
1346            BasisScaleFamily::ByVariableFactor,
1347            BasisScaleFamily::FactorSumToZero,
1348            BasisScaleFamily::BySmoothNumeric,
1349            BasisScaleFamily::BySmoothFactor,
1350            BasisScaleFamily::FactorSmoothFs,
1351            BasisScaleFamily::FactorSmoothSz,
1352        ];
1353        for family in invariant_families {
1354            let reference = build_local(&data, wrapper_basis(family, 1.0));
1355            for factor in [1e-9_f64, 1.0, 1e9] {
1356                let mut scaled = data.clone();
1357                scaled.column_mut(0).mapv_inplace(|value| factor * value);
1358                let actual = build_local(&scaled, wrapper_basis(family, factor));
1359                assert_local_geometry_close(
1360                    &actual,
1361                    &reference,
1362                    2e-8,
1363                    &format!("{family:?} factor={factor:e}"),
1364                );
1365            }
1366        }
1367
1368        // `bs="re"` is a random intercept+slope, so each level carries one
1369        // invariant intercept column and one degree-one slope column.
1370        let family = BasisScaleFamily::FactorSmoothRe;
1371        let reference = build_local(&data, wrapper_basis(family, 1.0));
1372        for factor in [1e-9_f64, 1.0, 1e9] {
1373            let mut scaled = data.clone();
1374            scaled.column_mut(0).mapv_inplace(|value| factor * value);
1375            let actual = build_local(&scaled, wrapper_basis(family, factor));
1376            let mut pulled_back = actual.design.to_dense();
1377            for slope_col in (1..pulled_back.ncols()).step_by(2) {
1378                pulled_back
1379                    .column_mut(slope_col)
1380                    .mapv_inplace(|value| value / factor);
1381            }
1382            assert_matrix_close(&pulled_back, &reference.design.to_dense(), 2e-9);
1383            assert_eq!(
1384                actual.active_penalties.len(),
1385                reference.active_penalties.len()
1386            );
1387            for (observed, target) in actual
1388                .active_penalties
1389                .iter()
1390                .zip(reference.active_penalties.iter())
1391            {
1392                assert_matrix_close(&observed.matrix, &target.matrix, 2e-10);
1393                assert_eq!(observed.nullity, target.nullity);
1394            }
1395            assert_eq!(
1396                joint_unpenalized_dim(actual.dim, &actual.active_penalties),
1397                0,
1398                "the combined random-intercept/slope penalty must be full rank"
1399            );
1400        }
1401    }
1402
1403    #[test]
1404    fn numeric_modulator_has_exact_degree_one_design_and_invariant_penalty_2315() {
1405        let data = array![
1406            [0.00, 0.4],
1407            [0.08, 0.7],
1408            [0.17, 1.1],
1409            [0.26, 0.8],
1410            [0.35, 1.4],
1411            [0.44, 0.6],
1412            [0.56, 1.2],
1413            [0.65, 0.9],
1414            [0.74, 1.5],
1415            [0.83, 0.5],
1416            [0.92, 1.3],
1417            [1.00, 1.0]
1418        ];
1419        for family in [
1420            BasisScaleFamily::ByVariableNumeric,
1421            BasisScaleFamily::BySmoothNumeric,
1422        ] {
1423            let reference = build_local(&data, wrapper_basis(family, 1.0));
1424            for factor in [1e-9_f64, 1.0, 1e9] {
1425                let mut scaled = data.clone();
1426                scaled.column_mut(1).mapv_inplace(|value| factor * value);
1427                let actual = build_local(&scaled, wrapper_basis(family, 1.0));
1428                let pulled_back = actual.design.to_dense().mapv(|value| value / factor);
1429                assert_matrix_close(&pulled_back, &reference.design.to_dense(), 2e-10);
1430                assert_eq!(
1431                    actual.active_penalties.len(),
1432                    reference.active_penalties.len()
1433                );
1434                for (observed, target) in actual
1435                    .active_penalties
1436                    .iter()
1437                    .zip(reference.active_penalties.iter())
1438                {
1439                    assert_matrix_close(&observed.matrix, &target.matrix, 2e-10);
1440                }
1441            }
1442        }
1443    }
1444
1445    #[test]
1446    fn euclidean_registry_frames_are_similarity_invariant_for_every_declared_family_2315() {
1447        let coordinates = array![
1448            [-0.9, 0.2],
1449            [-0.4, 0.8],
1450            [0.1, -0.6],
1451            [0.7, 0.4],
1452            [1.2, -0.1]
1453        ];
1454        let cases = [
1455            (BasisScaleFamily::ThinPlate, Some(0.7)),
1456            (BasisScaleFamily::Matern, Some(0.7)),
1457            (BasisScaleFamily::MeasureJet, Some(0.7)),
1458            (BasisScaleFamily::PureDuchon, None),
1459            (BasisScaleFamily::HybridDuchon, Some(0.7)),
1460        ];
1461        for (family, length_scale) in cases {
1462            let basis = zoo_basis(family);
1463            let contract = basis.scale_contract();
1464            assert_eq!(contract.family, family);
1465            let mut reference_centers = CenterStrategy::FarthestPoint { num_centers: 3 };
1466            let reference = contract
1467                .normalize_euclidean_frame(
1468                    coordinates.clone(),
1469                    None,
1470                    length_scale,
1471                    &mut reference_centers,
1472                )
1473                .expect("reference scale frame");
1474            for factor in [1e-9_f64, 1.0, 1e9] {
1475                let mut actual_centers = CenterStrategy::FarthestPoint { num_centers: 3 };
1476                let actual = contract
1477                    .normalize_euclidean_frame(
1478                        coordinates.mapv(|value| factor * value),
1479                        None,
1480                        length_scale.map(|ell| factor * ell),
1481                        &mut actual_centers,
1482                    )
1483                    .expect("rescaled frame");
1484                assert_matrix_close(&actual.coordinates, &reference.coordinates, 3e-12);
1485                match (actual.length_scale, reference.length_scale) {
1486                    (Some(observed), Some(target)) => {
1487                        let (observed, target) =
1488                            (observed.standardized_value(), target.standardized_value());
1489                        assert!(
1490                            (observed - target).abs() <= 3e-12 * (1.0 + target.abs()),
1491                            "{family:?} effective range changed at factor {factor}: {observed} vs {target}"
1492                        );
1493                    }
1494                    (None, None) => {}
1495                    pair => panic!("{family:?} changed optional range shape: {pair:?}"),
1496                }
1497            }
1498        }
1499    }
1500
1501    #[test]
1502    fn euclidean_frame_scales_nested_resolved_centers_once_2623() {
1503        let coordinates = array![[-1.2, 0.4], [-0.3, 1.1], [0.6, -0.8], [1.7, 0.2]];
1504        let original_centers = array![[-1.2, 0.4], [0.6, -0.8], [1.7, 0.2]];
1505        let mut strategy = CenterStrategy::DuchonSpectral {
1506            knots: Box::new(CenterStrategy::UserProvided(original_centers.clone())),
1507            basis: crate::basis::DuchonSpectralBasis::Fresh { rank: 3 },
1508        };
1509        let contract = zoo_basis(BasisScaleFamily::PureDuchon).scale_contract();
1510        let fresh = contract
1511            .normalize_euclidean_frame(coordinates.clone(), None, None, &mut strategy)
1512            .expect("fresh Euclidean frame");
1513        let expected_centers = original_centers.mapv(|value| value / fresh.input_scale.get());
1514        let CenterStrategy::DuchonSpectral { knots, .. } = &strategy else {
1515            panic!("Duchon spectral strategy changed shape");
1516        };
1517        let CenterStrategy::UserProvided(fresh_centers) = knots.as_ref() else {
1518            panic!("Duchon spectral knot strategy changed shape");
1519        };
1520        assert_matrix_close(fresh_centers, &expected_centers, 0.0);
1521
1522        let frozen_centers = fresh_centers.clone();
1523        let replay = contract
1524            .normalize_euclidean_frame(coordinates, Some(fresh.input_scale), None, &mut strategy)
1525            .expect("frozen Euclidean replay frame");
1526        assert_matrix_close(&replay.coordinates, &fresh.coordinates, 0.0);
1527        let CenterStrategy::DuchonSpectral { knots, .. } = &strategy else {
1528            panic!("Duchon spectral strategy changed shape on replay");
1529        };
1530        let CenterStrategy::UserProvided(replay_centers) = knots.as_ref() else {
1531            panic!("Duchon spectral knot strategy changed shape on replay");
1532        };
1533        assert_matrix_close(replay_centers, &frozen_centers, 0.0);
1534    }
1535
1536    fn euclidean_basis(family: BasisScaleFamily, factor: f64) -> SmoothBasisSpec {
1537        let centers = CenterStrategy::FarthestPoint { num_centers: 8 };
1538        match family {
1539            BasisScaleFamily::ThinPlate => SmoothBasisSpec::ThinPlate {
1540                feature_cols: vec![0, 1],
1541                spec: ThinPlateBasisSpec {
1542                    center_strategy: centers,
1543                    periodic: None,
1544                    length_scale: 0.55 * factor,
1545                    double_penalty: false,
1546                    identifiability: SpatialIdentifiability::None,
1547                    radial_reparam: None,
1548                },
1549                input_scale: None,
1550            },
1551            BasisScaleFamily::Matern => SmoothBasisSpec::Matern {
1552                feature_cols: vec![0, 1],
1553                spec: MaternBasisSpec {
1554                    center_strategy: centers,
1555                    periodic: None,
1556                    length_scale: crate::basis::MaternLengthScale::fixed(0.55 * factor),
1557                    nu: MaternNu::ThreeHalves,
1558                    include_intercept: false,
1559                    double_penalty: false,
1560                    identifiability: MaternIdentifiability::None,
1561                    aniso_log_scales: None,
1562                },
1563                input_scale: None,
1564            },
1565            BasisScaleFamily::MeasureJet => SmoothBasisSpec::MeasureJet {
1566                feature_cols: vec![0, 1],
1567                spec: MeasureJetBasisSpec {
1568                    center_strategy: centers,
1569                    order_s: 1.5,
1570                    alpha: 1.0,
1571                    tau0: 1e-3,
1572                    num_scales: 3,
1573                    length_scale: 0.55 * factor,
1574                    double_penalty: false,
1575                    learn_length_scale: false,
1576                    multiscale: false,
1577                    identifiability: MeasureJetIdentifiability::CenterSumToZero,
1578                    frozen_quadrature: None,
1579                },
1580                input_scale: None,
1581            },
1582            BasisScaleFamily::PureDuchon | BasisScaleFamily::HybridDuchon => {
1583                SmoothBasisSpec::Duchon {
1584                    feature_cols: vec![0, 1],
1585                    spec: DuchonBasisSpec {
1586                        center_strategy: centers,
1587                        periodic: None,
1588                        length_scale: (family == BasisScaleFamily::HybridDuchon)
1589                            .then_some(0.55 * factor),
1590                        power: if family == BasisScaleFamily::HybridDuchon {
1591                            1.0
1592                        } else {
1593                            0.0
1594                        },
1595                        nullspace_order: DuchonNullspaceOrder::Linear,
1596                        identifiability: SpatialIdentifiability::None,
1597                        aniso_log_scales: None,
1598                        operator_penalties: DuchonOperatorPenaltySpec::all_disabled(),
1599                        boundary: OneDimensionalBoundary::Open,
1600                        radial_reparam: None,
1601                    },
1602                    input_scale: None,
1603                }
1604            }
1605            BasisScaleFamily::ByVariableNumeric
1606            | BasisScaleFamily::ByVariableFactor
1607            | BasisScaleFamily::FactorSumToZero
1608            | BasisScaleFamily::OpenBSpline
1609            | BasisScaleFamily::CyclicBSpline
1610            | BasisScaleFamily::NaturalCubic
1611            | BasisScaleFamily::BySmoothNumeric
1612            | BasisScaleFamily::BySmoothFactor
1613            | BasisScaleFamily::FactorSmoothFs
1614            | BasisScaleFamily::FactorSmoothSz
1615            | BasisScaleFamily::FactorSmoothRe
1616            | BasisScaleFamily::SphereWahba
1617            | BasisScaleFamily::SphereHarmonic
1618            | BasisScaleFamily::ConstantCurvature
1619            | BasisScaleFamily::Pca
1620            | BasisScaleFamily::TensorBSpline => {
1621                panic!("{family:?} is not a Euclidean spatial fixture")
1622            }
1623        }
1624    }
1625
1626    #[test]
1627    fn every_euclidean_spatial_builder_obeys_its_registry_pullback_2315() {
1628        let data = array![
1629            [-0.95, -0.22],
1630            [-0.82, 0.31],
1631            [-0.63, -0.47],
1632            [-0.48, 0.62],
1633            [-0.27, -0.08],
1634            [-0.11, 0.41],
1635            [0.06, -0.55],
1636            [0.21, 0.17],
1637            [0.38, 0.73],
1638            [0.54, -0.31],
1639            [0.69, 0.49],
1640            [0.83, -0.66],
1641            [0.97, 0.04],
1642            [1.08, 0.58],
1643            [1.19, -0.39],
1644            [1.31, 0.26]
1645        ];
1646        for family in [
1647            BasisScaleFamily::ThinPlate,
1648            BasisScaleFamily::Matern,
1649            BasisScaleFamily::MeasureJet,
1650            BasisScaleFamily::PureDuchon,
1651            BasisScaleFamily::HybridDuchon,
1652        ] {
1653            let reference = build_local(&data, euclidean_basis(family, 1.0));
1654            for factor in [1e-9_f64, 1.0, 1e9] {
1655                let scaled = data.mapv(|value| factor * value);
1656                let actual = build_local(&scaled, euclidean_basis(family, factor));
1657                assert_local_geometry_close(
1658                    &actual,
1659                    &reference,
1660                    3e-7,
1661                    &format!("{family:?} factor={factor:e}"),
1662                );
1663            }
1664        }
1665    }
1666
1667    #[test]
1668    fn legacy_or_vector_isotropic_scale_wire_states_are_refused_2319() {
1669        let basis = zoo_basis(BasisScaleFamily::Matern);
1670        let mut legacy = serde_json::to_value(&basis).unwrap();
1671        let fields = legacy
1672            .get_mut("Matern")
1673            .and_then(serde_json::Value::as_object_mut)
1674            .unwrap();
1675        fields.remove("input_scale");
1676        fields.insert("input_scales".to_string(), serde_json::json!([1.0, 2.0]));
1677        assert!(serde_json::from_value::<SmoothBasisSpec>(legacy).is_err());
1678
1679        let mut vector = serde_json::to_value(&basis).unwrap();
1680        vector
1681            .get_mut("Matern")
1682            .and_then(serde_json::Value::as_object_mut)
1683            .unwrap()
1684            .insert("input_scale".to_string(), serde_json::json!([1.0, 1.0]));
1685        assert!(serde_json::from_value::<SmoothBasisSpec>(vector).is_err());
1686    }
1687
1688    #[test]
1689    fn frozen_spatial_spec_requires_the_realized_scalar_frame_2319() {
1690        let frozen = TermCollectionSpec {
1691            linear_terms: Vec::new(),
1692            random_effect_terms: Vec::new(),
1693            smooth_terms: vec![SmoothTermSpec {
1694                frozen_parametric_residualization: None,
1695                name: "matern".to_string(),
1696                basis: zoo_basis(BasisScaleFamily::Matern),
1697                shape: ShapeConstraint::None,
1698                joint_null_rotation: None,
1699            }],
1700        };
1701
1702        let error = frozen
1703            .validate_frozen("model payload")
1704            .expect_err("an unresolved isotropic frame is not frozen");
1705        assert!(error.contains("Matern input_scale is missing"), "{error}");
1706    }
1707}
gam_terms/smooth/scale_contract.rs

gam_terms/smooth/
scale_contract.rs
