✨ feat: nan handling (#59)

* ✨ feat: add nan implementation of m4 algorithm

* ✨ feat: add nan implementation of minmax algorithm

* ✨ feat: add nan implementation of minmaxlttb algorithm

* 💩 feat: update lib script to incorporate nan-handling functions

* 🚧 feat: add new nan downsampler

* ✅ tests: add new nan functions to rust mod tests

* 🎨 chore: format code

* ✨ feat: expose new nan downsamplers to api

* ✅ tests: update tsdownsample tests to support nan downsamplers

* 🎨 chore: format code

* ✅ tests: add test for nan downsamplers

* ✨ feat: add python counterparts of Rust downsamplers

* ✅ tests: re-enable commented out tests

* 🎨 chore: format code

* 🔥 chore: remove commented code

* 📝 docs: update README.md

* 📝 docs: update NaN descriptions

* 🧹 remove threaded

* 🎉 cleanup code

* 🙈 fix typo in NaNMinMaxDownsampler

* 🕵️ benchmark NaN downsamplers

* 🧹

* 🧹

* 🧹 limit duplicate code

* 🙈 fix linting

* 🧹

---------

Co-authored-by: Jeroen Van Der Donckt <[email protected]>
Co-authored-by: jvdd <[email protected]>

Author: NielsPraet

README.md

@@ -6,13 +6,14 @@
 [![CodeQL](https://github.com/predict-idlab/tsdownsample/actions/workflows/codeql.yml/badge.svg)](https://github.com/predict-idlab/tsdownsample/actions/workflows/codeql.yml)
 [![Testing](https://github.com/predict-idlab/tsdownsample/actions/workflows/ci-downsample_rs.yml/badge.svg)](https://github.com/predict-idlab/tsdownsample/actions/workflows/ci-downsample_rs.yml)
 [![Testing](https://github.com/predict-idlab/tsdownsample/actions/workflows/ci-tsdownsample.yml/badge.svg)](https://github.com/predict-idlab/tsdownsample/actions/workflows/ci-tsdownsample.yml)
+
 <!-- TODO: codecov -->
 
 Extremely fast **time series downsampling 📈** for visualization, written in Rust.
 
 ## Features ✨
 
-* **Fast**: written in rust with PyO3 bindings
+- **Fast**: written in rust with PyO3 bindings
   - leverages optimized [argminmax](https://github.com/jvdd/argminmax) - which is SIMD accelerated with runtime feature detection
   - scales linearly with the number of data points
   <!-- TODO check if it scales sublinearly -->
@@ -25,21 +26,21 @@ Extremely fast **time series downsampling 📈** for visualization, written in R
       </blockquote>
       In Rust - which is a compiled language - there is no GIL, so CPU-bound tasks can be parallelized (with <a href="https://github.com/rayon-rs/rayon">Rayon</a>) with little to no overhead.
     </details>
-* **Efficient**: memory efficient
+- **Efficient**: memory efficient
   - works on views of the data (no copies)
   - no intermediate data structures are created
-* **Flexible**: works on any type of data
-    - supported datatypes are 
-      - for `x`: `f32`, `f64`, `i16`, `i32`, `i64`, `u16`, `u32`, `u64`, `datetime64`, `timedelta64`
-      - for `y`: `f16`, `f32`, `f64`, `i8`, `i16`, `i32`, `i64`, `u8`, `u16`, `u32`, `u64`, `datetime64`, `timedelta64`, `bool`
+- **Flexible**: works on any type of data
+  - supported datatypes are
+    - for `x`: `f32`, `f64`, `i16`, `i32`, `i64`, `u16`, `u32`, `u64`, `datetime64`, `timedelta64`
+    - for `y`: `f16`, `f32`, `f64`, `i8`, `i16`, `i32`, `i64`, `u8`, `u16`, `u32`, `u64`, `datetime64`, `timedelta64`, `bool`
     <details>
       <summary><i>!! 🚀 <code>f16</code> <a href="https://github.com/jvdd/argminmax">argminmax</a> is 200-300x faster than numpy</i></summary>
       In contrast with all other data types above, <code>f16</code> is *not* hardware supported (i.e., no instructions for f16) by most modern CPUs!! <br>
       🐌 Programming languages facilitate support for this datatype by either (i) upcasting to <u>f32</u> or (ii) using a software implementation. <br>
       💡 As for argminmax, only comparisons are needed - and thus no arithmetic operations - creating a <u>symmetrical ordinal mapping from <code>f16</code> to <code>i16</code></u> is sufficient. This mapping allows to use the hardware supported scalar and SIMD <code>i16</code> instructions - while not producing any memory overhead 🎉 <br>
       <i>More details are described in <a href="https://github.com/jvdd/argminmax/pull/1">argminmax PR #1</a>.</i>
     </details>
-* **Easy to use**: simple & flexible API
+- **Easy to use**: simple & flexible API
 
 ## Install
 
@@ -83,6 +84,7 @@ downsample([x], y, n_out, **kwargs) -> ndarray[uint64]
 ```
 
 **Arguments**:
+
 - `x` is optional
 - `x` and `y` are both positional arguments
 - `n_out` is a mandatory keyword argument that defines the number of output values<sup>*</sup>
@@ -93,7 +95,8 @@ downsample([x], y, n_out, **kwargs) -> ndarray[uint64]
 
 **Returns**: a `ndarray[uint64]` of indices that can be used to index the original data.
 
-<sup>*</sup><i>When there are gaps in the time series, fewer than `n_out` indices may be returned.</i>
+<sup>\*</sup><i>When there are gaps in the time series, fewer than `n_out` indices may be returned.</i>
+
 ### Downsampling algorithms 📈
 
 The following downsampling algorithms (classes) are implemented:
@@ -107,12 +110,28 @@ The following downsampling algorithms (classes) are implemented:
 
 <sup>*</sup><i>Default value for `minmax_ratio` is 4, which is empirically proven to be a good default. More details here: https://arxiv.org/abs/2305.00332</i>
 
+### Handling NaNs
+
+This library supports two `NaN`-policies:
+
+1. Omit `NaN`s (`NaN`s are ignored during downsampling).
+2. Return index of first `NaN` once there is at least one present in the bin of the considered data.
+
+|             Omit `NaN`s | Return `NaN`s              |
+| ----------------------: | :------------------------- |
+|     `MinMaxDownsampler` | `NaNMinMaxDownsampler`     |
+|         `M4Downsampler` | `NaNM4Downsampler`         |
+| `MinMaxLTTBDownsampler` | `NaNMinMaxLTTBDownsampler` |
+|       `LTTBDownsampler` |                            |
+
+> Note that NaNs are not supported for `x`-data.
 
 ## Limitations & assumptions 🚨
 
 Assumes;
+
 1. `x`-data is (non-strictly) monotonic increasing (i.e., sorted)
-2. no `NaNs` in the data
+2. no `NaN`s in `x`-data
 
 ---
 

downsample_rs/src/m4.rs

@@ -1,4 +1,4 @@
-use argminmax::ArgMinMax;
+use argminmax::{ArgMinMax, NaNArgMinMax};
 use num_traits::{AsPrimitive, FromPrimitive};
 use rayon::iter::IndexedParallelIterator;
 use rayon::prelude::*;
@@ -13,55 +13,82 @@ use super::POOL;
 
 // ----------- WITH X
 
-pub fn m4_with_x<Tx, Ty>(x: &[Tx], arr: &[Ty], n_out: usize) -> Vec<usize>
-where
-    for<'a> &'a [Ty]: ArgMinMax,
-    Tx: Num + FromPrimitive + AsPrimitive<f64>,
-    Ty: Copy + PartialOrd,
-{
-    assert_eq!(n_out % 4, 0);
-    let bin_idx_iterator = get_equidistant_bin_idx_iterator(x, n_out / 4);
-    m4_generic_with_x(arr, bin_idx_iterator, n_out, |arr| arr.argminmax())
+macro_rules! m4_with_x {
+    ($func_name:ident, $trait:path, $f_argminmax:expr) => {
+        pub fn $func_name<Tx, Ty>(x: &[Tx], arr: &[Ty], n_out: usize) -> Vec<usize>
+        where
+            for<'a> &'a [Ty]: $trait,
+            Tx: Num + FromPrimitive + AsPrimitive<f64>,
+            Ty: Copy + PartialOrd,
+        {
+            assert_eq!(n_out % 4, 0);
+            let bin_idx_iterator = get_equidistant_bin_idx_iterator(x, n_out / 4);
+            m4_generic_with_x(arr, bin_idx_iterator, n_out, $f_argminmax)
+        }
+    };
 }
 
+m4_with_x!(m4_with_x, ArgMinMax, |arr| arr.argminmax());
+m4_with_x!(m4_with_x_nan, NaNArgMinMax, |arr| arr.nanargminmax());
+
 // ----------- WITHOUT X
 
-pub fn m4_without_x<T: Copy + PartialOrd>(arr: &[T], n_out: usize) -> Vec<usize>
-where
-    for<'a> &'a [T]: ArgMinMax,
-{
-    assert_eq!(n_out % 4, 0);
-    m4_generic(arr, n_out, |arr| arr.argminmax())
+macro_rules! m4_without_x {
+    ($func_name:ident, $trait:path, $f_argminmax:expr) => {
+        pub fn $func_name<T: Copy + PartialOrd>(arr: &[T], n_out: usize) -> Vec<usize>
+        where
+            for<'a> &'a [T]: $trait,
+        {
+            assert_eq!(n_out % 4, 0);
+            m4_generic(arr, n_out, $f_argminmax)
+        }
+    };
 }
 
+m4_without_x!(m4_without_x, ArgMinMax, |arr| arr.argminmax());
+m4_without_x!(m4_without_x_nan, NaNArgMinMax, |arr| arr.nanargminmax());
+
 // ------------------------------------- PARALLEL --------------------------------------
 
 // ----------- WITH X
 
-pub fn m4_with_x_parallel<Tx, Ty>(x: &[Tx], arr: &[Ty], n_out: usize) -> Vec<usize>
-where
-    for<'a> &'a [Ty]: ArgMinMax,
-    Tx: Num + FromPrimitive + AsPrimitive<f64> + Send + Sync,
-    Ty: Copy + PartialOrd + Send + Sync,
-{
-    assert_eq!(n_out % 4, 0);
-    let bin_idx_iterator = get_equidistant_bin_idx_iterator_parallel(x, n_out / 4);
-    m4_generic_with_x_parallel(arr, bin_idx_iterator, n_out, |arr| arr.argminmax())
+macro_rules! m4_with_x_parallel {
+    ($func_name:ident, $trait:path, $f_argminmax:expr) => {
+        pub fn $func_name<Tx, Ty>(x: &[Tx], arr: &[Ty], n_out: usize) -> Vec<usize>
+        where
+            for<'a> &'a [Ty]: $trait,
+            Tx: Num + FromPrimitive + AsPrimitive<f64> + Send + Sync,
+            Ty: Copy + PartialOrd + Send + Sync,
+        {
+            assert_eq!(n_out % 4, 0);
+            let bin_idx_iterator = get_equidistant_bin_idx_iterator_parallel(x, n_out / 4);
+            m4_generic_with_x_parallel(arr, bin_idx_iterator, n_out, $f_argminmax)
+        }
+    };
 }
 
+m4_with_x_parallel!(m4_with_x_parallel, ArgMinMax, |arr| arr.argminmax());
+m4_with_x_parallel!(m4_with_x_parallel_nan, NaNArgMinMax, |arr| arr
+    .nanargminmax());
+
 // ----------- WITHOUT X
 
-pub fn m4_without_x_parallel<T: Copy + PartialOrd + Send + Sync>(
-    arr: &[T],
-    n_out: usize,
-) -> Vec<usize>
-where
-    for<'a> &'a [T]: ArgMinMax,
-{
-    assert_eq!(n_out % 4, 0);
-    m4_generic_parallel(arr, n_out, |arr| arr.argminmax())
+macro_rules! m4_without_x_parallel {
+    ($func_name:ident, $trait:path, $f_argminmax:expr) => {
+        pub fn $func_name<T: Copy + PartialOrd + Send + Sync>(arr: &[T], n_out: usize) -> Vec<usize>
+        where
+            for<'a> &'a [T]: $trait,
+        {
+            assert_eq!(n_out % 4, 0);
+            m4_generic_parallel(arr, n_out, $f_argminmax)
+        }
+    };
 }
 
+m4_without_x_parallel!(m4_without_x_parallel, ArgMinMax, |arr| arr.argminmax());
+m4_without_x_parallel!(m4_without_x_parallel_nan, NaNArgMinMax, |arr| arr
+    .nanargminmax());
+
 // TODO: check for duplicate data in the output array
 // -> In the current implementation we always add 4 datapoints per bin (if of
 //    course the bin has >= 4 datapoints). However, the argmin and argmax might

downsample_rs/src/minmax.rs

@@ -1,7 +1,7 @@
 use rayon::iter::IndexedParallelIterator;
 use rayon::prelude::*;
 
-use argminmax::ArgMinMax;
+use argminmax::{ArgMinMax, NaNArgMinMax};
 use num_traits::{AsPrimitive, FromPrimitive};
 
 use super::searchsorted::{
@@ -14,55 +14,83 @@ use super::POOL;
 
 // ----------- WITH X
 
-pub fn min_max_with_x<Tx, Ty>(x: &[Tx], arr: &[Ty], n_out: usize) -> Vec<usize>
-where
-    for<'a> &'a [Ty]: ArgMinMax,
-    Tx: Num + FromPrimitive + AsPrimitive<f64>,
-    Ty: Copy + PartialOrd,
-{
-    assert_eq!(n_out % 2, 0);
-    let bin_idx_iterator = get_equidistant_bin_idx_iterator(x, n_out / 2);
-    min_max_generic_with_x(arr, bin_idx_iterator, n_out, |arr| arr.argminmax())
+macro_rules! min_max_with_x {
+    ($func_name:ident, $trait:path, $f_argminmax:expr) => {
+        pub fn $func_name<Tx, Ty>(x: &[Tx], arr: &[Ty], n_out: usize) -> Vec<usize>
+        where
+            for<'a> &'a [Ty]: $trait,
+            Tx: Num + FromPrimitive + AsPrimitive<f64>,
+            Ty: Copy + PartialOrd,
+        {
+            assert_eq!(n_out % 2, 0);
+            let bin_idx_iterator = get_equidistant_bin_idx_iterator(x, n_out / 2);
+            min_max_generic_with_x(arr, bin_idx_iterator, n_out, $f_argminmax)
+        }
+    };
 }
 
+min_max_with_x!(min_max_with_x, ArgMinMax, |arr| arr.argminmax());
+min_max_with_x!(min_max_with_x_nan, NaNArgMinMax, |arr| arr.nanargminmax());
+
 // ----------- WITHOUT X
 
-pub fn min_max_without_x<T: Copy + PartialOrd>(arr: &[T], n_out: usize) -> Vec<usize>
-where
-    for<'a> &'a [T]: ArgMinMax,
-{
-    assert_eq!(n_out % 2, 0);
-    min_max_generic(arr, n_out, |arr| arr.argminmax())
+macro_rules! min_max_without_x {
+    ($func_name:ident, $trait:path, $f_argminmax:expr) => {
+        pub fn $func_name<T: Copy + PartialOrd>(arr: &[T], n_out: usize) -> Vec<usize>
+        where
+            for<'a> &'a [T]: $trait,
+        {
+            assert_eq!(n_out % 2, 0);
+            min_max_generic(arr, n_out, $f_argminmax)
+        }
+    };
 }
 
+min_max_without_x!(min_max_without_x, ArgMinMax, |arr| arr.argminmax());
+min_max_without_x!(min_max_without_x_nan, NaNArgMinMax, |arr| arr
+    .nanargminmax());
+
 // ------------------------------------- PARALLEL --------------------------------------
 
 // ----------- WITH X
 
-pub fn min_max_with_x_parallel<Tx, Ty>(x: &[Tx], arr: &[Ty], n_out: usize) -> Vec<usize>
-where
-    for<'a> &'a [Ty]: ArgMinMax,
-    Tx: Num + FromPrimitive + AsPrimitive<f64> + Send + Sync,
-    Ty: Copy + PartialOrd + Send + Sync,
-{
-    assert_eq!(n_out % 2, 0);
-    let bin_idx_iterator = get_equidistant_bin_idx_iterator_parallel(x, n_out / 2);
-    min_max_generic_with_x_parallel(arr, bin_idx_iterator, n_out, |arr| arr.argminmax())
+macro_rules! min_max_with_x_parallel {
+    ($func_name:ident, $trait:path, $f_argminmax:expr) => {
+        pub fn $func_name<Tx, Ty>(x: &[Tx], arr: &[Ty], n_out: usize) -> Vec<usize>
+        where
+            for<'a> &'a [Ty]: $trait,
+            Tx: Num + FromPrimitive + AsPrimitive<f64> + Send + Sync,
+            Ty: Copy + PartialOrd + Send + Sync,
+        {
+            assert_eq!(n_out % 2, 0);
+            let bin_idx_iterator = get_equidistant_bin_idx_iterator_parallel(x, n_out / 2);
+            min_max_generic_with_x_parallel(arr, bin_idx_iterator, n_out, $f_argminmax)
+        }
+    };
 }
 
+min_max_with_x_parallel!(min_max_with_x_parallel, ArgMinMax, |arr| arr.argminmax());
+min_max_with_x_parallel!(min_max_with_x_parallel_nan, NaNArgMinMax, |arr| arr
+    .nanargminmax());
+
 // ----------- WITHOUT X
 
-pub fn min_max_without_x_parallel<T: Copy + PartialOrd + Send + Sync>(
-    arr: &[T],
-    n_out: usize,
-) -> Vec<usize>
-where
-    for<'a> &'a [T]: ArgMinMax,
-{
-    assert_eq!(n_out % 2, 0);
-    min_max_generic_parallel(arr, n_out, |arr| arr.argminmax())
+macro_rules! min_max_without_x_parallel {
+    ($func_name:ident, $trait:path, $f_argminmax:expr) => {
+        pub fn $func_name<T: Copy + PartialOrd + Send + Sync>(arr: &[T], n_out: usize) -> Vec<usize>
+        where
+            for<'a> &'a [T]: $trait,
+        {
+            assert_eq!(n_out % 2, 0);
+            min_max_generic_parallel(arr, n_out, $f_argminmax)
+        }
+    };
 }
 
+min_max_without_x_parallel!(min_max_without_x_parallel, ArgMinMax, |arr| arr.argminmax());
+min_max_without_x_parallel!(min_max_without_x_parallel_nan, NaNArgMinMax, |arr| arr
+    .nanargminmax());
+
 // ----------------------------------- GENERICS ------------------------------------
 
 // --------------------- WITHOUT X