zsfm_tabfm/ensemble/
scalers.rs

1//! Per-column feature scalers, ported from `tabfm/src/classifier_and_regressor.py`'s
2//! `PreprocessingPipeline`, applied in this exact order: `CustomStandardScaler` -> one of 5
3//! optional normalizers (if not `"none"`) -> `OutlierRemover` (last, not second — verified
4//! against `PreprocessingPipeline.fit`). Each scaler exposes `fit`/`transform` mirroring
5//! sklearn's split so the same fitted state (from training columns) can be applied to held-out
6//! columns.
7
8/// `CustomStandardScaler`: `clip((x - mean) / (std + eps), -100, 100)`.
9pub struct CustomStandardScaler {
10    mean: f64,
11    scale: f64,
12}
13
14impl CustomStandardScaler {
15    const EPSILON: f64 = 1e-6;
16    const CLIP: f64 = 100.0;
17
18    pub fn fit(x: &[f64]) -> Self {
19        let mean = mean(x);
20        let scale = std_dev(x, mean, 0) + Self::EPSILON;
21        CustomStandardScaler { mean, scale }
22    }
23
24    pub fn transform(&self, x: &[f64]) -> Vec<f64> {
25        x.iter()
26            .map(|&v| ((v - self.mean) / self.scale).clamp(-Self::CLIP, Self::CLIP))
27            .collect()
28    }
29}
30
31/// `OutlierRemover` (`threshold=4.0` default): two-pass mean/std (outliers masked before the
32/// second pass), then a smooth log-based soft clip (NOT a hard clip) at the recomputed bounds.
33pub struct OutlierRemover {
34    lower_bound: f64,
35    upper_bound: f64,
36}
37
38impl OutlierRemover {
39    pub fn fit(x: &[f64], threshold: f64) -> Self {
40        let mean1 = mean(x);
41        let std1 = std_dev(x, mean1, 1).max(1e-6);
42        let lower1 = mean1 - threshold * std1;
43        let upper1 = mean1 + threshold * std1;
44
45        let clean: Vec<f64> = x
46            .iter()
47            .copied()
48            .filter(|&v| v >= lower1 && v <= upper1)
49            .collect();
50        let (mean2, std2) = if clean.is_empty() {
51            (mean1, std1)
52        } else {
53            let m = mean(&clean);
54            (m, std_dev(&clean, m, 1).max(1e-6))
55        };
56
57        OutlierRemover {
58            lower_bound: mean2 - threshold * std2,
59            upper_bound: mean2 + threshold * std2,
60        }
61    }
62
63    pub fn transform(&self, x: &[f64]) -> Vec<f64> {
64        x.iter()
65            .map(|&v| {
66                let v = (-((v.abs()).ln_1p()) + self.lower_bound).max(v);
67                (v.abs().ln_1p() + self.upper_bound).min(v)
68            })
69            .collect()
70    }
71}
72
73/// Plain `sklearn.preprocessing.StandardScaler` (ddof=0, no epsilon/clip) — used once, globally,
74/// on the raw regression target (`y_scaler_` in `TabFMRegressor`), separate from the per-member
75/// `CustomStandardScaler` applied to features.
76pub struct StandardScaler {
77    mean: f64,
78    scale: f64,
79}
80
81impl StandardScaler {
82    pub fn fit(x: &[f64]) -> Self {
83        let mean = mean(x);
84        let std = std_dev(x, mean, 0);
85        StandardScaler { mean, scale: if std == 0.0 { 1.0 } else { std } }
86    }
87
88    pub fn transform(&self, x: &[f64]) -> Vec<f64> {
89        x.iter().map(|&v| (v - self.mean) / self.scale).collect()
90    }
91
92    pub fn inverse_transform(&self, x: &[f64]) -> Vec<f64> {
93        x.iter().map(|&v| v * self.scale + self.mean).collect()
94    }
95
96    pub fn inverse_transform_scalar(&self, v: f64) -> f64 {
97        v * self.scale + self.mean
98    }
99}
100
101/// `RobustScaler(unit_variance=True)`: `(x - median) / (IQR / 1.349...)`, where the divisor
102/// makes the scale consistent with a standard-normal's std (`1.349... = Φ⁻¹(0.75) - Φ⁻¹(0.25)`).
103pub struct RobustScaler {
104    median: f64,
105    scale: f64,
106}
107
108impl RobustScaler {
109    pub fn fit(x: &[f64]) -> Self {
110        let mut sorted = x.to_vec();
111        sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
112        let median = percentile_linear(&sorted, 50.0);
113        let q25 = percentile_linear(&sorted, 25.0);
114        let q75 = percentile_linear(&sorted, 75.0);
115        let iqr = q75 - q25;
116        let norm_factor = norm_ppf(0.75) - norm_ppf(0.25); // ~1.3489795...
117        let scale = if iqr == 0.0 { 1.0 } else { iqr / norm_factor };
118        RobustScaler { median, scale }
119    }
120
121    pub fn transform(&self, x: &[f64]) -> Vec<f64> {
122        x.iter().map(|&v| (v - self.median) / self.scale).collect()
123    }
124}
125
126/// `QuantileTransformer(output_distribution="normal")`: empirical CDF (via linear-interpolated
127/// percentiles, sklearn's default `n_quantiles=1000` capped at the sample count) mapped through
128/// the inverse standard-normal CDF.
129pub struct QuantileTransformer {
130    references: Vec<f64>, // uniform grid in [0,1], length n_quantiles
131    quantiles: Vec<f64>,  // data values at each reference quantile, monotonic non-decreasing
132}
133
134impl QuantileTransformer {
135    const BOUNDS_THRESHOLD: f64 = 1e-7;
136
137    pub fn fit(x: &[f64]) -> Self {
138        let n_quantiles = 1000usize.min(x.len().max(1));
139        let mut sorted = x.to_vec();
140        sorted.sort_by(|a, b| a.partial_cmp(b).unwrap());
141
142        let references: Vec<f64> = (0..n_quantiles)
143            .map(|i| i as f64 / (n_quantiles - 1).max(1) as f64)
144            .collect();
145        let mut quantiles: Vec<f64> =
146            references.iter().map(|&r| percentile_linear(&sorted, r * 100.0)).collect();
147        // Enforce monotonicity (sklearn does this to guard against numerical noise).
148        for i in 1..quantiles.len() {
149            if quantiles[i] < quantiles[i - 1] {
150                quantiles[i] = quantiles[i - 1];
151            }
152        }
153        QuantileTransformer { references, quantiles }
154    }
155
156    pub fn transform(&self, x: &[f64]) -> Vec<f64> {
157        x.iter()
158            .map(|&v| {
159                let u = interp_monotonic(v, &self.quantiles, &self.references)
160                    .clamp(Self::BOUNDS_THRESHOLD, 1.0 - Self::BOUNDS_THRESHOLD);
161                norm_ppf(u)
162            })
163            .collect()
164    }
165}
166
167/// `PowerTransformer(method="yeo-johnson", standardize=True)`: per-feature MLE-fit `lambda`
168/// (via a ported `scipy.optimize.brent`), then standardize the transformed values.
169pub struct PowerTransformer {
170    lambda: f64,
171    post_mean: f64,
172    post_std: f64,
173}
174
175impl PowerTransformer {
176    pub fn fit(x: &[f64]) -> Self {
177        let lambda = yeo_johnson_optimize(x);
178        let transformed: Vec<f64> = x.iter().map(|&v| yeo_johnson_transform(v, lambda)).collect();
179        let post_mean = mean(&transformed);
180        let post_std = std_dev(&transformed, post_mean, 0).max(1e-300);
181        PowerTransformer { lambda, post_mean, post_std }
182    }
183
184    pub fn transform(&self, x: &[f64]) -> Vec<f64> {
185        x.iter()
186            .map(|&v| (yeo_johnson_transform(v, self.lambda) - self.post_mean) / self.post_std)
187            .collect()
188    }
189}
190
191/// One member's full `PreprocessingPipeline` for a single column: fit on `train_col`, apply to
192/// both `train_col` and `test_col`. Order: `CustomStandardScaler` -> normalizer (if not
193/// `NormMethod::None`) -> `OutlierRemover`.
194pub fn apply_pipeline(
195    train_col: &[f64],
196    test_col: &[f64],
197    norm_method: super::config_gen::NormMethod,
198    outlier_threshold: f64,
199) -> (Vec<f64>, Vec<f64>) {
200    use super::config_gen::NormMethod;
201
202    let scaler = CustomStandardScaler::fit(train_col);
203    let mut train = scaler.transform(train_col);
204    let mut test = scaler.transform(test_col);
205
206    match norm_method {
207        NormMethod::None => {}
208        NormMethod::Power => {
209            let n = PowerTransformer::fit(&train);
210            train = n.transform(&train);
211            test = n.transform(&test);
212        }
213        NormMethod::Quantile => {
214            let n = QuantileTransformer::fit(&train);
215            train = n.transform(&train);
216            test = n.transform(&test);
217        }
218        NormMethod::QuantileRtdl => {
219            // Noise injection uses our own seeded RNG (not bit-identical to NumPy's
220            // default_rng/PCG64) — a documented gap, see the plan's scoping note; the
221            // downstream QuantileTransformer + StandardScaler math is otherwise exact.
222            let noisy = inject_rtdl_noise(&train);
223            let n = QuantileTransformer::fit(&noisy);
224            let train_q = n.transform(&noisy);
225            let test_q = n.transform(&test);
226            let std = StandardScaler::fit(&train_q);
227            train = std.transform(&train_q);
228            test = std.transform(&test_q);
229        }
230        NormMethod::Robust => {
231            let n = RobustScaler::fit(&train);
232            train = n.transform(&train);
233            test = n.transform(&test);
234        }
235    }
236
237    let outlier = OutlierRemover::fit(&train, outlier_threshold);
238    (outlier.transform(&train), outlier.transform(&test))
239}
240
241/// `RTDLQuantileTransformer`'s noise-injection step (`noise=1e-3` default): `x + (noise /
242/// max(std(x), noise)) * standard_normal(shape)`. Uses a simple seeded LCG-derived Gaussian
243/// (Box-Muller), not NumPy's PCG64 — see `apply_pipeline`'s doc comment.
244fn inject_rtdl_noise(x: &[f64]) -> Vec<f64> {
245    let std = std_dev(x, mean(x), 0);
246    let noise = 1e-3;
247    let noise_std = noise / std.max(noise);
248    let mut state = 0x2545F4914F6CDD1Du64;
249    let mut next_u64 = move || {
250        state ^= state << 13;
251        state ^= state >> 7;
252        state ^= state << 17;
253        state
254    };
255    x.iter()
256        .map(|&v| {
257            let u1 = (next_u64() >> 11) as f64 / (1u64 << 53) as f64;
258            let u2 = (next_u64() >> 11) as f64 / (1u64 << 53) as f64;
259            let z = (-2.0 * u1.max(1e-300).ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos();
260            v + noise_std * z
261        })
262        .collect()
263}
264
265fn yeo_johnson_transform(x: f64, lambda: f64) -> f64 {
266    const EPS: f64 = 1e-6;
267    if x >= 0.0 {
268        if lambda.abs() > EPS {
269            ((x + 1.0).powf(lambda) - 1.0) / lambda
270        } else {
271            (x + 1.0).ln()
272        }
273    } else if (lambda - 2.0).abs() > EPS {
274        -((-x + 1.0).powf(2.0 - lambda) - 1.0) / (2.0 - lambda)
275    } else {
276        -(-x + 1.0).ln()
277    }
278}
279
280/// Negative log-likelihood of the Yeo-Johnson transform at `lambda`, matching sklearn's
281/// `PowerTransformer._yeo_johnson_optimize`'s objective (maximized there; minimized here).
282fn yeo_johnson_neg_log_likelihood(x: &[f64], lambda: f64) -> f64 {
283    let n = x.len() as f64;
284    let transformed: Vec<f64> = x.iter().map(|&v| yeo_johnson_transform(v, lambda)).collect();
285    let m = mean(&transformed);
286    let var = variance(&transformed, m, 0);
287    if var < f64::MIN_POSITIVE {
288        return f64::INFINITY;
289    }
290    let mut loglike = -n / 2.0 * var.ln();
291    let sum_term: f64 = x.iter().map(|&v| v.signum() * v.abs().ln_1p()).sum();
292    loglike += (lambda - 1.0) * sum_term;
293    -loglike
294}
295
296fn yeo_johnson_optimize(x: &[f64]) -> f64 {
297    brent_minimize(|lambda| yeo_johnson_neg_log_likelihood(x, lambda), -2.0, 2.0)
298}
299
300/// A 1-D bracket-then-Brent minimizer, matching the shape of `scipy.optimize.brent`: first grows
301/// an `(a, b, c)` triplet bracketing a minimum from the two starting points (golden-ratio
302/// expansion), then refines with Brent's method (golden section + parabolic interpolation).
303fn brent_minimize(f: impl Fn(f64) -> f64, xa0: f64, xb0: f64) -> f64 {
304    const GOLD: f64 = 1.618_034;
305    const GLIMIT: f64 = 100.0;
306    const TINY: f64 = 1e-21;
307
308    let (mut ax, mut bx) = (xa0, xb0);
309    let (mut fa, mut fb) = (f(ax), f(bx));
310    if fa < fb {
311        std::mem::swap(&mut ax, &mut bx);
312        std::mem::swap(&mut fa, &mut fb);
313    }
314    let mut cx = bx + GOLD * (bx - ax);
315    let mut fc = f(cx);
316
317    while fc < fb {
318        let r = (bx - ax) * (fb - fc);
319        let q = (bx - cx) * (fb - fa);
320        let denom = 2.0 * (q - r).abs().max(TINY) * (q - r).signum();
321        let mut u = bx - ((bx - cx) * q - (bx - ax) * r) / denom;
322        let ulim = bx + GLIMIT * (cx - bx);
323
324        let fu;
325        if (bx - u) * (u - cx) > 0.0 {
326            fu = f(u);
327            if fu < fc {
328                ax = bx;
329                bx = u;
330                fa = fb;
331                fb = fu;
332                break;
333            } else if fu > fb {
334                cx = u;
335                fc = fu;
336                break;
337            }
338            u = cx + GOLD * (cx - bx);
339            let fu2 = f(u);
340            ax = bx; bx = cx; cx = u;
341            fa = fb; fb = fc; fc = fu2;
342        } else if (cx - u) * (u - ulim) > 0.0 {
343            fu = f(u);
344            if fu < fc {
345                bx = cx; cx = u; let u2 = cx + GOLD * (cx - bx);
346                fb = fc; fc = fu; let fu2 = f(u2);
347                ax = bx; bx = cx; cx = u2; fa = fb; fb = fc; fc = fu2;
348            } else {
349                ax = bx; bx = cx; cx = u;
350                fa = fb; fb = fc; fc = fu;
351            }
352        } else if (u - ulim) * (ulim - cx) >= 0.0 {
353            u = ulim;
354            fu = f(u);
355            ax = bx; bx = cx; cx = u;
356            fa = fb; fb = fc; fc = fu;
357        } else {
358            u = cx + GOLD * (cx - bx);
359            fu = f(u);
360            ax = bx; bx = cx; cx = u;
361            fa = fb; fb = fc; fc = fu;
362        }
363        if (cx - ax).abs() > 1e6 {
364            break; // safety valve against runaway brackets on pathological inputs
365        }
366    }
367
368    // Ensure ax < cx for Brent's method proper.
369    let (mut a, mut c) = if ax < cx { (ax, cx) } else { (cx, ax) };
370    let mut x = bx;
371    let mut w = bx;
372    let mut v = bx;
373    let mut fx = f(x);
374    let mut fw = fx;
375    let mut fv = fx;
376    let mut d = 0.0f64;
377    let mut e = 0.0f64;
378    const CGOLD: f64 = 0.381_966;
379    const ZEPS: f64 = 1e-12;
380    const TOL: f64 = 1e-8;
381
382    for _ in 0..100 {
383        let xm = 0.5 * (a + c);
384        let tol1 = TOL * x.abs() + ZEPS;
385        let tol2 = 2.0 * tol1;
386        if (x - xm).abs() <= tol2 - 0.5 * (c - a) {
387            break;
388        }
389        let mut use_golden = true;
390        if e.abs() > tol1 {
391            let r = (x - w) * (fx - fv);
392            let mut q = (x - v) * (fx - fw);
393            let mut p = (x - v) * q - (x - w) * r;
394            q = 2.0 * (q - r);
395            if q > 0.0 {
396                p = -p;
397            }
398            q = q.abs();
399            let etemp = e;
400            e = d;
401            if p.abs() < (0.5 * q * etemp).abs() && p > q * (a - x) && p < q * (c - x) {
402                d = p / q;
403                let u = x + d;
404                if u - a < tol2 || c - u < tol2 {
405                    d = if xm - x >= 0.0 { tol1 } else { -tol1 };
406                }
407                use_golden = false;
408            }
409        }
410        if use_golden {
411            e = if x >= xm { a - x } else { c - x };
412            d = CGOLD * e;
413        }
414        let u = if d.abs() >= tol1 { x + d } else { x + if d >= 0.0 { tol1 } else { -tol1 } };
415        let fu = f(u);
416        if fu <= fx {
417            if u >= x {
418                a = x;
419            } else {
420                c = x;
421            }
422            v = w; fv = fw;
423            w = x; fw = fx;
424            x = u; fx = fu;
425        } else {
426            if u < x {
427                a = u;
428            } else {
429                c = u;
430            }
431            if fu <= fw || w == x {
432                v = w; fv = fw;
433                w = u; fw = fu;
434            } else if fu <= fv || v == x || v == w {
435                v = u; fv = fu;
436            }
437        }
438    }
439    x
440}
441
442// ---------------------------------------------------------------------------
443// Shared numeric helpers
444// ---------------------------------------------------------------------------
445
446fn mean(x: &[f64]) -> f64 {
447    if x.is_empty() { 0.0 } else { x.iter().sum::<f64>() / x.len() as f64 }
448}
449
450fn variance(x: &[f64], mean: f64, ddof: usize) -> f64 {
451    let n = x.len();
452    if n <= ddof {
453        return 0.0;
454    }
455    let sum_sq: f64 = x.iter().map(|&v| (v - mean).powi(2)).sum();
456    sum_sq / (n - ddof) as f64
457}
458
459fn std_dev(x: &[f64], mean: f64, ddof: usize) -> f64 {
460    variance(x, mean, ddof).sqrt()
461}
462
463/// NumPy's default ("linear") percentile interpolation on an already-sorted slice.
464pub fn percentile_linear(sorted: &[f64], p: f64) -> f64 {
465    let n = sorted.len();
466    if n == 0 {
467        return f64::NAN;
468    }
469    if n == 1 {
470        return sorted[0];
471    }
472    let rank = (p / 100.0) * (n - 1) as f64;
473    let lo = rank.floor() as usize;
474    let hi = rank.ceil() as usize;
475    if lo == hi {
476        sorted[lo]
477    } else {
478        let frac = rank - lo as f64;
479        sorted[lo] * (1.0 - frac) + sorted[hi] * frac
480    }
481}
482
483/// `np.interp(x, xp, fp)`-equivalent: linear interpolation with clamping outside the domain,
484/// for a monotonic non-decreasing `xp`.
485fn interp_monotonic(x: f64, xp: &[f64], fp: &[f64]) -> f64 {
486    let n = xp.len();
487    if x <= xp[0] {
488        return fp[0];
489    }
490    if x >= xp[n - 1] {
491        return fp[n - 1];
492    }
493    // Binary search for the interval containing x.
494    let idx = xp.partition_point(|&v| v <= x);
495    let (x0, x1) = (xp[idx - 1], xp[idx]);
496    let (y0, y1) = (fp[idx - 1], fp[idx]);
497    if x1 == x0 {
498        y0
499    } else {
500        y0 + (y1 - y0) * (x - x0) / (x1 - x0)
501    }
502}
503
504/// Inverse standard-normal CDF (`scipy.stats.norm.ppf`), via Acklam's rational approximation
505/// with one Halley's-method refinement step (accurate to ~1e-9).
506pub fn norm_ppf(p: f64) -> f64 {
507    if p <= 0.0 {
508        return f64::NEG_INFINITY;
509    }
510    if p >= 1.0 {
511        return f64::INFINITY;
512    }
513    const A: [f64; 6] = [
514        -3.969_683_028_665_376e+01, 2.209_460_984_245_205e+02, -2.759_285_104_469_687e+02,
515        1.383_577_518_672_690e+02, -3.066_479_806_614_716e+01, 2.506_628_277_459_239e+00,
516    ];
517    const B: [f64; 5] = [
518        -5.447_609_879_822_406e+01, 1.615_858_368_580_409e+02, -1.556_989_798_598_866e+02,
519        6.680_131_188_771_972e+01, -1.328_068_155_288_572e+01,
520    ];
521    const C: [f64; 6] = [
522        -7.784_894_002_430_293e-03, -3.223_964_580_411_365e-01, -2.400_758_277_161_838e+00,
523        -2.549_732_539_343_734e+00, 4.374_664_141_464_968e+00, 2.938_163_982_698_783e+00,
524    ];
525    const D: [f64; 4] = [
526        7.784_695_709_041_462e-03, 3.224_671_290_700_398e-01, 2.445_134_137_142_996e+00,
527        3.754_408_661_907_416e+00,
528    ];
529    const P_LOW: f64 = 0.02425;
530    let p_high = 1.0 - P_LOW;
531
532    let mut x = if p < P_LOW {
533        let q = (-2.0 * p.ln()).sqrt();
534        (((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
535            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
536    } else if p <= p_high {
537        let q = p - 0.5;
538        let r = q * q;
539        (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5]) * q
540            / (((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0)
541    } else {
542        let q = (-2.0 * (1.0 - p).ln()).sqrt();
543        -(((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
544            / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
545    };
546
547    // One Halley's-method refinement using the standard normal CDF/PDF.
548    let e = 0.5 * erfc(-x / std::f64::consts::SQRT_2) - p;
549    let u = e * (2.0 * std::f64::consts::PI).sqrt() * (x * x / 2.0).exp();
550    x -= u / (1.0 + x * u / 2.0);
551    x
552}
553
554fn erfc(x: f64) -> f64 {
555    1.0 - erf(x)
556}
557
558/// Abramowitz-Stegun 7.1.26 rational approximation to `erf` (max error ~1.5e-7) — sufficient
559/// precision for one Halley refinement step above.
560fn erf(x: f64) -> f64 {
561    let sign = if x < 0.0 { -1.0 } else { 1.0 };
562    let x = x.abs();
563    const A1: f64 = 0.254_829_592;
564    const A2: f64 = -0.284_496_736;
565    const A3: f64 = 1.421_413_741;
566    const A4: f64 = -1.453_152_027;
567    const A5: f64 = 1.061_405_429;
568    const P: f64 = 0.327_591_1;
569    let t = 1.0 / (1.0 + P * x);
570    let y = 1.0 - (((((A5 * t + A4) * t) + A3) * t + A2) * t + A1) * t * (-x * x).exp();
571    sign * y
572}
573
574#[cfg(test)]
575mod tests {
576    use super::*;
577
578    #[test]
579    fn test_custom_standard_scaler() {
580        let x = vec![1.0, 2.0, 3.0, 4.0, 5.0];
581        let s = CustomStandardScaler::fit(&x);
582        let t = s.transform(&x);
583        // mean=3, std (ddof=0) = sqrt(2) ~ 1.41421356
584        assert!((t[0] - (1.0 - 3.0) / (2f64.sqrt() + 1e-6)).abs() < 1e-9);
585        assert!((t[2]).abs() < 1e-9); // middle value maps to ~0
586    }
587
588    #[test]
589    fn test_robust_scaler_matches_scipy_norm_factor() {
590        let x: Vec<f64> = (0..101).map(|i| i as f64).collect();
591        let s = RobustScaler::fit(&x);
592        // median of 0..100 is 50; q25=25, q75=75, iqr=50
593        assert!((s.median - 50.0).abs() < 1e-9);
594    }
595
596    #[test]
597    fn test_norm_ppf_known_values() {
598        assert!((norm_ppf(0.5)).abs() < 1e-6);
599        assert!((norm_ppf(0.975) - 1.959_963_985).abs() < 1e-4);
600        assert!((norm_ppf(0.025) + 1.959_963_985).abs() < 1e-4);
601    }
602
603    #[test]
604    fn test_yeo_johnson_identity_at_lambda_1() {
605        assert!((yeo_johnson_transform(5.0, 1.0) - 5.0).abs() < 1e-9);
606        assert!((yeo_johnson_transform(-3.0, 1.0) - (-3.0)).abs() < 1e-9);
607    }
608
609    #[test]
610    fn test_power_transformer_matches_sklearn() {
611        // Generated via: sklearn.preprocessing.PowerTransformer(method="yeo-johnson",
612        // standardize=True) fit on np.random.default_rng(0).uniform(-2, 2, 20).
613        let x = vec![
614            0.5478467492858172, -0.9208531449445188, -1.8361059042552212, -1.9338894578858836,
615            1.2530809568010897, 1.6510223091108869, 0.42654310306871945, 0.9179862439359936,
616            0.17449996586169148, 1.740289695151073, 1.2634142164861286, -1.9890459993194076,
617            1.4296171063502774, -1.8656576987781426, 0.9186217857197763, -1.297377517589764,
618            1.4527156893995463, 0.16584488099636685, -0.8011524378504609, -0.3092511152093662,
619        ];
620        let expected_lambda = 1.345222766603304;
621        let expected_transformed = [
622            0.2719479783598148, -0.8318522439585222, -1.3650247006836644, -1.4181615227275621,
623            0.9606112117567647, 1.3853857310641091, 0.16297351213481462, 0.6223628894658091,
624            -0.05314385725391679, 1.483863071047962, 0.9713347980371447, -1.4478643497510573,
625            1.1461058305187737, -1.3811492793940139, 0.6229863176690302, -1.0599639970952783,
626            1.1707302052556954, -0.06030208166470946, -0.7561721820272899, -0.42466733075390517,
627        ];
628        let pt = PowerTransformer::fit(&x);
629        assert!(
630            (pt.lambda - expected_lambda).abs() < 1e-4,
631            "lambda {} vs expected {}",
632            pt.lambda,
633            expected_lambda
634        );
635        let got = pt.transform(&x);
636        for (g, e) in got.iter().zip(expected_transformed.iter()) {
637            assert!((g - e).abs() < 1e-3, "got {g} vs expected {e}");
638        }
639    }
640
641    #[test]
642    fn test_percentile_linear() {
643        let sorted = vec![1.0, 2.0, 3.0, 4.0, 5.0];
644        assert!((percentile_linear(&sorted, 50.0) - 3.0).abs() < 1e-9);
645        assert!((percentile_linear(&sorted, 0.0) - 1.0).abs() < 1e-9);
646        assert!((percentile_linear(&sorted, 100.0) - 5.0).abs() < 1e-9);
647    }
648
649    #[test]
650    fn test_quantile_transformer_matches_sklearn() {
651        // Generated via: sklearn.preprocessing.QuantileTransformer(output_distribution="normal",
652        // random_state=0) fit on np.random.default_rng(0).uniform(-2, 2, 20). n_quantiles is
653        // capped to n_samples=20 by sklearn since 1000 > 20.
654        let x = vec![
655            0.5478467492858172, -0.9208531449445188, -1.8361059042552212, -1.9338894578858836,
656            1.2530809568010897, 1.6510223091108869, 0.42654310306871945, 0.9179862439359936,
657            0.17449996586169148, 1.740289695151073, 1.2634142164861286, -1.9890459993194076,
658            1.4296171063502774, -1.8656576987781426, 0.9186217857197763, -1.297377517589764,
659            1.4527156893995463, 0.16584488099636685, -0.8011524378504609, -0.3092511152093662,
660        ];
661        let expected = [
662            0.199201324789267, -0.6336400007797011, -1.003147967662534, -1.6198562586382699,
663            0.633640000779701, 1.6198562586382697, 0.0660118123758406, 0.33603814037182306,
664            -0.06601181237584074, 5.19933758270342, 0.8045963803603002, -5.199337582605575,
665            1.0031479676625337, -1.2521195202652193, 0.47950565333094985, -0.8045963803603002,
666            1.2521195202652189, -0.199201324789267, -0.47950565333095013, -0.33603814037182317,
667        ];
668        let qt = QuantileTransformer::fit(&x);
669        let got = qt.transform(&x);
670        for (g, e) in got.iter().zip(expected.iter()) {
671            assert!((g - e).abs() < 1e-3, "got {g} vs expected {e}");
672        }
673    }
674}
zsfm_tabfm/ensemble/scalers.rs

zsfm_tabfm/ensemble/
scalers.rs
