baselines 0.1.0

Baseline correction algorithms for signals, spectra, and images
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
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

//! Asymmetrically reweighted penalized least-squares smoothing.

use crate::Result;
use crate::fit::{Fit, FitReport};
use crate::whittaker::engine::{
    Reweighter, WhittakerParams, WhittakerWorkspace, fit_alloc, fit_into, relative_change,
};
use crate::workspace::logistic;

/// Parameters for [`arpls`].
#[derive(Debug, Clone, Copy, Default, PartialEq)]
pub struct ArPlsParams {
    /// Shared Whittaker parameters.
    pub whittaker: WhittakerParams,
}

/// Fits an arPLS baseline.
///
/// # References
///
/// - J. Baek et al., "Baseline correction using asymmetrically reweighted
///   penalized least squares smoothing", *Analyst*, 2015.
/// - `pybaselines.Baseline.arpls` is used as a behavioral reference.
pub fn arpls(y: &[f64], params: ArPlsParams) -> Result<Fit> {
    params.whittaker.validate()?;
    fit_alloc(y, params.whittaker, ArPlsWeights)
}

/// Fits an arPLS baseline into an existing output buffer.
pub fn arpls_into(
    y: &[f64],
    params: ArPlsParams,
    baseline: &mut [f64],
    workspace: &mut WhittakerWorkspace,
) -> Result<FitReport> {
    params.whittaker.validate()?;
    fit_into(y, params.whittaker, ArPlsWeights, baseline, workspace)
}

#[derive(Debug, Clone, Copy)]
struct ArPlsWeights;

impl Reweighter for ArPlsWeights {
    fn initialize(&self, _y: &[f64], weights: &mut [f64]) {
        weights.fill(1.0);
    }

    fn update(&self, y: &[f64], baseline: &[f64], weights: &mut [f64], _iter: usize) -> f64 {
        let previous = weights.to_vec();
        let negative: Vec<f64> = y
            .iter()
            .zip(baseline)
            .map(|(observed, fitted)| observed - fitted)
            .filter(|residual| *residual < 0.0)
            .collect();

        if negative.is_empty() {
            weights.fill(1.0);
            return relative_change(&previous, weights);
        }

        let mean = negative.iter().sum::<f64>() / negative.len() as f64;
        let variance = negative
            .iter()
            .map(|value| {
                let centered = value - mean;
                centered * centered
            })
            .sum::<f64>()
            / negative.len() as f64;
        let std = variance.sqrt().max(f64::EPSILON);

        for ((weight, observed), fitted) in weights.iter_mut().zip(y).zip(baseline) {
            let residual = observed - fitted;
            let exponent = 2.0 * (residual - (2.0 * std - mean)) / std;
            *weight = 1.0 - logistic(exponent);
        }

        relative_change(&previous, weights)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::whittaker::asls::{AslsParams, asls};

    #[test]
    fn whittaker_methods_return_finite_baselines() {
        let y: Vec<f64> = (0..100)
            .map(|i| {
                let x = i as f64 / 99.0;
                0.5 + 0.2 * x + (-(x - 0.45).powi(2) / 0.002).exp()
            })
            .collect();

        let asls_fit = asls(&y, AslsParams::default()).unwrap();
        let arpls_fit = arpls(&y, ArPlsParams::default()).unwrap();

        assert!(asls_fit.baseline.iter().all(|value| value.is_finite()));
        assert!(arpls_fit.baseline.iter().all(|value| value.is_finite()));
    }

    #[test]
    fn non_positive_tolerance_forces_max_iterations() {
        let y: Vec<f64> = (0..100)
            .map(|i| {
                let x = i as f64 / 99.0;
                0.5 + 0.2 * x + (-(x - 0.45).powi(2) / 0.002).exp()
            })
            .collect();

        let fit = asls(
            &y,
            AslsParams {
                whittaker: crate::whittaker::WhittakerParams {
                    tol: -1.0,
                    max_iter: 3,
                    ..crate::whittaker::WhittakerParams::default()
                },
                p: 0.01,
            },
        )
        .unwrap();

        assert_eq!(fit.report.iterations, 3);
        assert!(!fit.report.converged);
    }
}