Repository: AnBowell/whittaker-eilers (language: python)

Issue:
Maybe add batch processing?
Thank you for developing this package! I am currently using it with Python and encountered a limitation when trying to process a batch of signals. The `whittaker_smoother.smooth()` function only accepts 1-D arrays. Maybe add such functionality to process `(batch_size, signal_len)` input? 


Discussion / hints:
I like the idea! 

I think if I were to do this, I'd end up adding some form of multithreading to it too as otherwise it would just be the equivalent of looping over the arrays. I'll have a think about how it could be structured! 

Fix patch:
diff --git a/Cargo.toml b/Cargo.toml
index f77bf59..6fcc0c1 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "whittaker-eilers"
-version = "0.1.3"
+version = "0.2.0"
 edition = "2021"
 
 description = "A sparse matrix implementation of Whittaker-Eilers smoothing and interpolation"
@@ -16,12 +16,15 @@ readme = "README.md"
 
 
 [dependencies]
-sprs = "0.11.1"
+
+sprs = "0.11.2"
 sprs-ldl = "0.10.0"
-nalgebra = "0.32.3"
+nalgebra = "0.33.2"
+rayon = { version = "1.10.0", optional = true }
+
 
 [dev-dependencies]
-plotly = "0.8.4"
+plotly = "0.10.0"
 approx = "0.5.1"
 criterion = "0.5.1"
 rand = "0.8.5"
@@ -34,3 +37,6 @@ harness = false
 
 [workspace]
 members = ["whittaker-eilers-py"]
+
+[features]
+rayon = ["dep:rayon"]
diff --git a/README.md b/README.md
index feb55a8..70a37b1 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,7 @@
 [![Downloads](https://static.pepy.tech/badge/whittaker-eilers)](https://pepy.tech/project/whittaker-eilers) [![Downloads](https://static.pepy.tech/badge/whittaker-eilers/month)](https://pepy.tech/project/whittaker-eilers)
+
 # Whittaker-Eilers Smoothing and Interpolation
+
 **The Whittaker-Eilers smoother is the perfect smoother.** It offers extremely quick, efficient smoothing with built-in interpolation via weights on each measurement. This crate provides a sparse-matrix implementation for additional speed and memory efficiency and can handle both equally and unequally spaced measurements.
 
 ---
@@ -10,16 +12,19 @@ whittaker-eilers = "0.1.3"
 ```
 
 ## Usage
+
 To start smoothing and interpolating data, create a reusable WhittakerSmoother struct via the `new` function. You'll only need to recreate this struct if the length or sampling rate of your data changes.
 
 ### Equally spaced data
+
 This is the fastest smoothing option. It smooths equally spaced y measurements using two tunable parameters, `lambda` (2e4) and the smoother `order` (2). The larger the lambda, the smoother the data.
+
 ```rust
 use whittaker_eilers::WhittakerSmoother;
 
 let data_to_smooth = vec![1.1, 1.9, 3.1, 3.91, 5.0, 6.02, 7.01, 7.7, 9.0, 10.0];
 
-let whittaker_smoother = 
+let whittaker_smoother =
             WhittakerSmoother::new(2e4, 2, data_to_smooth.len(), None, None)
             .unwrap();
 
@@ -27,17 +32,17 @@ let smoothed_data = whittaker_smoother.smooth(&data_to_smooth).unwrap();
 println!("Smoothed data: {:?}", smoothed_data);
 ```
 
+### Non-equally spaced data
 
+If you wish to smooth unequally spaced data, you need to provide an `x_input` with the sample times/positions.
 
-### Non-equally spaced data
-If you wish to smooth unequally spaced data, you need to provide an `x_input` with the sample times/positions. 
 ```rust
 use whittaker_eilers::WhittakerSmoother;
 
 let x_input = vec![1.1, 1.9, 3.1, 3.91, 5.0, 6.02, 7.01, 7.7, 9.0, 10.0];
 let data_to_smooth = vec![1.1, 1.9, 3.1, 3.91, 5.0, 6.02, 7.01, 7.7, 9.0, 10.0];
 
-let whittaker_smoother = 
+let whittaker_smoother =
             WhittakerSmoother::new(2e4, 2, data_to_smooth.len(), Some(&x_input), None)
             .unwrap();
 
@@ -48,7 +53,9 @@ println!("Smoothed data: {:?}", smoothed_data);
 ```
 
 ### Weighted data & Interpolation
-Each measurement can then be weighted to trust some measurements more than others. Setting `weights` to 0 for measurements will lead to interpolation. 
+
+Each measurement can then be weighted to trust some measurements more than others. Setting `weights` to 0 for measurements will lead to interpolation.
+
 ```rust
 use whittaker_eilers::WhittakerSmoother;
 
@@ -67,18 +74,17 @@ println!("Smoothed data: {:?}", smoothed_data);
 
 ```
 
-
 ### Smoothing with cross validation
-With this package, you can also calculate the cross validation error alongside the smoothed series. This shouldn't really be used in production where speed is necessary though!
 
+With this package, you can also calculate the cross validation error alongside the smoothed series. This shouldn't really be used in production where speed is necessary though!
 
-```rust 
+```rust
 use whittaker_eilers::WhittakerSmoother;
 
 let x_input = vec![1.1, 1.9, 3.1, 3.91, 5.0, 6.02, 7.01, 7.7, 9.0, 10.0 ,11.0, 12.0, 13.0];
 let data_to_smooth = vec![1.1, 1.9, 3.1, 3.91, 5.0, 6.02, 7.01, 7.7, 9.0, 10.0, 11.0, 12.0, 13.0];
 
-let whittaker_smoother = 
+let whittaker_smoother =
             WhittakerSmoother::new(2e4, 2, data_to_smooth.len(), Some(&x_input), None)
             .unwrap();
 
@@ -88,17 +94,17 @@ println!("Result: {:?}", smoothed_and_cross_validated);
 ```
 
 ### Automatic smoothing
+
 Smoothing data requires a choice of Lambda. This can be done using visual inspection or by finding the lambda
 which results in the lowest cross validation error. The `smooth_optimal` function runs the smoother for a variety of lambdas and returns the results with the ability to retrieve the optimal one.
 
-
-```rust 
+```rust
 use whittaker_eilers::WhittakerSmoother;
 
 let x_input = vec![1.1, 1.9, 3.1, 3.91, 5.0, 6.02, 7.01, 7.7, 9.0, 10.0 ,11.0, 12.0, 13.0];
 let data_to_smooth = vec![1.1, 1.9, 3.1, 3.91, 5.0, 6.02, 7.01, 7.7, 9.0, 10.0, 11.0, 12.0, 13.0];
 
-let mut whittaker_smoother = 
+let mut whittaker_smoother =
             WhittakerSmoother::new(2e4, 2, data_to_smooth.len(), Some(&x_input), None)
             .unwrap();
 
@@ -110,26 +116,48 @@ println!("Optimal result: {:?}", results.get_optimal());
 
 ```
 
+### Parallel Smoothing
+
+Using the feature `rayon` enables this package to smooth multiple data series in parallel. While y, the data to be smoothed can vary, x, weights, and the other smoothing parameters cannot be varied across the different data series.
+
+```rust
+use whittaker_eilers::WhittakerSmoother;
+
+let x_input = vec![1.1, 1.9, 3.1, 3.91, 5.0, 6.02, 7.01, 7.7, 9.0, 10.0 ,11.0, 12.0, 13.0];
+let data_to_smooth = vec![1.1, 1.9, 3.1, 3.91, 5.0, 6.02, 7.01, 7.7, 9.0, 10.0, 11.0, 12.0, 13.0];
 
----
 
+let y_inputs: Vec<Vec<f64>> = (0..10).map(|_| data_to_smooth.clone()).collect();
+
+let smoother = WhittakerSmoother::new(2e4, 2, y_inputs[0].len(), None, None).unwrap();
+
+let results: Vec<Vec<f64>> = smoother
+    .smooth_parallel(&y_inputs)
+    .into_iter()
+    .map(|x| x.unwrap())
+    .collect();
+
+println!("Result: {:?}", results);
+
+```
+
+---
 
 You can use these methods in combination with each other for instance, interpolating measurements without providing an x input. For more advanced examples of usage take a look at the examples, tests, and benches in the [Github](https://github.com/AnBowell/whittaker-eilers) repository. Here's an image of some smoothed data from an example:
 
 <img src="/examples/images/smoothed_data.png" alt="Time-series smoothed by Whittaker-Eilers method" width="800" />
 
-
 ## Further Reading
+
 If you'd like to see a more detailed run through of the library, check out this [Medium post](https://medium.com/towards-data-science/the-perfect-way-to-smooth-your-noisy-data-4f3fe6b44440). Within it, I run through examples and benchmarks against other smoothing methods.
 
 ## Future Features
+
 - Scatter plot smoothing
 - Generic typing
 
-
-
-
 ## References
+
 The algorithm implemented here mirrors a 2003 implementation by Paul H. C. Eilers in Matlab. I've included scripts and data from the original paper in the tests for this crate. The original paper and code can be found here:
 
 [A Perfect Smoother](https://pubs.acs.org/doi/10.1021/ac034173t)
diff --git a/benches/whittaker.rs b/benches/whittaker.rs
index 828c75c..16a9c3b 100644
--- a/benches/whittaker.rs
+++ b/benches/whittaker.rs
@@ -7,6 +7,17 @@ fn new_y_whittaker(y: &Vec<f64>) -> Vec<f64> {
         .smooth(y)
         .unwrap()
 }
+#[cfg(feature = "rayon")]
+fn new_parallel_y_whittaker(y: &Vec<Vec<f64>>) -> Vec<Vec<f64>> {
+    let smoother = WhittakerSmoother::new(2e4, 2, y[0].len(), None, None).unwrap();
+
+    smoother
+        .smooth_parallel(y)
+        .into_iter()
+        .map(|x| x.unwrap())
+        .collect()
+}
+
 fn new_y_whittaker_cross_validate(y: &Vec<f64>) -> CrossValidationResult {
     WhittakerSmoother::new(2e4, 2, y.len(), None, None)
         .unwrap()
@@ -35,9 +46,18 @@ fn new_x_y_weights_whittaker(x: &Vec<f64>, y: &Vec<f64>, weights: &Vec<f64>) ->
 
 fn criterion_benchmark(c: &mut Criterion) {
     let wood_data_vec: Vec<f64> = WOOD_DATASET.to_vec();
+
     let wood_x_vec: Vec<f64> = (0..wood_data_vec.len()).map(|x| x as f64).collect();
+
     let weights = vec![1.0; wood_data_vec.len()];
 
+    let repeat_wood_data = (0..100000)
+        .into_iter()
+        .map(|_| wood_data_vec.clone())
+        .collect();
+
+    // let repeat_wood_data = [wood_data_vec.as_slice(); 100000];
+
     let reusable_smoother = WhittakerSmoother::new(
         2e4,
         2,
@@ -46,9 +66,11 @@ fn criterion_benchmark(c: &mut Criterion) {
         Some(&weights),
     )
     .unwrap();
+
     c.bench_function("Whittaker Wood Y Only", |b| {
         b.iter(|| new_y_whittaker(black_box(&wood_data_vec)))
     });
+
     c.bench_function("Whittaker Wood Cross validate", |b| {
         b.iter(|| new_y_whittaker_cross_validate(black_box(&wood_data_vec)))
     });
@@ -67,9 +89,15 @@ fn criterion_benchmark(c: &mut Criterion) {
             )
         })
     });
+
     c.bench_function("Whittaker Wood X, Y, and weights reused", |b| {
         b.iter(|| reusable_smoother.smooth(&wood_data_vec).unwrap())
     });
+
+    #[cfg(feature = "rayon")]
+    c.bench_function("Whittaker Wood Y Only Parallel", |b| {
+        b.iter(|| new_parallel_y_whittaker(black_box(&repeat_wood_data)))
+    });
 }
 criterion_group!(benches, criterion_benchmark);
 criterion_main!(benches);
diff --git a/examples/smooth_and_interpolate.rs b/examples/smooth_and_interpolate.rs
index 17e6584..b5135c8 100644
--- a/examples/smooth_and_interpolate.rs
+++ b/examples/smooth_and_interpolate.rs
@@ -20,7 +20,7 @@ fn basic_smooth(x_input: &Vec<f64>, y_input: &Vec<f64>, lambda: f64, order: usiz
     plot.add_trace(raw_points);
     plot.add_trace(smoothed_points);
 
-    let layout = Layout::new().title("Basic smoothing of equally spaced Y".into());
+    let layout = Layout::new().title("Basic smoothing of equally spaced Y");
     plot.set_layout(layout);
 
     plot.show();
@@ -49,7 +49,7 @@ fn smooth_with_x(x_input_with_noise: &Vec<f64>, y_input: &Vec<f64>, lambda: f64,
     plot.add_trace(raw_points);
     plot.add_trace(smoothed_points);
 
-    let layout = Layout::new().title("Basic smoothing of Y with arbitrary X".into());
+    let layout = Layout::new().title("Basic smoothing of Y with arbitrary X");
     plot.set_layout(layout);
     plot.show();
 }
@@ -83,7 +83,7 @@ fn smooth_with_weights(
     plot.add_trace(raw_points);
     plot.add_trace(smoothed_points);
 
-    let layout = Layout::new().title("Basic smoothing of Y with arbitrary X and weights".into());
+    let layout = Layout::new().title("Basic smoothing of Y with arbitrary X and weights");
     plot.set_layout(layout);
     plot.show();
 }
@@ -112,7 +112,7 @@ fn smooth_and_interpolate(
     plot.add_trace(raw_points);
     plot.add_trace(smoothed_points);
 
-    let layout = Layout::new().title("Smoothing and Interpolation using weights".into());
+    let layout = Layout::new().title("Smoothing and Interpolation using weights");
     plot.set_layout(layout);
     plot.show();
 }
@@ -138,7 +138,7 @@ fn smooth_cross_validate(x_input: &Vec<f64>, y_input: &Vec<f64>, weights: &Vec<f
     plot.add_trace(raw_points);
     plot.add_trace(smoothed_points);
 
-    let layout = Layout::new().title("Whittaker optimal cross validation".into());
+    let layout = Layout::new().title("Whittaker optimal cross validation");
     plot.set_layout(layout);
     plot.show();
 
@@ -161,7 +161,7 @@ fn smooth_cross_validate(x_input: &Vec<f64>, y_input: &Vec<f64>, weights: &Vec<f
 
     plot.add_trace(smoothed_points);
 
-    let layout = Layout::new().title("Cross Validation Error".into());
+    let layout = Layout::new().title("Cross Validation Error");
     plot.set_layout(layout);
     plot.show();
 }
diff --git a/src/whittaker_smoother.rs b/src/whittaker_smoother.rs
index 0305ba5..de7be5a 100644
--- a/src/whittaker_smoother.rs
+++ b/src/whittaker_smoother.rs
@@ -8,6 +8,9 @@ use sprs::SymmetryCheck::CheckSymmetry;
 use sprs::{CsMat, CsMatView};
 use sprs_ldl::{Ldl, LdlNumeric};
 
+#[cfg(feature = "rayon")]
+use rayon::iter::{IntoParallelRefIterator, ParallelIterator};
+
 /// Whitaker-Eilers Smoother and Interpolator
 ///
 /// The smoother must be created via [WhittakerSmoother::new()] and once created, can be reused to smooth multiple sets of data as
@@ -245,6 +248,29 @@ impl WhittakerSmoother {
         };
     }
 
+    /// Run parallel Whittaker-Eilers smoothing and interpolation for multiple data series
+    ///
+    /// Convenience function to smooth many series in parallel. It should only be used when many series need to be smoothed with the same length, values of x,
+    /// weights, order, and lambda. If any parameters need to be changed, then the core Whittaker struct would need to be
+    /// mutated which makes parallel processing tough without cloning it repeatedly. At this point, it would be easier for the user to handle the parallelization!
+    ///
+    /// # Arguments
+    /// * `y_inputs`: A vec or slice of vecs or slices which are to be smoothed and interpolated by the Whittaker-Eilers smoother.
+    ///
+    /// # Returns:
+    /// The smoothed and interpolated data as a vector of results.
+    #[cfg(feature = "rayon")]
+    pub fn smooth_parallel<'a, T>(&self, y_inputs: &'a [T]) -> Vec<Result<Vec<f64>, WhittakerError>>
+    where
+        [T]: IntoParallelRefIterator<'a>,
+        <[T] as IntoParallelRefIterator<'a>>::Item: AsRef<[f64]>,
+    {
+        y_inputs
+            .par_iter()
+            .map(|y_input| self.smooth(y_input.as_ref()))
+            .collect()
+    }
+
     /// Run Whittaker-Eilers smoothing, interpolation and cross validation.
     ///
     /// This function will run the smoother and assess the cross validation error on the result. This is defined in Eilers'
diff --git a/whittaker-eilers-py/Cargo.toml b/whittaker-eilers-py/Cargo.toml
index ad34e6b..aac11aa 100644
--- a/whittaker-eilers-py/Cargo.toml
+++ b/whittaker-eilers-py/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "whittaker-eilers-py"
-version = "0.1.3"
+version = "0.2.0"
 edition = "2021"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
@@ -11,9 +11,11 @@ name = "whittaker_eilers"
 crate-type = ["cdylib"]
 
 [dependencies]
-whittaker-eilers-rs = { package = "whittaker-eilers", path = "../" }
+whittaker-eilers-rs = { package = "whittaker-eilers", path = "../", features = [
+    "rayon",
+] }
 
 [dependencies.pyo3]
-version = "0.19.2"
+version = "0.23.2"
 # "abi3-py37" tells pyo3 (and maturin) to build using the stable ABI with minimum Python version 3.7
 features = ["abi3-py37"]
diff --git a/whittaker-eilers-py/README.md b/whittaker-eilers-py/README.md
index fca6fcc..3ef047f 100644
--- a/whittaker-eilers-py/README.md
+++ b/whittaker-eilers-py/README.md
@@ -1,4 +1,5 @@
 # Whittaker-Eilers Smoothing and Interpolation
+
 **The Whittaker-Eilers smoother is the perfect smoother.** It offers extremely quick, efficient smoothing with built-in interpolation via weights on each measurement. This package provides a sparse-matrix implementation for additional speed and memory efficiency and can handle both equally and unequally spaced measurements. This package was originally written in Rust so additional examples, tests, and benchmarks are also available in addition to it being super speedy. The API is almost identical.
 
 ---
@@ -8,10 +9,13 @@ pip install whittaker-eilers
 ```
 
 ## Usage
+
 To start smoothing and interpolating data, create a reusable WhittakerSmoother class. You'll only need to recreate this class if the length or sampling rate of your data changes.
 
 ### Equally spaced data
+
 This is the fastest smoothing option. It smooths equally spaced y measurements using two tunable parameters, `lambda` (2e4) and the smoother `order` (2). The larger the lambda, the smoother the data.
+
 ```python
 from whittaker_eilers import WhittakerSmoother
 
@@ -24,10 +28,10 @@ smoothed_data = whittaker_smoother.smooth(data_to_smooth)
 print("Smoothed data: {}".format(smoothed_data))
 ```
 
+### Non-equally spaced data
 
+If you wish to smooth unequally spaced data, you need to provide an `x_input` with the sample times/positions.
 
-### Non-equally spaced data
-If you wish to smooth unequally spaced data, you need to provide an `x_input` with the sample times/positions. 
 ```python
 from whittaker_eilers import WhittakerSmoother
 
@@ -46,7 +50,9 @@ print("Smoothed non-equally spaced data: {}".format(smoothed_data))
 ```
 
 ### Weighted data & Interpolation
-Each measurement can then be weighted to trust some measurements more than others. Setting `weights` to 0 for measurements will lead to interpolation. 
+
+Each measurement can then be weighted to trust some measurements more than others. Setting `weights` to 0 for measurements will lead to interpolation.
+
 ```python
 from whittaker_eilers import WhittakerSmoother
 
@@ -70,8 +76,8 @@ print("Smoothed and interpolated weighted data: {}".format(smoothed_data))
 ```
 
 ### Smoothing with cross validation
-With this package, you can also calculate the cross validation error alongside the smoothed series. This shouldn't really be used in production where speed is necessary though!
 
+With this package, you can also calculate the cross validation error alongside the smoothed series. This shouldn't really be used in production where speed is necessary though!
 
 ```python
 from whittaker_eilers import WhittakerSmoother
@@ -95,11 +101,11 @@ print("Smoothed :{}".format(smoothed_data_with_cross_validation.get_smoothed()))
 ```
 
 ### Automatic smoothing
+
 Smoothing data requires a choice of Lambda. This can be done using visual inspection or by finding the lambda
 which results in the lowest cross validation error. The `smooth_optimal` function runs the smoother for a variety of lambdas and returns the results with the ability to retrieve the optimal one.
 
-
-```python 
+```python
 from whittaker_eilers import WhittakerSmoother
 
 x_input = [1.1, 1.9, 3.1, 3.91, 5.0, 6.02, 7.01, 7.7, 9.0, 10.0]
@@ -117,12 +123,29 @@ print("Optimally smoothed data: {}".format(optimal_smooth.get_optimal().get_smoo
 
 ```
 
+### Parallel Smoothing
+
+With this package you're able to smooth multiple data series in parallel. While y, the data to be smoothed can vary, x, weights, and the other smoothing parameters cannot be varied across the different data series.
+
+```python
+y = [[1.2, 4.5, 6.5, 3.2, 1.23, 23.23, 231.23]] * 10
+
+whittaker_smoother = WhittakerSmoother(lmbda=20, order=2, data_length=len(y[0]))
+
+result = whittaker_smoother.smooth_parallel(y)
+
+print("Smoothed data: {}".format(result))
+
+```
+
 ---
+
 You can use these methods in combination with each other for instance, interpolating measurements without providing an x input. For more advanced examples of usage take a look at the examples, tests, and benches in the [Github](https://github.com/AnBowell/whittaker-eilers) repository. Here's an image of some smoothed data from an example:
 
 <img src="/examples/images/smoothed_data.png" alt="Time-series smoothed by Whittaker-Eilers method" width="800" />
 
 ### Other methods
+
 ```python
 from whittaker_eilers import WhittakerSmoother
 x_input = [1.1, 1.9, 3.1, 3.91, 5.0, 6.02, 7.01, 7.7, 9.0, 10.0]
@@ -145,14 +168,18 @@ whittaker_smoother.update_weights([0.5] * len(x_input))
 whittaker_smoother.update_order(3)
 whittaker_smoother.update_lambda(4321.0)
 ```
+
 ## Further Reading
+
 If you'd like to see a more detailed run through of the library, check out this [Medium post](https://medium.com/towards-data-science/the-perfect-way-to-smooth-your-noisy-data-4f3fe6b44440). Within it, I run through examples and benchmarks against other smoothing methods.
 
 ## Future Features
+
 - Scatter plot smoothing
 - Generic typing
 
 ## References
+
 The algorithm implemented here mirrors a 2003 implementation by Paul H. C. Eilers in Matlab. I've included scripts and data from the original paper in the tests for this package. The original paper and code can be found here:
 
 [A Perfect Smoother](https://pubs.acs.org/doi/10.1021/ac034173t)
diff --git a/whittaker-eilers-py/examples/parallel_smoothing.py b/whittaker-eilers-py/examples/parallel_smoothing.py
new file mode 100644
index 0000000..56a80ff
--- /dev/null
+++ b/whittaker-eilers-py/examples/parallel_smoothing.py
@@ -0,0 +1,81 @@
+"""Demonstrates smoothing many series in parallel"""
+
+from whittaker_eilers import WhittakerSmoother
+from time import perf_counter
+import numpy as np
+
+
+NUMBER_OF_SERIES = 100000
+
+
+def main():
+
+    data = np.loadtxt("graph.txt", skiprows=5)
+    years = data[:, 0]
+    temp_anom = data[:, 1]
+
+    smooth_only_y(temp_anom)
+    smooth_y_with_x(years, temp_anom)
+    smooth_y_with_x_and_weights(years, temp_anom)
+
+
+def smooth_only_y(y):
+    whittaker_smoother = WhittakerSmoother(lmbda=20, order=2, data_length=len(y))
+
+    many_inputs = [y] * NUMBER_OF_SERIES
+
+    start_time = perf_counter()
+    result = whittaker_smoother.smooth_parallel(many_inputs)
+    end_time = perf_counter()
+
+    print(
+        "Took: {:.2f}s to process {} time-series".format(
+            end_time - start_time, len(result)
+        )
+    )
+
+
+def smooth_y_with_x(x, y):
+    whittaker_smoother = WhittakerSmoother(
+        lmbda=90000, order=3, data_length=len(y), x_input=x
+    )
+
+    many_inputs = [y] * NUMBER_OF_SERIES
+
+    many_inputs[-1] = many_inputs[-1][
+        100:
+    ].tolist()  # Creates an error as it's of a different length! Represented as None...
+
+    start_time = perf_counter()
+    result = whittaker_smoother.smooth_parallel(many_inputs)
+    end_time = perf_counter()
+    print(
+        "Took: {:.2f}s to process {} time-series with x inputs".format(
+            end_time - start_time, len(result)
+        )
+    )
+
+
+def smooth_y_with_x_and_weights(x, y):
+
+    weights = np.full(x.size, 1.0)
+    weights[::5] = 0.0
+
+    whittaker_smoother = WhittakerSmoother(
+        lmbda=20, order=2, data_length=len(y), x_input=x, weights=weights.tolist()
+    )
+    many_inputs = [y] * NUMBER_OF_SERIES
+
+    start_time = perf_counter()
+    result = whittaker_smoother.smooth_parallel(many_inputs)
+    end_time = perf_counter()
+
+    print(
+        "Took: {:.2f}s to process {} time-series with x inputs & weights".format(
+            end_time - start_time, len(result)
+        )
+    )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/whittaker-eilers-py/pyproject.toml b/whittaker-eilers-py/pyproject.toml
index 2a05a61..cf5a5c3 100644
--- a/whittaker-eilers-py/pyproject.toml
+++ b/whittaker-eilers-py/pyproject.toml
@@ -1,5 +1,5 @@
 [build-system]
-requires = ["maturin>=1.0,<2.0"]
+requires = ["maturin==1.6.0"]
 build-backend = "maturin"
 
 [tool.maturin]
diff --git a/whittaker-eilers-py/src/lib.rs b/whittaker-eilers-py/src/lib.rs
index b791862..b5c31cf 100644
--- a/whittaker-eilers-py/src/lib.rs
+++ b/whittaker-eilers-py/src/lib.rs
@@ -2,11 +2,15 @@ mod cross_validation;
 mod errors;
 mod whittaker_smoother;
 
-use pyo3::{pymodule, types::PyModule, PyResult, Python};
+use pyo3::{
+    pymodule,
+    types::{PyModule, PyModuleMethods},
+    Bound, PyResult,
+};
 use whittaker_smoother::WhittakerSmoother;
 
 #[pymodule]
-fn whittaker_eilers(_py: Python, m: &PyModule) -> PyResult<()> {
+fn whittaker_eilers(m: &Bound<'_, PyModule>) -> PyResult<()> {
     #![doc = include_str!("../README.md")]
     m.add_class::<WhittakerSmoother>()?;
     Ok(())
diff --git a/whittaker-eilers-py/src/whittaker_smoother.rs b/whittaker-eilers-py/src/whittaker_smoother.rs
index ae54f67..7744884 100644
--- a/whittaker-eilers-py/src/whittaker_smoother.rs
+++ b/whittaker-eilers-py/src/whittaker_smoother.rs
@@ -30,7 +30,7 @@ pub struct WhittakerSmoother(WhittakerSmootherRs);
 #[pymethods]
 impl WhittakerSmoother {
     #[new]
-    // #[pyo3(signature = (lmbda, order, data_length, x_input, weights), text_signature = "(lmbda, order, data_length, x_input, weights)")]
+    #[pyo3(signature = (lmbda, order, data_length, x_input=None, weights=None))]
     pub fn __init__(
         lmbda: f64, // Lambda is a key word in python
         order: usize,
@@ -111,6 +111,29 @@ impl WhittakerSmoother {
         self.0.smooth(&y_vals).map_err(map_err_to_py)
     }
 
+    /// Run parallel Whittaker-Eilers smoothing and interpolation for multiple data series
+    ///
+    /// Convenience function to smooth many series in parallel. It should only be used when many series need to be smoothed with the same length, values of x,
+    /// weights, order, and lambda. If any parameters need to be changed, then the core Whittaker struct/class would need to be
+    /// mutated which makes parallel processing tough without cloning it repeatedly. At this point, it would be easier for the user to handle the parallelization!
+    ///
+    /// Will return None for any time series which result in errors. Run erroneous inputs through the standard smoothing function to diagnose errors.
+    ///
+    /// Parameters
+    /// ----------
+    /// vals_y : The list of lists of values which are to be smoothed and interpolated by the Whittaker-Eilers smoother.
+    ///
+    /// Returns
+    /// -------
+    /// The smoothed and interpolated data for each list. None when the smoothing fails.
+    pub fn smooth_parallel(&self, y_val_series: Vec<Vec<f64>>) -> Vec<Option<Vec<f64>>> {
+        self.0
+            .smooth_parallel(&y_val_series)
+            .into_iter()
+            .map(|e| e.ok()) // Can't return a pyresult inside a vec so have opted for an option, otherwise you could smooth 999999 series and have an error crash it on the very last one.
+            .collect()
+    }
+
     /// Run Whittaker-Eilers smoothing, interpolation and cross validation.
     ///
     /// This function will run the smoother and assess the cross validation error on the result. This is defined in Eiler's
diff --git a/whittaker-eilers-py/whittaker_eilers.pyi b/whittaker-eilers-py/whittaker_eilers.pyi
index 0c21e07..242ead0 100644
--- a/whittaker-eilers-py/whittaker_eilers.pyi
+++ b/whittaker-eilers-py/whittaker_eilers.pyi
@@ -35,8 +35,8 @@ class WhittakerSmoother:
         lmbda: float,
         order: int,
         data_length: int,
-        x_input: Optional[List] = None,
-        weights: Optional[List] = None,
+        x_input: Optional[List[float]] = None,
+        weights: Optional[List[float]] = None,
     ) -> None: ...
     def get_order(self) -> int:
         """Retrieve the smoother's current order."""
@@ -44,9 +44,11 @@ class WhittakerSmoother:
     def get_lambda(self) -> float:
         """Retrieve the smoother's current lambda."""
         ...
+
     def get_data_length(self) -> int:
         """Retrieve the smoother's current length."""
         ...
+
     def update_weights(self, weights: List) -> None:
         """Updates the weights of the data to be smoothed.
         The length of weights should be equal to that of the data you are to smooth. The values of the weights should fall between 0 and 1.
@@ -56,6 +58,7 @@ class WhittakerSmoother:
         weights : The weights of the measurements to be smoothed. The smaller the weight the more the measurement will be ignored. Setting a weight to 0 results in interpolation.
         """
         ...
+
     def update_order(self, order: int) -> None:
         """Updates the order of the Whittaker-Eilers smoother.
 
@@ -65,6 +68,7 @@ class WhittakerSmoother:
         ----------
         order : The order to smooth."""
         ...
+
     def update_lambda(self, lmbda: float) -> None:
         """Updates the smoothing constant `lambda` of the Whittaker-Eilers smoother.
 
@@ -75,7 +79,8 @@ class WhittakerSmoother:
         lmbda : The smoothing constant of the Whittaker-Eilers smoother.
         """
         ...
-    def smooth(self, y_vals: List[float]) -> List:
+
+    def smooth(self, y_vals: List[float]) -> List[float]:
         """Run Whittaker-Eilers smoothing and interpolation.
 
         This function actually runs the solver which results in the smoothed data. If you just wish to continuously smooth
@@ -91,6 +96,27 @@ class WhittakerSmoother:
         The smoothed and interpolated data."""
         ...
 
+    def smooth_parallel(
+        self, y_val_series: List[List[float]]
+    ) -> List[List[float] | None]:
+        """Run parallel Whittaker-Eilers smoothing and interpolation for multiple data series
+
+        Convenience function to smooth many series in parallel. It should only be used when many series need to be smoothed with the same length, values of x,
+        weights, order, and lambda. If any parameters need to be changed, then the core Whittaker struct/class would need to be
+        mutated which makes parallel processing tough without cloning it repeatedly. At this point, it would be easier for the user to handle the parallelization!
+
+        Will return None for any time series which result in errors. Run erroneous inputs through the standard smoothing function to diagnose errors.
+
+        Parameters
+        ----------
+        vals_y : The list of lists of values which are to be smoothed and interpolated by the Whittaker-Eilers smoother.
+
+        Returns
+        -------
+        The smoothed and interpolated data for each list. None when the smoothing fails.
+        """
+        ...
+
     def smooth_and_cross_validate(self, y_vals: List[float]) -> CrossValidationResult:
         """Run Whittaker-Eilers smoothing, interpolation and cross validation.
 


Test patch:
diff --git a/tests/validation.rs b/tests/validation.rs
index 1161b53..d7385a4 100644
--- a/tests/validation.rs
+++ b/tests/validation.rs
@@ -42,6 +42,50 @@ fn validate_standard_whittaker() {
     }
 }
 
+#[test]
+#[cfg(feature = "rayon")]
+fn validate_standard_parallel_whittaker() {
+    let input_data = read_input_to_vecs();
+    let original_result_order_2 = read_output_to_vec("tests/data/output/output_only_y_2e4_2.csv");
+
+    let mut whittaker_smoother =
+        WhittakerSmoother::new(2e4, 2, input_data.y.len(), None, None).unwrap();
+
+    let multiple_series_input = [input_data.y.as_slice(); 5];
+
+    let rust_whittaker_out_vec_res = whittaker_smoother.smooth_parallel(&multiple_series_input);
+
+    let rust_whittaker_out_vec: Vec<&Vec<f64>> = rust_whittaker_out_vec_res
+        .iter()
+        .map(|x| x.as_ref().unwrap())
+        .collect();
+
+    for rust_whittaker_out in rust_whittaker_out_vec {
+        assert_eq!(original_result_order_2.len(), rust_whittaker_out.len());
+        for i in 0..original_result_order_2.len() {
+            assert_relative_eq!(
+                original_result_order_2[i],
+                rust_whittaker_out[i],
+                epsilon = 1e-8
+            )
+        }
+    }
+    let original_result_order_3 = read_output_to_vec("tests/data/output/output_only_y_2e4_3.csv");
+
+    whittaker_smoother.update_order(3).unwrap();
+
+    let rust_whittaker_out = whittaker_smoother.smooth(&input_data.y).unwrap();
+
+    assert_eq!(original_result_order_3.len(), rust_whittaker_out.len());
+    for i in 0..original_result_order_3.len() {
+        assert_relative_eq!(
+            original_result_order_3[i],
+            rust_whittaker_out[i],
+            epsilon = 1e-8
+        )
+    }
+}
+
 #[test]
 fn validated_weighted_whittaker() {
     let input_data = read_input_to_vecs();