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

gam_terms/smooth/
scale_contract.rs
