light_curve_feature/data/
data_sample.rs

1use crate::data::sorted_array::SortedArray;
2use crate::float_trait::Float;
3use crate::types::CowArray1;
4
5use conv::prelude::*;
6use ndarray::{Array1, ArrayView1, Zip, s};
7
8/// A [`TimeSeries`] component
9#[derive(Clone, Debug)]
10pub struct DataSample<'a, T>
11where
12    T: Float,
13{
14    pub sample: CowArray1<'a, T>,
15    sorted: Option<SortedArray<T>>,
16    min: Option<T>,
17    max: Option<T>,
18    mean: Option<T>,
19    median: Option<T>,
20    std: Option<T>,
21    std2: Option<T>,
22}
23
24impl<'a, T> PartialEq for DataSample<'a, T>
25where
26    T: Float,
27{
28    fn eq(&self, other: &Self) -> bool {
29        self.sample == other.sample
30    }
31}
32
33macro_rules! data_sample_getter {
34    ($attr: ident, $getter: ident, $func: expr, $method_sorted: ident) => {
35        // This lint is false-positive in macros
36        // https://github.com/rust-lang/rust-clippy/issues/1553
37        #[allow(clippy::redundant_closure_call)]
38        pub fn $getter(&mut self) -> T {
39            match self.$attr {
40                Some(x) => x,
41                None => {
42                    self.$attr = Some(match self.sorted.as_ref() {
43                        Some(sorted) => sorted.$method_sorted(),
44                        None => $func(self),
45                    });
46                    self.$attr.unwrap()
47                }
48            }
49        }
50    };
51    ($attr: ident, $getter: ident, $func: expr) => {
52        // This lint is false-positive in macros
53        // https://github.com/rust-lang/rust-clippy/issues/1553
54        #[allow(clippy::redundant_closure_call)]
55        pub fn $getter(&mut self) -> T {
56            match self.$attr {
57                Some(x) => x,
58                None => {
59                    self.$attr = Some($func(self));
60                    self.$attr.unwrap()
61                }
62            }
63        }
64    };
65}
66
67impl<'a, T> DataSample<'a, T>
68where
69    T: Float,
70{
71    pub fn new(sample: CowArray1<'a, T>) -> Self {
72        Self {
73            sample,
74            sorted: None,
75            min: None,
76            max: None,
77            mean: None,
78            median: None,
79            std: None,
80            std2: None,
81        }
82    }
83
84    pub fn as_slice(&mut self) -> &[T] {
85        if !self.sample.is_standard_layout() {
86            let owned: Array1<_> = self.sample.iter().copied().collect::<Vec<_>>().into();
87            self.sample = owned.into();
88        }
89        self.sample.as_slice().unwrap()
90    }
91
92    pub fn get_sorted(&mut self) -> &SortedArray<T> {
93        if self.sorted.is_none() {
94            self.sorted = Some(self.sample.to_vec().into());
95        }
96        self.sorted.as_ref().unwrap()
97    }
98
99    fn set_min_max(&mut self) {
100        let (min, max) =
101            self.sample
102                .slice(s![1..])
103                .fold((self.sample[0], self.sample[0]), |(min, max), &x| {
104                    if x > max {
105                        (min, x)
106                    } else if x < min {
107                        (x, max)
108                    } else {
109                        (min, max)
110                    }
111                });
112        self.min = Some(min);
113        self.max = Some(max);
114    }
115
116    data_sample_getter!(
117        min,
118        get_min,
119        |ds: &mut DataSample<'a, T>| {
120            ds.set_min_max();
121            ds.min.unwrap()
122        },
123        minimum
124    );
125    data_sample_getter!(
126        max,
127        get_max,
128        |ds: &mut DataSample<'a, T>| {
129            ds.set_min_max();
130            ds.max.unwrap()
131        },
132        maximum
133    );
134    data_sample_getter!(mean, get_mean, |ds: &mut DataSample<'a, T>| {
135        ds.sample.mean().expect("time series must be non-empty")
136    });
137    data_sample_getter!(median, get_median, |ds: &mut DataSample<'a, T>| {
138        ds.get_sorted().median()
139    });
140    data_sample_getter!(std, get_std, |ds: &mut DataSample<'a, T>| {
141        ds.get_std2().sqrt()
142    });
143    data_sample_getter!(std2, get_std2, |ds: &mut DataSample<'a, T>| {
144        // Benchmarks show that it is faster than `ndarray::ArrayBase::var(T::one)`
145        let mean = ds.get_mean();
146        ds.sample
147            .fold(T::zero(), |sum, &x| sum + (x - mean).powi(2))
148            / (ds.sample.len() - 1).approx().unwrap()
149    });
150
151    pub fn signal_to_noise(&mut self, value: T) -> T {
152        if self.get_std().is_zero() {
153            T::zero()
154        } else {
155            (value - self.get_mean()) / self.get_std()
156        }
157    }
158
159    /// Returns true if all values are equal. Always true for zero- or one- length
160    pub fn is_all_same(&self) -> bool {
161        if self.sample.is_empty() {
162            return true;
163        }
164        if self.max.is_some() && self.max == self.min {
165            return true;
166        }
167        if self.std2 == Some(T::zero()) {
168            return true;
169        }
170        if let Some(sorted) = &self.sorted {
171            return sorted[0] == sorted[sorted.len() - 1];
172        }
173        let x0 = self.sample[0];
174        // all() returns true for the empty slice, i.e. single-point time series
175        Zip::from(self.sample.slice(s![1..])).all(|&x| x == x0)
176    }
177}
178
179impl<'a, T> From<SortedArray<T>> for DataSample<'a, T>
180where
181    T: Float,
182{
183    fn from(sorted: SortedArray<T>) -> Self {
184        let sample = sorted.0.clone().into();
185        Self {
186            sample,
187            sorted: Some(sorted),
188            min: None,
189            max: None,
190            median: None,
191            mean: None,
192            std: None,
193            std2: None,
194        }
195    }
196}
197
198impl<'a, T, Slice: ?Sized> From<&'a Slice> for DataSample<'a, T>
199where
200    T: Float,
201    Slice: AsRef<[T]>,
202{
203    fn from(s: &'a Slice) -> Self {
204        ArrayView1::from(s).into()
205    }
206}
207
208impl<'a, T> From<Vec<T>> for DataSample<'a, T>
209where
210    T: Float,
211{
212    fn from(v: Vec<T>) -> Self {
213        Array1::from(v).into()
214    }
215}
216
217impl<'a, T> From<ArrayView1<'a, T>> for DataSample<'a, T>
218where
219    T: Float,
220{
221    fn from(a: ArrayView1<'a, T>) -> Self {
222        Self::new(a.into())
223    }
224}
225
226impl<'a, T> From<Array1<T>> for DataSample<'a, T>
227where
228    T: Float,
229{
230    fn from(a: Array1<T>) -> Self {
231        Self::new(a.into())
232    }
233}
234
235impl<'a, T> From<CowArray1<'a, T>> for DataSample<'a, T>
236where
237    T: Float,
238{
239    fn from(a: CowArray1<'a, T>) -> Self {
240        Self::new(a)
241    }
242}
243
244#[cfg(test)]
245#[allow(clippy::unreadable_literal)]
246#[allow(clippy::excessive_precision)]
247mod tests {
248    use super::*;
249
250    use approx::assert_relative_eq;
251
252    macro_rules! data_sample_test {
253        ($name: ident, $method: ident, $desired: literal, $x: tt $(,)?) => {
254            #[test]
255            fn $name() {
256                let x = $x;
257                let desired = $desired;
258
259                let mut ds: DataSample<_> = DataSample::from(&x);
260                assert_relative_eq!(ds.$method(), desired, epsilon = 1e-6);
261                assert_relative_eq!(ds.$method(), desired, epsilon = 1e-6);
262
263                let mut ds: DataSample<_> = DataSample::from(&x);
264                ds.get_sorted();
265                assert_relative_eq!(ds.$method(), desired, epsilon = 1e-6);
266                assert_relative_eq!(ds.$method(), desired, epsilon = 1e-6);
267            }
268        };
269    }
270
271    data_sample_test!(
272        data_sample_min,
273        get_min,
274        -7.79420906,
275        [3.92948846, 3.28436964, 6.73375373, -7.79420906, -7.23407407],
276    );
277
278    data_sample_test!(
279        data_sample_max,
280        get_max,
281        6.73375373,
282        [3.92948846, 3.28436964, 6.73375373, -7.79420906, -7.23407407],
283    );
284
285    data_sample_test!(
286        data_sample_mean,
287        get_mean,
288        -0.21613426,
289        [3.92948846, 3.28436964, 6.73375373, -7.79420906, -7.23407407],
290    );
291
292    data_sample_test!(
293        data_sample_median_odd,
294        get_median,
295        3.28436964,
296        [3.92948846, 3.28436964, 6.73375373, -7.79420906, -7.23407407],
297    );
298
299    data_sample_test!(
300        data_sample_median_even,
301        get_median,
302        5.655794743124782,
303        [
304            9.47981408, 3.86815751, 9.90299294, -2.986894, 7.44343197, 1.52751816
305        ],
306    );
307
308    data_sample_test!(
309        data_sample_std,
310        get_std,
311        6.7900544035968435,
312        [3.92948846, 3.28436964, 6.73375373, -7.79420906, -7.23407407],
313    );
314
315    /// https://github.com/light-curve/light-curve-feature/issues/95
316    #[test]
317    fn std2_overflow() {
318        const N: usize = (1 << 24) + 2;
319        // Such a large integer cannot be represented as a float32
320        let x = Array1::linspace(0.0_f32, 1.0, N);
321        let mut ds = DataSample::new(x.into());
322        // This should not panic
323        let _std2 = ds.get_std2();
324    }
325
326    // signal_to_noise measures how many std deviations a value is from the mean
327    #[test]
328    fn signal_to_noise_correct() {
329        // mean=3, std of [1,2,3,4,5] = sqrt(2.5) ≈ 1.5811
330        let x = [1.0_f64, 2.0, 3.0, 4.0, 5.0];
331        let mut ds: DataSample<f64> = DataSample::from(&x);
332        let snr = ds.signal_to_noise(5.0);
333        let expected = (5.0 - 3.0) / f64::sqrt(2.5);
334        assert_relative_eq!(snr, expected, epsilon = 1e-10);
335    }
336
337    // signal_to_noise returns 0 for flat data (zero std) regardless of input value
338    #[test]
339    fn signal_to_noise_zero_for_flat_data() {
340        let x = [7.0_f64; 5];
341        let mut ds: DataSample<f64> = DataSample::from(&x);
342        assert_eq!(ds.signal_to_noise(42.0), 0.0);
343    }
344
345    // is_all_same is true for empty array
346    #[test]
347    fn is_all_same_true_for_empty() {
348        let x: &[f64] = &[];
349        let ds: DataSample<f64> = DataSample::from(x);
350        assert!(ds.is_all_same());
351    }
352
353    // is_all_same is true for single-element array
354    #[test]
355    fn is_all_same_true_for_single_element() {
356        let x = [std::f64::consts::PI];
357        let ds: DataSample<f64> = DataSample::from(&x);
358        assert!(ds.is_all_same());
359    }
360
361    // is_all_same is true for constant array
362    #[test]
363    fn is_all_same_true_for_constant() {
364        let x = [5.0_f64; 10];
365        let ds: DataSample<f64> = DataSample::from(&x);
366        assert!(ds.is_all_same());
367    }
368
369    // is_all_same is false for array with at least two distinct values
370    #[test]
371    fn is_all_same_false_for_varying() {
372        let x = [1.0_f64, 1.0, 2.0, 1.0];
373        let ds: DataSample<f64> = DataSample::from(&x);
374        assert!(!ds.is_all_same());
375    }
376
377    // get_sorted produces sorted order regardless of input order
378    #[test]
379    fn get_sorted_is_monotonically_nondecreasing() {
380        let x = [5.0_f64, 1.0, 3.0, 2.0, 4.0];
381        let mut ds: DataSample<f64> = DataSample::from(&x);
382        let sorted = ds.get_sorted();
383        for i in 1..sorted.len() {
384            assert!(
385                sorted[i - 1] <= sorted[i],
386                "sorted array not monotone at index {i}"
387            );
388        }
389    }
390
391    // get_sorted is idempotent: min/max/median queried on sorted data give same results
392    #[test]
393    fn sorted_and_unsorted_statistics_agree() {
394        let x = [5.0_f64, 1.0, 3.0, 2.0, 4.0];
395
396        let mut ds_unsorted: DataSample<f64> = DataSample::from(&x);
397        let mut ds_sorted: DataSample<f64> = DataSample::from(&x);
398        ds_sorted.get_sorted(); // pre-populate sorted cache
399
400        assert_relative_eq!(ds_unsorted.get_min(), ds_sorted.get_min(), epsilon = 1e-10);
401        assert_relative_eq!(ds_unsorted.get_max(), ds_sorted.get_max(), epsilon = 1e-10);
402        assert_relative_eq!(
403            ds_unsorted.get_mean(),
404            ds_sorted.get_mean(),
405            epsilon = 1e-10
406        );
407        assert_relative_eq!(
408            ds_unsorted.get_median(),
409            ds_sorted.get_median(),
410            epsilon = 1e-10
411        );
412        assert_relative_eq!(ds_unsorted.get_std(), ds_sorted.get_std(), epsilon = 1e-10);
413    }
414}
light_curve_feature/data/data_sample.rs

light_curve_feature/data/
data_sample.rs
