5 files changed, 143 insertions, 4 deletions

diff --git a/CHANGES.md b/CHANGES.md
index 2788c52..0ce7699 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -1,3 +1,8 @@
+0.14.1.0
+--------
+
+- convolution
+
 0.14.0.0
 --------
 
diff --git a/hmatrix.cabal b/hmatrix.cabal
index 4d3beb6..0ef5cd7 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.14.0.1
+Version:            0.14.1.0
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -112,6 +112,7 @@ library
                         Numeric.LinearAlgebra.LAPACK,
                         Numeric.LinearAlgebra.Algorithms,
                         Numeric.LinearAlgebra.Util,
+                        Numeric.LinearAlgebra.Util.Convolution,
                         Graphics.Plot,
                         Data.Packed.ST,
                         Data.Packed.Development
diff --git a/lib/Numeric/LinearAlgebra/Util/Convolution.hs b/lib/Numeric/LinearAlgebra/Util/Convolution.hs
new file mode 100644
index 0000000..b64b169
--- /dev/null
+++ b/lib/Numeric/LinearAlgebra/Util/Convolution.hs
@@ -0,0 +1,118 @@
+{-# LANGUAGE FlexibleContexts #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra.Util.Convolution
+Copyright   :  (c) Alberto Ruiz 2012
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+
+-}
+-----------------------------------------------------------------------------
+
+module Numeric.LinearAlgebra.Util.Convolution(
+   -- * 1D
+   corr, conv,
+   -- * 2D
+   corr2, conv2,
+   -- * Misc
+   separable,  corrMin
+) where
+
+import Numeric.LinearAlgebra
+
+
+vectSS :: Element t => Int -> Vector t -> Matrix t
+vectSS n v = fromRows [ subVector k n v | k <- [0 .. dim v - n] ]
+
+
+corr :: Product t => Vector t -- ^ kernel
+                  -> Vector t -- ^ source
+                  -> Vector t
+{- ^ correlation
+
+@\> (fromList[1,2,3]) (fromList [1..10])
+fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0]@
+
+-}
+corr ker v | dim ker <= dim v = vectSS (dim ker) v <> ker
+           | otherwise = error $ "corr: dim kernel ("++show (dim ker)++") > dim vector ("++show (dim v)++")"
+
+
+conv :: (Product t, Num t) => Vector t -> Vector t -> Vector t
+{- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product)
+
+@\> conv (fromList[1,1]) (fromList [-1,1])
+fromList [-1.0,0.0,1.0]@
+
+-}
+conv ker v = corr ker' v'
+  where
+    ker' = (flatten.fliprl.asRow) ker
+    v' | dim ker > 1 = join [z,v,z]
+       | otherwise   = v
+    z = constant 0 (dim ker -1)
+
+corrMin :: (Container Vector t, RealElement t, Product t)
+        => Vector t
+        -> Vector t
+        -> Vector t
+-- ^ similar to 'corr', using 'min' instead of (*)
+corrMin ker v = minEvery ss (asRow ker) <> ones
+  where
+    minEvery a b = cond a b a a b
+    ss = vectSS (dim ker) v
+    ones = konst' 1 (dim ker)
+
+
+
+matSS :: Element t => Int -> Matrix t -> [Matrix t]
+matSS dr m = map (reshape c) [ subVector (k*c) n v | k <- [0 .. r - dr] ]
+  where
+    v = flatten m
+    c = cols m
+    r = rows m
+    n = dr*c
+
+
+corr2 :: Product a => Matrix a -> Matrix a -> Matrix a
+-- ^ 2D correlation
+corr2 ker mat = dims
+              . concatMap (map ((<.> ker') . flatten) . matSS c . trans)
+              . matSS r $ mat
+  where
+    r = rows ker
+    c = cols ker
+    ker' = flatten (trans ker)
+    rr = rows mat - r + 1
+    rc = cols mat - c + 1
+    dims | rr > 0 && rc > 0 = (rr >< rc)
+         | otherwise = error $ "corr2: dim kernel ("++sz ker++") > dim matrix ("++sz mat++")"
+    sz m = show (rows m)++"x"++show (cols m)
+
+conv2 :: (Num a, Product a) => Matrix a -> Matrix a -> Matrix a
+-- ^ 2D convolution
+conv2 k m = corr2 (fliprl . flipud $ k) pm
+  where
+    pm | r == 0 && c == 0 = m
+       | r == 0           = fromBlocks [[z3,m,z3]]
+       |           c == 0 = fromBlocks [[z2],[m],[z2]]
+       | otherwise        = fromBlocks [[z1,z2,z1]
+                                       ,[z3, m,z3]
+                                       ,[z1,z2,z1]]
+    r = rows k - 1
+    c = cols k - 1
+    h = rows m
+    w = cols m
+    z1 = konst' 0 (r,c)
+    z2 = konst' 0 (r,w)
+    z3 = konst' 0 (h,c)
+
+-- TODO: could be simplified using future empty arrays
+
+
+separable :: Element t => (Vector t -> Vector t) -> Matrix t -> Matrix t
+-- ^ 2D process implemented as separated 1D processes by rows and columns.
+separable f = fromColumns . map f . toColumns . fromRows . map f . toRows
+
diff --git a/packages/tests/hmatrix-tests.cabal b/packages/tests/hmatrix-tests.cabal
index 1c6443e..10b7c80 100644
--- a/packages/tests/hmatrix-tests.cabal
+++ b/packages/tests/hmatrix-tests.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix-tests
-Version:            0.2
+Version:            0.3
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -21,7 +21,7 @@ extra-source-files: CHANGES
 library
 
     Build-Depends:      base >= 4 && < 5,
-                        hmatrix >= 0.14,
+                        hmatrix >= 0.14.1,
                         QuickCheck >= 2, HUnit, random
 
     hs-source-dirs:     src
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index fd66767..4e9a521 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -42,6 +42,7 @@ import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr)
 import Control.Arrow((***))
 import Debug.Trace
 import Control.Monad(when)
+import Numeric.LinearAlgebra.Util.Convolution
 
 import Data.Packed.ST
 
@@ -421,7 +422,20 @@ accumTest = utest "accum" ok
          &&
          toList (flatten x) == [1,0,0,0,1,0,0,0,1] 
 
----------------------------------------------------------------------
+--------------------------------------------------------------------------------
+
+convolutionTest = utest "convolution" ok
+  where
+--    a = fromList [1..10]               :: Vector Double
+    b = fromList [1..3]                :: Vector Double
+    c = (5><7) [1..]                   :: Matrix Double
+--    d = (3><3) [0,-1,0,-1,4,-1,0,-1,0] :: Matrix Double
+    ok =  separable (corr b) c == corr2 (outer b b) c
+       && separable (conv b) c == conv2 (outer b b) c
+
+--------------------------------------------------------------------------------
+
+
 
 -- | All tests must pass with a maximum dimension of about 20
 --  (some tests may fail with bigger sizes due to precision loss).
@@ -596,6 +610,7 @@ runTests n = do
         , condTest
         , conformTest
         , accumTest
+        , convolutionTest
         ]
     when (errors c + failures c > 0) exitFailure
     return ()