hmatrix 0.14.0.1 → 0.14.1.0

raw patch · 4 files changed

+180/−6 lines, 4 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• CHANGES.md +7/−0
• hmatrix.cabal +3/−2
• lib/Numeric/LinearAlgebra/Util.hs +56/−4
• lib/Numeric/LinearAlgebra/Util/Convolution.hs +114/−0

@@ -1,3 +1,10 @@+0.14.1.0+--------++- In Numeric.LinearAlgebra.Util:+    convolution: corr, conv, corr2, conv2, separable, corrMin+    kronecker: vec, vech, dup, vtrans+ 0.14.0.0 -------- 

@@ -1,5 +1,5 @@ Name:               hmatrix-Version:            0.14.0.1+Version:            0.14.1.0 License:            GPL License-file:       LICENSE Author:             Alberto Ruiz@@ -128,7 +128,8 @@                         Numeric.IO,                         Numeric.Chain,                         Numeric.Vector,-                        Numeric.Matrix+                        Numeric.Matrix,+                        Numeric.LinearAlgebra.Util.Convolution      C-sources:          lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,                         lib/Numeric/GSL/gsl-aux.c

@@ -1,3 +1,4 @@+{-# LANGUAGE FlexibleContexts #-} ----------------------------------------------------------------------------- {- | Module      :  Numeric.LinearAlgebra.Util@@ -11,6 +12,7 @@ -----------------------------------------------------------------------------  module Numeric.LinearAlgebra.Util(+    -- * Convenience functions for real elements     disp,     zeros, ones,     diagl,@@ -19,11 +21,28 @@     (&),(!), (#),     rand, randn,     cross,-    norm+    norm,+    -- * Convolution+    -- ** 1D+    corr, conv, corrMin,+    -- ** 2D+    corr2, conv2, separable,+    -- * Tools for the Kronecker product+    --+    -- | (see A. Fusiello, A matter of notation: Several uses of the Kronecker product in+    --  3d computer vision, Pattern Recognition Letters 28 (15) (2007) 2127-2132)++    --+    -- | @`vec` (a \<> x \<> b) == ('trans' b ` 'kronecker' ` a) \<> 'vec' x@+    vec,+    vech,+    dup,+    vtrans ) where -import Numeric.LinearAlgebra+import Numeric.LinearAlgebra hiding (i) import System.Random(randomIO)+import Numeric.LinearAlgebra.Util.Convolution   disp :: Int -> Matrix Double -> IO ()@@ -87,7 +106,7 @@ col = asColumn . fromList  cross :: Vector Double -> Vector Double -> Vector Double--- ^ cross product of dimension 3 real vectors+-- ^ cross product (for three-element real vectors) cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3]           | otherwise = error $ "cross ("++show x++") ("++show y++")"   where@@ -98,7 +117,40 @@     z3 = x1*y2-x2*y1  norm :: Vector Double -> Double--- ^ 2-norm of real vectors+-- ^ 2-norm of real vector norm = pnorm PNorm2 +--------------------------------------------------------------------------------++vec :: Element t => Matrix t -> Vector t+-- ^ stacking of columns+vec = flatten . trans+++vech :: Element t => Matrix t -> Vector t+-- ^ half-vectorization (of the lower triangular part)+vech m = join . zipWith f [0..] . toColumns $ m+  where+    f k v = subVector k (dim v - k) v+++dup :: (Num t, Num (Vector t), Element t) => Int -> Matrix t+-- ^ duplication matrix (@'dup' k \<> 'vech' m == 'vec' m@, for symmetric m of 'dim' k)+dup k = trans $ fromRows $ map f es+  where+    rs = zip [0..] (toRows (ident (k^(2::Int))))+    es = [(i,j) | j <- [0..k-1], i <- [0..k-1], i>=j ]+    f (i,j) | i == j = g (k*j + i)+            | otherwise = g (k*j + i) + g (k*i + j)+    g j = v+      where+        Just v = lookup j rs+++vtrans :: Element t => Int -> Matrix t -> Matrix t+-- ^ generalized \"vector\" transposition: @'vtrans' 1 == 'trans'@, and @'vtrans' ('rows' m) m == 'asColumn' ('vec' m)@+vtrans p m | r == 0 = fromBlocks . map (map asColumn . takesV (replicate q p)) . toColumns $ m+           | otherwise = error $ "vtrans " ++ show p ++ " of matrix with " ++ show (rows m) ++ " rows"+  where+    (q,r) = divMod (rows m) p 

@@ -0,0 +1,114 @@+{-# 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(+   corr, conv, corrMin,+   corr2, conv2, separable+) 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++>>> corr (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+-- ^ matrix computation implemented as separated vector operations by rows and columns.+separable f = fromColumns . map f . toColumns . fromRows . map f . toRows+