gam 0.3.116

Generalized penalized likelihood engine
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81

//! Per-axis input standardization and length-scale compensation helpers for
//! the spatial smooth arm.
//!
//! Pure numeric helpers relocated verbatim from `smooth.rs` (issue #780
//! decomposition): per-column variance scales, in-place standardization, the
//! geometric-mean scale, and the kernel length-scale compensation maps that
//! keep the Matérn/Duchon/thin-plate range in original coordinates after
//! standardization. No behavior change — bodies are byte-identical and the
//! parent re-imports each name so every call site is unchanged.

use ndarray::{Array2, ArrayView2};

/// Compute per-column standard deviations for multivariate spatial inputs (d > 1).
/// Returns `None` when d == 1 (standardization unnecessary) or when the caller
/// already supplies frozen scales (prediction path).
pub(crate) fn compute_spatial_input_scales(x: ArrayView2<'_, f64>) -> Option<Vec<f64>> {
    let d = x.ncols();
    if d <= 1 {
        return None;
    }
    let n = x.nrows() as f64;
    if n < 2.0 {
        return None;
    }
    let mut scales = Vec::with_capacity(d);
    for j in 0..d {
        let col = x.column(j);
        let mean = col.sum() / n;
        let var = col.iter().map(|&v| (v - mean).powi(2)).sum::<f64>() / (n - 1.0);
        scales.push(var.sqrt().max(1e-12));
    }
    Some(scales)
}

/// Apply per-column standardization to a data matrix using precomputed scales.
pub(crate) fn apply_input_standardization(x: &mut Array2<f64>, scales: &[f64]) {
    for j in 0..x.ncols() {
        let inv = 1.0 / scales[j];
        x.column_mut(j).mapv_inplace(|v| v * inv);
    }
}

/// Geometric mean of strictly positive scales: `(∏ s_a)^(1/d)`.
///
/// Computed via log-sum-divide to avoid overflow / underflow when d is large
/// or when individual scales are small. The Matérn / Duchon / thin-plate
/// auto-standardization paths use this to compensate the user's
/// `length_scale` so the kernel range remains expressed in *original* data
/// coordinates after per-axis division by σ_a:
///
///   ‖x_std − c_std‖ / L_eff with L_eff = L_user / σ_geom
///
/// matches `‖x − c‖ / L_user` exactly for uniform σ_a (= σ_geom) and reduces
/// to the natural anisotropic-Mahalanobis preconditioning when σ_a vary —
/// the convention σ_geom = (∏σ_a)^(1/d) preserves the kernel volume scale.
fn geometric_mean_scale(scales: &[f64]) -> f64 {
    if scales.is_empty() {
        return 1.0;
    }
    let log_mean: f64 = scales.iter().map(|&s| s.ln()).sum::<f64>() / scales.len() as f64;
    log_mean.exp()
}

pub(crate) fn compensate_length_scale_for_standardization(
    length_scale: f64,
    scales: &[f64],
) -> f64 {
    let sigma_geom = geometric_mean_scale(scales);
    if sigma_geom > 0.0 && sigma_geom.is_finite() {
        length_scale / sigma_geom
    } else {
        length_scale
    }
}

pub(crate) fn compensate_optional_length_scale_for_standardization(
    length_scale: Option<f64>,
    scales: &[f64],
) -> Option<f64> {
    length_scale.map(|l| compensate_length_scale_for_standardization(l, scales))
}