diff --git a/src/linalg/convolution.rs b/src/linalg/convolution.rs new file mode 100644 index 00000000..b121b34a --- /dev/null +++ b/src/linalg/convolution.rs @@ -0,0 +1,129 @@ +use base::allocator::Allocator; +use base::default_allocator::DefaultAllocator; +use base::dimension::{Dim, DimAdd, DimDiff, DimMax, DimMaximum, DimSub, DimSum}; +use std::cmp; +use storage::Storage; +use {zero, Real, Vector, VectorN, U1}; + +impl> Vector { + /// Returns the convolution of the target vector and a kernel + /// + /// # Arguments + /// + /// * `kernel` - A Vector with size > 0 + /// + /// # Errors + /// Inputs must statisfy `vector.len() >= kernel.len() > 0`. + /// + pub fn convolve_full( + &self, + kernel: Vector, + ) -> VectorN, U1>> + where + D1: DimAdd, + D2: DimAdd>, + DimSum: DimSub, + S2: Storage, + DefaultAllocator: Allocator, U1>>, + { + let vec = self.len(); + let ker = kernel.len(); + + if ker == 0 || ker > vec { + panic!("convolve_full expects `self.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); + } + + let result_len = self.data.shape().0.add(kernel.data.shape().0).sub(U1); + let mut conv = VectorN::zeros_generic(result_len, U1); + + for i in 0..(vec + ker - 1) { + let u_i = if i > vec { i - ker } else { 0 }; + let u_f = cmp::min(i, vec - 1); + + if u_i == u_f { + conv[i] += self[u_i] * kernel[(i - u_i)]; + } else { + for u in u_i..(u_f + 1) { + if i - u < ker { + conv[i] += self[u] * kernel[(i - u)]; + } + } + } + } + conv + } + /// Returns the convolution of the target vector and a kernel + /// The output convolution consists only of those elements that do not rely on the zero-padding. + /// # Arguments + /// + /// * `kernel` - A Vector with size > 0 + /// + /// + /// # Errors + /// Inputs must statisfy `self.len() >= kernel.len() > 0`. + /// + pub fn convolve_valid(&self, kernel: Vector, + ) -> VectorN, D2>> + where + D1: DimAdd, + D2: Dim, + DimSum: DimSub, + S2: Storage, + DefaultAllocator: Allocator, D2>>, + { + let vec = self.len(); + let ker = kernel.len(); + + if ker == 0 || ker > vec { + panic!("convolve_valid expects `self.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); + } + + let result_len = self.data.shape().0.add(U1).sub(kernel.data.shape().0); + let mut conv = VectorN::zeros_generic(result_len, U1); + + for i in 0..(vec - ker + 1) { + for j in 0..ker { + conv[i] += self[i + j] * kernel[ker - j - 1]; + } + } + conv + } + + /// Returns the convolution of the targetvector and a kernel + /// The output convolution is the same size as vector, centered with respect to the ‘full’ output. + /// # Arguments + /// + /// * `kernel` - A Vector with size > 0 + /// + /// # Errors + /// Inputs must statisfy `self.len() >= kernel.len() > 0`. + pub fn convolve_same(&self, kernel: Vector) -> VectorN> + where + D1: DimMax, + D2: DimMax>, + S2: Storage, + DefaultAllocator: Allocator>, + { + let vec = self.len(); + let ker = kernel.len(); + + if ker == 0 || ker > vec { + panic!("convolve_same expects `self.len() >= kernel.len() > 0`, received {} and {} respectively.",vec,ker); + } + + let result_len = self.data.shape().0.max(kernel.data.shape().0); + let mut conv = VectorN::zeros_generic(result_len, U1); + + for i in 0..vec { + for j in 0..ker { + let val = if i + j < 1 || i + j >= vec + 1 { + zero::() + } else { + self[i + j - 1] + }; + conv[i] += val * kernel[ker - j - 1]; + } + } + conv + } +} diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index 4418b283..b6a9e8d8 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -17,6 +17,7 @@ mod solve; mod svd; mod symmetric_eigen; mod symmetric_tridiagonal; +mod convolution; //// FIXME: Not complete enough for publishing. //// This handles only cases where each eigenvalue has multiplicity one. @@ -33,3 +34,4 @@ pub use self::schur::*; pub use self::svd::*; pub use self::symmetric_eigen::*; pub use self::symmetric_tridiagonal::*; +pub use self::convolution::*; diff --git a/tests/linalg/convolution.rs b/tests/linalg/convolution.rs new file mode 100644 index 00000000..b0d57f72 --- /dev/null +++ b/tests/linalg/convolution.rs @@ -0,0 +1,119 @@ +use na::{Vector2,Vector3,Vector4,Vector5,DVector}; +use std::panic; + +// +// Should mimic calculations in Python's scipy library +// >>>from scipy.signal import convolve +// + +// >>> convolve([1,2,3,4],[1,2],"same") +// array([ 1, 4, 7, 10]) +#[test] +fn convolve_same_check(){ + // Static Tests + let actual_s = Vector4::from_vec(vec![1.0,4.0,7.0,10.0]); + let expected_s = Vector4::new(1.0,2.0,3.0,4.0).convolve_same(Vector2::new(1.0,2.0)); + + assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7)); + + // Dynamic Tests + let actual_d = DVector::from_vec(vec![1.0,4.0,7.0,10.0]); + let expected_d = DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_same(DVector::from_vec(vec![1.0,2.0])); + + assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); + + // Panic Tests + // These really only apply to dynamic sized vectors + assert!( + panic::catch_unwind(|| { + DVector::from_vec(vec![1.0,2.0]).convolve_same(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); + + assert!( + panic::catch_unwind(|| { + DVector::::from_vec(vec![]).convolve_same(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); + + assert!( + panic::catch_unwind(|| { + DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_same(DVector::::from_vec(vec![])); + }).is_err() + ); +} + +// >>> convolve([1,2,3,4],[1,2],"full") +// array([ 1, 4, 7, 10, 8]) +#[test] +fn convolve_full_check(){ + // Static Tests + let actual_s = Vector5::new(1.0,4.0,7.0,10.0,8.0); + let expected_s = Vector4::new(1.0,2.0,3.0,4.0).convolve_full(Vector2::new(1.0,2.0)); + + assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7)); + + // Dynamic Tests + let actual_d = DVector::from_vec(vec![1.0,4.0,7.0,10.0,8.0]); + let expected_d = DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_full(DVector::from_vec(vec![1.0,2.0])); + + assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); + + // Panic Tests + // These really only apply to dynamic sized vectors + assert!( + panic::catch_unwind(|| { + DVector::from_vec(vec![1.0,2.0]).convolve_full(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); + + assert!( + panic::catch_unwind(|| { + DVector::::from_vec(vec![]).convolve_full(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); + + assert!( + panic::catch_unwind(|| { + DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_full(DVector::::from_vec(vec![])); + }).is_err() + ); +} + +// >>> convolve([1,2,3,4],[1,2],"valid") +// array([ 4, 7, 10]) +#[test] +fn convolve_valid_check(){ + // Static Tests + let actual_s = Vector3::from_vec(vec![4.0,7.0,10.0]); + let expected_s = Vector4::new(1.0,2.0,3.0,4.0).convolve_valid( Vector2::new(1.0,2.0)); + + assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7)); + + // Dynamic Tests + let actual_d = DVector::from_vec(vec![4.0,7.0,10.0]); + let expected_d = DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_valid(DVector::from_vec(vec![1.0,2.0])); + + assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); + + // Panic Tests + // These really only apply to dynamic sized vectors + assert!( + panic::catch_unwind(|| { + DVector::from_vec(vec![1.0,2.0]).convolve_valid(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); + + assert!( + panic::catch_unwind(|| { + DVector::::from_vec(vec![]).convolve_valid(DVector::from_vec(vec![1.0,2.0,3.0,4.0])); + }).is_err() + ); + + assert!( + panic::catch_unwind(|| { + DVector::from_vec(vec![1.0,2.0,3.0,4.0]).convolve_valid(DVector::::from_vec(vec![])); + }).is_err() + ); + +} \ No newline at end of file diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs index 74a5e03c..4e0bf2eb 100644 --- a/tests/linalg/mod.rs +++ b/tests/linalg/mod.rs @@ -11,3 +11,4 @@ mod real_schur; mod solve; mod svd; mod tridiagonal; +mod convolution; \ No newline at end of file