diff --git a/sci-rs-core/Cargo.toml b/sci-rs-core/Cargo.toml new file mode 100644 index 00000000..26fe7017 --- /dev/null +++ b/sci-rs-core/Cargo.toml @@ -0,0 +1,30 @@ +[package] +name = "sci-rs-core" +version = "0.0.0" +edition = "2021" +authors = ["Jacob Trueb "] +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 } diff --git a/sci-rs-core/src/lib.rs b/sci-rs-core/src/lib.rs new file mode 100644 index 00000000..f52e09a4 --- /dev/null +++ b/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 = core::result::Result; + +/// 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; diff --git a/sci-rs-core/src/num_rs/convolve/mod.rs b/sci-rs-core/src/num_rs/convolve/mod.rs new file mode 100644 index 00000000..2d80c7f7 --- /dev/null +++ b/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(a: ArrayView1, v: ArrayView1, mode: ConvolveMode) -> Result> +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); + } +} diff --git a/sci-rs-core/src/num_rs/convolve/ndarray_conv_binds.rs b/sci-rs-core/src/num_rs/convolve/ndarray_conv_binds.rs new file mode 100644 index 00000000..b0f38c79 --- /dev/null +++ b/sci-rs-core/src/num_rs/convolve/ndarray_conv_binds.rs @@ -0,0 +1,12 @@ +use super::ConvolveMode; +use ndarray_conv::ConvMode; + +impl From for ConvMode { + fn from(value: ConvolveMode) -> Self { + match value { + ConvolveMode::Full => ConvMode::Full, + ConvolveMode::Same => ConvMode::Same, + ConvolveMode::Valid => ConvMode::Valid, + } + } +} diff --git a/sci-rs-core/src/num_rs/mod.rs b/sci-rs-core/src/num_rs/mod.rs new file mode 100644 index 00000000..e6e9679c --- /dev/null +++ b/sci-rs-core/src/num_rs/mod.rs @@ -0,0 +1,4 @@ +#[cfg(feature = "alloc")] +mod convolve; +#[cfg(feature = "alloc")] +pub use convolve::*; diff --git a/sci-rs/Cargo.toml b/sci-rs/Cargo.toml index 697e7933..99647819 100644 --- a/sci-rs/Cargo.toml +++ b/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" diff --git a/sci-rs/src/signal/convolve.rs b/sci-rs/src/signal/convolve.rs index b527805e..7e0486b5 100644 --- a/sci-rs/src/signal/convolve.rs +++ b/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. ///