Merge pull request #539 from npapapietro/Implement_convolution_#520

Implement convolution #520

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<N: Real, D1: Dim, S1: Storage<N, D1>> Vector<N, D1, S1> {
+    /// 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<D2, S2>(
+        &self,
+        kernel: Vector<N, D2, S2>,
+    ) -> VectorN<N, DimDiff<DimSum<D1, D2>, U1>>
+    where
+        D1: DimAdd<D2>,
+        D2: DimAdd<D1, Output = DimSum<D1, D2>>,
+        DimSum<D1, D2>: DimSub<U1>,
+        S2: Storage<N, D2>,
+        DefaultAllocator: Allocator<N, DimDiff<DimSum<D1, D2>, 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<D2, S2>(&self, kernel: Vector<N, D2, S2>,
+    ) -> VectorN<N, DimDiff<DimSum<D1, U1>, D2>>
+    where
+        D1: DimAdd<U1>,
+        D2: Dim,
+        DimSum<D1, U1>: DimSub<D2>,
+        S2: Storage<N, D2>,
+        DefaultAllocator: Allocator<N, DimDiff<DimSum<D1, U1>, 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<D2, S2>(&self, kernel: Vector<N, D2, S2>) -> VectorN<N, DimMaximum<D1, D2>>
+    where
+        D1: DimMax<D2>,
+        D2: DimMax<D1, Output = DimMaximum<D1, D2>>,
+        S2: Storage<N, D2>,
+        DefaultAllocator: Allocator<N, DimMaximum<D1, D2>>,
+    {
+        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::<N>()
+                } else {
+                    self[i + j - 1]
+                };
+                conv[i] += val * kernel[ker - j - 1];
+            }
+        }
+        conv
+    }
+}

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::*;

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::<f32>::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::<f32>::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::<f32>::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::<f32>::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::<f32>::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::<f32>::from_vec(vec![]));
+        }).is_err()
+    );
+
+}

tests/linalg/mod.rs

@@ -11,3 +11,4 @@ mod real_schur;
 mod solve;
 mod svd;
 mod tridiagonal;
+mod convolution;