numpy-esque linear convolve

sci-rs-core/Cargo.toml

@@ -0,0 +1,30 @@
+[package]
+name = "sci-rs-core"
+version = "0.0.0"
+edition = "2021"
+authors = ["Jacob Trueb <[email protected]>"]
+description = "Core library for sci-rs internals."
+license = "MIT OR Apache-2.0"
+repository = "https://github.com/qsib-cbie/sci-rs.git"
+homepage = "https://github.com/qsib-cbie/sci-rs.git"
+readme = "../README.md"
+keywords = ["scipy", "dsp", "signal", "filter", "design"]
+categories = ["science", "mathematics", "no-std", "embedded"]
+
+
+[package.metadata.docs.rs]
+all-features = true
+
+[features]
+default = ['alloc']
+
+# Allow allocating vecs, matrices, etc.
+alloc = []
+
+# Enable FFT and standard library features
+std = ['alloc']
+
+[dependencies]
+ndarray = { version = "0.16.1", default-features = false }
+ndarray-conv = { version = "0.5.0" }
+num-traits = { version = "0.2.15", default-features = false }

sci-rs-core/src/lib.rs

@@ -0,0 +1,79 @@
+//! Core library for sci-rs.
+
+#![cfg_attr(not(feature = "std"), no_std)]
+
+#[cfg(feature = "alloc")]
+extern crate alloc;
+#[cfg(feature = "alloc")]
+use alloc::format;
+
+use core::{error, fmt};
+
+pub type Result<T> = core::result::Result<T, Error>;
+
+/// Errors raised whilst running sci-rs.
+#[derive(Debug, PartialEq, Eq)]
+pub enum Error {
+    /// Argument parsed into function were invalid.
+    #[cfg(feature = "alloc")]
+    InvalidArg {
+        /// The invalid arg
+        arg: alloc::string::String,
+        /// Explaining why arg is invalid.
+        reason: alloc::string::String,
+    },
+    /// Argument parsed into function were invalid.
+    #[cfg(not(feature = "alloc"))]
+    InvalidArg,
+    /// Two or more optional arguments passed into functions conflict.
+    #[cfg(feature = "alloc")]
+    ConflictArg {
+        /// Explaining what arg is invalid.
+        reason: alloc::string::String,
+    },
+    /// Two or more optional arguments passed into functions conflict.
+    #[cfg(not(feature = "alloc"))]
+    ConflictArg,
+    /// Errors raised by [ndarray_conv::Error]
+    #[cfg(feature = "alloc")]
+    Conv { reason: alloc::string::String },
+    /// Errors raised by [ndarray_conv::Error]
+    #[cfg(not(feature = "alloc"))]
+    Conv,
+}
+
+impl fmt::Display for Error {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(
+            f,
+            "{}",
+            match self {
+                #[cfg(feature = "alloc")]
+                Error::InvalidArg { arg, reason } =>
+                    format!("Invalid Argument on arg = {} with reason = {}", arg, reason),
+                #[cfg(not(feature = "alloc"))]
+                Error::InvalidArg =>
+                    "There were invalid arguments. Reasons not shown without `alloc` feature.",
+                #[cfg(feature = "alloc")]
+                Error::ConflictArg { reason } =>
+                    format!("Conflicting Arguments with reason = {}", reason),
+                #[cfg(not(feature = "alloc"))]
+                Error::ConflictArg =>
+                    "There were conflicting arguments. Reasons not shown without `alloc` feature.",
+                #[cfg(feature = "alloc")]
+                Error::Conv { reason } => format!(
+                    "An error occurred during the convolution from ndarray_conv with reason {}.",
+                    reason
+                ),
+                #[cfg(not(feature = "alloc"))]
+                Error::Conv => "An error occurred during the convolution from ndarray_conv. Reasons not shown without `alloc` feature.",
+            }
+        )
+    }
+}
+
+impl error::Error for Error {}
+
+/// Collection of numpy-like functions for use by sci-rs.
+/// Provide behaviour parity against Numpy, even if the types are not identical.
+pub mod num_rs;

sci-rs-core/src/num_rs/convolve/mod.rs

@@ -0,0 +1,135 @@
+mod ndarray_conv_binds;
+
+use crate::{Error, Result};
+use alloc::string::ToString;
+use ndarray::{Array1, ArrayView1};
+use ndarray_conv::{ConvExt, PaddingMode};
+
+/// Convolution mode determines behavior near edges and output size
+pub enum ConvolveMode {
+    /// Full convolution, output size is `in1.len() + in2.len() - 1`
+    Full,
+    /// Valid convolution, output size is `max(in1.len(), in2.len()) - min(in1.len(), in2.len()) + 1`
+    Valid,
+    /// Same convolution, output size is `in1.len()`
+    Same,
+}
+
+/// Best effort parallel behaviour with numpy's convolve method. We take `v` as the convolution
+/// kernel.
+///
+/// Returns the discrete, linear convolution of two one-dimensional sequences.
+///
+/// # Parameters
+/// * `a` : (N,) [[array_like]]([ndarray::Array1])  
+///   Signal to be (linearly) convolved.
+/// * `v` : (M,) [[array_like]]([ndarray::Array1])  
+///   Second one-dimensional input array.
+/// * `mode` : [ConvolveMode]  
+///   [ConvolveMode::Full]:  
+///   By default, mode is 'full'.  This returns the convolution at each point of overlap, with an
+///   output shape of (N+M-1,). At the end-points of the convolution, the signals do not overlap
+///   completely, and boundary effects may be seen.
+///
+///   [ConvolveMode::Same]:  
+///   Mode 'same' returns output of length ``max(M, N)``.  Boundary effects are still visible.
+///
+///   [ConvolveMode::Valid]:  
+///   Mode 'valid' returns output of length ``max(M, N) - min(M, N) + 1``.  The convolution
+///   product is only given for points where the signals overlap completely.  Values outside the
+///   signal boundary have no effect.
+///
+/// # Panics
+/// We assume that `v` is shorter than `a`.
+///
+/// # Examples
+/// With [ConvolveMode::Full]:
+/// ```
+/// use ndarray::array;
+/// use sci_rs_core::num_rs::{ConvolveMode, convolve};
+///
+/// let a = array![1., 2., 3.];
+/// let v = array![0., 1., 0.5];
+///
+/// let expected = array![0., 1., 2.5, 4., 1.5];
+/// let result = convolve((&a).into(), (&v).into(), ConvolveMode::Full).unwrap();
+/// assert_eq!(result, expected);
+/// ```
+/// With [ConvolveMode::Same]:
+/// ```
+/// use ndarray::array;
+/// use sci_rs_core::num_rs::{ConvolveMode, convolve};
+///
+/// let a = array![1., 2., 3.];
+/// let v = array![0., 1., 0.5];
+///
+/// let expected = array![1., 2.5, 4.];
+/// let result = convolve((&a).into(), (&v).into(), ConvolveMode::Same).unwrap();
+/// assert_eq!(result, expected);
+/// ```
+/// With [ConvolveMode::Same]:
+/// ```
+/// use ndarray::array;
+/// use sci_rs_core::num_rs::{ConvolveMode, convolve};
+///
+/// let a = array![1., 2., 3.];
+/// let v = array![0., 1., 0.5];
+///
+/// let expected = array![2.5];
+/// let result = convolve((&a).into(), (&v).into(), ConvolveMode::Valid).unwrap();
+/// assert_eq!(result, expected);
+/// ```
+pub fn convolve<T>(a: ArrayView1<T>, v: ArrayView1<T>, mode: ConvolveMode) -> Result<Array1<T>>
+where
+    T: num_traits::NumAssign + core::marker::Copy,
+{
+    // Convolve
+    let result = a.conv(&v, mode.into(), PaddingMode::Zeros);
+    #[cfg(feature = "alloc")]
+    {
+        result.map_err(|e| Error::Conv {
+            reason: e.to_string(),
+        })
+    }
+    #[cfg(not(feature = "alloc"))]
+    {
+        result.map_err({ Error::Conv })
+    }
+}
+
+#[cfg(test)]
+mod linear_convolve {
+    use super::*;
+    use alloc::vec;
+    use ndarray::array;
+
+    #[test]
+    fn full() {
+        let a = array![1., 2., 3.];
+        let v = array![0., 1., 0.5];
+
+        let expected = array![0., 1., 2.5, 4., 1.5];
+        let result = convolve((&a).into(), (&v).into(), ConvolveMode::Full).unwrap();
+        assert_eq!(result, expected);
+    }
+
+    #[test]
+    fn same() {
+        let a = array![1., 2., 3.];
+        let v = array![0., 1., 0.5];
+
+        let expected = array![1., 2.5, 4.];
+        let result = convolve((&a).into(), (&v).into(), ConvolveMode::Same).unwrap();
+        assert_eq!(result, expected);
+    }
+
+    #[test]
+    fn valid() {
+        let a = array![1., 2., 3.];
+        let v = array![0., 1., 0.5];
+
+        let expected = array![2.5];
+        let result = convolve((&a).into(), (&v).into(), ConvolveMode::Valid).unwrap();
+        assert_eq!(result, expected);
+    }
+}

sci-rs-core/src/num_rs/convolve/ndarray_conv_binds.rs

@@ -0,0 +1,12 @@
+use super::ConvolveMode;
+use ndarray_conv::ConvMode;
+
+impl<const N: usize> From<ConvolveMode> for ConvMode<N> {
+    fn from(value: ConvolveMode) -> Self {
+        match value {
+            ConvolveMode::Full => ConvMode::Full,
+            ConvolveMode::Same => ConvMode::Same,
+            ConvolveMode::Valid => ConvMode::Valid,
+        }
+    }
+}

sci-rs-core/src/num_rs/mod.rs

@@ -0,0 +1,4 @@
+#[cfg(feature = "alloc")]
+mod convolve;
+#[cfg(feature = "alloc")]
+pub use convolve::*;

sci-rs/Cargo.toml

@@ -19,10 +19,10 @@ all-features = true
 default = ['alloc']
 
 # Allow allocating vecs, matrices, etc.
-alloc = ['nalgebra/alloc', 'nalgebra/libm', 'kalmanfilt/alloc']
+alloc = ['nalgebra/alloc', 'nalgebra/libm', 'kalmanfilt/alloc', 'sci-rs-core/alloc']
 
 # Enable FFT and standard library features
-std = ['nalgebra/std', 'nalgebra/macros', 'rustfft', 'alloc']
+std = ['nalgebra/std', 'nalgebra/macros', 'rustfft', 'alloc','sci-rs-core/std']
 
 # Enable debug plotting through python system calls
 plot = ['std']
@@ -36,6 +36,7 @@ lstsq = { version = "0.6.0", default-features = false }
 rustfft = { version = "6.2.0", optional = true }
 kalmanfilt = { version = "0.3.0", default-features = false }
 gaussfilt = { version = "0.1.3", default-features = false }
+sci-rs-core = { path = "../sci-rs-core", default-features = false }
 
 [dev-dependencies]
 approx = "0.5.1"

sci-rs/src/signal/convolve.rs

@@ -2,15 +2,7 @@ use nalgebra::Complex;
 use num_traits::{Float, FromPrimitive, Signed, Zero};
 use rustfft::{FftNum, FftPlanner};
 
-/// Convolution mode determines behavior near edges and output size
-pub enum ConvolveMode {
-    /// Full convolution, output size is `in1.len() + in2.len() - 1`
-    Full,
-    /// Valid convolution, output size is `max(in1.len(), in2.len()) - min(in1.len(), in2.len()) + 1`
-    Valid,
-    /// Same convolution, output size is `in1.len()`
-    Same,
-}
+pub use sci_rs_core::num_rs::ConvolveMode;
 
 /// Performs FFT-based convolution on two slices of floating point values.
 ///