linear algebra moved

6 files changed, 588 insertions, 4 deletions

diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal
index 9a84b7e..3f12dad 100644
--- a/packages/base/hmatrix-base.cabal
+++ b/packages/base/hmatrix-base.cabal
@@ -38,11 +38,11 @@ library
                         Numeric.Conversion
                         Numeric.LinearAlgebra.LAPACK
                         Numeric.LinearAlgebra.Algorithms
-                        Data.Packed.Numeric
-                        Numeric.Vectorized
-                        Numeric.Vector
-                        Numeric.Matrix
                         Numeric.Container
+                        Numeric.LinearAlgebra
+                        Numeric.LinearAlgebra.Util
+                        Numeric.LinearAlgebra.Util.Convolution
+                        Data.Packed.IO
                         
     other-modules:      Data.Packed.Internal,
                         Data.Packed.Internal.Common,
@@ -50,6 +50,11 @@ library
                         Data.Packed.Internal.Vector,
                         Data.Packed.Internal.Matrix
                         Numeric.Chain
+                        Numeric.Vectorized
+                        Numeric.Vector
+                        Numeric.Matrix
+                        Data.Packed.Numeric
+
 
     C-sources:          src/C/lapack-aux.c
                         src/C/vector-aux.c
diff --git a/packages/base/src/Data/Packed/IO.hs b/packages/base/src/Data/Packed/IO.hs
new file mode 100644
index 0000000..dbb2943
--- /dev/null
+++ b/packages/base/src/Data/Packed/IO.hs
@@ -0,0 +1,141 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.IO
+-- Copyright   :  (c) Alberto Ruiz 2010
+-- License     :  BSD3
+--
+-- Maintainer  :  Alberto Ruiz
+-- Stability   :  provisional
+--
+-- Display, formatting and IO functions for numeric 'Vector' and 'Matrix'
+--
+-----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
+module Data.Packed.IO (
+    dispf, disps, dispcf, vecdisp, latexFormat, format,
+    readMatrix, fromArray2D
+) where
+
+import Data.Packed
+import Data.Packed.Development
+import Text.Printf(printf)
+import Data.List(intersperse)
+import Data.Complex
+
+{- | Creates a string from a matrix given a separator and a function to show each entry. Using
+this function the user can easily define any desired display function:
+
+@import Text.Printf(printf)@
+
+@disp = putStr . format \"  \" (printf \"%.2f\")@
+
+-}
+format :: (Element t) => String -> (t -> String) -> Matrix t -> String
+format sep f m = table sep . map (map f) . toLists $ m
+
+{- | Show a matrix with \"autoscaling\" and a given number of decimal places.
+
+>>> putStr . disps 2 $ 120 * (3><4) [1..]
+3x4  E3
+ 0.12  0.24  0.36  0.48
+ 0.60  0.72  0.84  0.96
+ 1.08  1.20  1.32  1.44
+
+-}
+disps :: Int -> Matrix Double -> String
+disps d x = sdims x ++ "  " ++ formatScaled d x
+
+{- | Show a matrix with a given number of decimal places.
+
+>>> dispf 2 (1/3 + ident 3)
+"3x3\n1.33  0.33  0.33\n0.33  1.33  0.33\n0.33  0.33  1.33\n"
+
+>>> putStr . dispf 2 $ (3><4)[1,1.5..]
+3x4
+1.00  1.50  2.00  2.50
+3.00  3.50  4.00  4.50
+5.00  5.50  6.00  6.50
+
+>>> putStr . unlines . tail . lines . dispf 2 . asRow $ linspace 10 (0,1)
+0.00  0.11  0.22  0.33  0.44  0.56  0.67  0.78  0.89  1.00
+
+-}
+dispf :: Int -> Matrix Double -> String
+dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x
+
+sdims x = show (rows x) ++ "x" ++ show (cols x)
+
+formatFixed d x = format "  " (printf ("%."++show d++"f")) $ x
+
+isInt = all lookslikeInt . toList . flatten
+
+formatScaled dec t = "E"++show o++"\n" ++ ss
+    where ss = format " " (printf fmt. g) t
+          g x | o >= 0    = x/10^(o::Int)
+              | otherwise = x*10^(-o)
+          o | rows t == 0 || cols t == 0 = 0
+            | otherwise = floor $ maximum $ map (logBase 10 . abs) $ toList $ flatten t
+          fmt = '%':show (dec+3) ++ '.':show dec ++"f"
+
+{- | Show a vector using a function for showing matrices.
+
+>>> putStr . vecdisp (dispf 2) $ linspace 10 (0,1)
+10 |> 0.00  0.11  0.22  0.33  0.44  0.56  0.67  0.78  0.89  1.00
+
+-}
+vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String
+vecdisp f v
+    = ((show (dim v) ++ " |> ") ++) . (++"\n")
+    . unwords . lines .  tail . dropWhile (not . (`elem` " \n"))
+    . f . trans . reshape 1
+    $ v
+
+{- | Tool to display matrices with latex syntax.
+
+>>>  latexFormat "bmatrix" (dispf 2 $ ident 2)
+"\\begin{bmatrix}\n1  &  0\n\\\\\n0  &  1\n\\end{bmatrix}"
+
+-}
+latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc.
+            -> String -- ^ Formatted matrix, with elements separated by spaces and newlines
+            -> String
+latexFormat del tab = "\\begin{"++del++"}\n" ++ f tab ++ "\\end{"++del++"}"
+    where f = unlines . intersperse "\\\\" . map unwords . map (intersperse " & " . words) . tail . lines
+
+-- | Pretty print a complex number with at most n decimal digits.
+showComplex :: Int -> Complex Double -> String
+showComplex d (a:+b)
+    | isZero a && isZero b = "0"
+    | isZero b = sa
+    | isZero a && isOne b = s2++"i"
+    | isZero a = sb++"i"
+    | isOne b = sa++s3++"i"
+    | otherwise = sa++s1++sb++"i"
+  where
+    sa = shcr d a
+    sb = shcr d b
+    s1 = if b<0 then "" else "+"
+    s2 = if b<0 then "-" else ""
+    s3 = if b<0 then "-" else "+"
+
+shcr d a | lookslikeInt a = printf "%.0f" a
+         | otherwise      = printf ("%."++show d++"f") a
+
+
+lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
+   where shx = show x
+
+isZero x = show x `elem` ["0.0","-0.0"]
+isOne  x = show x `elem` ["1.0","-1.0"]
+
+-- | Pretty print a complex matrix with at most n decimal digits.
+dispcf :: Int -> Matrix (Complex Double) -> String
+dispcf d m = sdims m ++ "\n" ++ format "  " (showComplex d) m
+
+--------------------------------------------------------------------
+
+-- | reads a matrix from a string containing a table of numbers.
+readMatrix :: String -> Matrix Double
+readMatrix = fromLists . map (map read). map words . filter (not.null) . lines
+
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
new file mode 100644
index 0000000..1db860c
--- /dev/null
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -0,0 +1,30 @@
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra
+Copyright   :  (c) Alberto Ruiz 2006-10
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+This module reexports all normally required functions for Linear Algebra applications.
+
+It also provides instances of standard classes 'Show', 'Read', 'Eq',
+'Num', 'Fractional', and 'Floating' for 'Vector' and 'Matrix'.
+In arithmetic operations one-component vectors and matrices automatically
+expand to match the dimensions of the other operand.
+
+-}
+-----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
+module Numeric.LinearAlgebra (
+    module Numeric.Container,
+    module Numeric.LinearAlgebra.Algorithms
+) where
+
+import Numeric.Container
+import Numeric.LinearAlgebra.Algorithms
+import Numeric.Matrix()
+import Numeric.Vector()
diff --git a/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs b/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
index af0c744..9cb67d4 100644
--- a/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs
@@ -10,6 +10,8 @@
 -- Functional interface to selected LAPACK functions (<http://www.netlib.org/lapack>).
 --
 -----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
 
 module Numeric.LinearAlgebra.LAPACK (
     -- * Matrix product
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs
new file mode 100644
index 0000000..533c54b
--- /dev/null
+++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs
@@ -0,0 +1,291 @@
+{-# LANGUAGE FlexibleContexts #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra.Util
+Copyright   :  (c) Alberto Ruiz 2013
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+
+-}
+-----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
+module Numeric.LinearAlgebra.Util(
+
+    -- * Convenience functions
+    size, disp,
+    zeros, ones,
+    diagl,
+    row,
+    col,
+    (&), (¦), (——), (#),
+    (?), (¿),
+    cross,
+    norm,
+    unitary,
+    mt,
+    pairwiseD2,
+    meanCov,
+    rowOuters,
+    null1,
+    null1sym,
+    -- * 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,
+{-    -- * Plot
+    mplot,
+    plot, parametricPlot,
+    splot, mesh, meshdom,
+    matrixToPGM, imshow,
+    gnuplotX, gnuplotpdf, gnuplotWin
+-}
+) where
+
+import Numeric.Container
+import Data.Packed.IO
+import Numeric.LinearAlgebra.Algorithms hiding (i)
+import Numeric.Matrix()
+import Numeric.Vector()
+
+import Numeric.LinearAlgebra.Util.Convolution
+--import Graphics.Plot
+
+
+{- | print a real matrix with given number of digits after the decimal point
+
+>>> disp 5 $ ident 2 / 3
+2x2
+0.33333  0.00000
+0.00000  0.33333
+
+-}
+disp :: Int -> Matrix Double -> IO ()
+
+disp n = putStrLn . dispf n
+
+
+{- | create a real diagonal matrix from a list
+
+>>> diagl [1,2,3]
+(3><3)
+ [ 1.0, 0.0, 0.0
+ , 0.0, 2.0, 0.0
+ , 0.0, 0.0, 3.0 ]
+
+-}
+diagl :: [Double] -> Matrix Double
+diagl = diag . fromList
+
+-- | a real matrix of zeros
+zeros :: Int -- ^ rows
+      -> Int -- ^ columns
+      -> Matrix Double
+zeros r c = konst 0 (r,c)
+
+-- | a real matrix of ones
+ones :: Int -- ^ rows
+     -> Int -- ^ columns
+     -> Matrix Double
+ones r c = konst 1 (r,c)
+
+-- | concatenation of real vectors
+infixl 3 &
+(&) :: Vector Double -> Vector Double -> Vector Double
+a & b = vjoin [a,b]
+
+{- | horizontal concatenation of real matrices
+
+ (unicode 0x00a6, broken bar)
+
+>>> ident 3 ¦ konst 7 (3,4)
+(3><7)
+ [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0
+ , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0
+ , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ]
+
+-}
+infixl 3 ¦
+(¦) :: Matrix Double -> Matrix Double -> Matrix Double
+a ¦ b = fromBlocks [[a,b]]
+
+-- | vertical concatenation of real matrices
+--
+-- (unicode 0x2014, em dash)
+(——) :: Matrix Double -> Matrix Double -> Matrix Double
+infixl 2 ——
+a —— b = fromBlocks [[a],[b]]
+
+(#) :: Matrix Double -> Matrix Double -> Matrix Double
+infixl 2 #
+a # b = fromBlocks [[a],[b]]
+
+-- | create a single row real matrix from a list
+row :: [Double] -> Matrix Double
+row = asRow . fromList
+
+-- | create a single column real matrix from a list
+col :: [Double] -> Matrix Double
+col = asColumn . fromList
+
+{- | extract rows
+
+>>> (20><4) [1..] ? [2,1,1]
+(3><4)
+ [ 9.0, 10.0, 11.0, 12.0
+ , 5.0,  6.0,  7.0,  8.0
+ , 5.0,  6.0,  7.0,  8.0 ]
+
+-}
+infixl 9 ?
+(?) :: Element t => Matrix t -> [Int] -> Matrix t
+(?) = flip extractRows
+
+{- | extract columns
+
+(unicode 0x00bf, inverted question mark, Alt-Gr ?)
+
+>>> (3><4) [1..] ¿ [3,0]
+(3><2)
+ [  4.0, 1.0
+ ,  8.0, 5.0
+ , 12.0, 9.0 ]
+
+-}
+infixl 9 ¿
+(¿) :: Element t => Matrix t -> [Int] -> Matrix t
+(¿)= flip extractColumns
+
+
+cross :: Vector Double -> Vector Double -> Vector Double
+-- ^ 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
+    [x1,x2,x3] = toList x
+    [y1,y2,y3] = toList y
+    z1 = x2*y3-x3*y2
+    z2 = x3*y1-x1*y3
+    z3 = x1*y2-x2*y1
+
+norm :: Vector Double -> Double
+-- ^ 2-norm of real vector
+norm = pnorm PNorm2
+
+
+-- | Obtains a vector in the same direction with 2-norm=1
+unitary :: Vector Double -> Vector Double
+unitary v = v / scalar (norm v)
+
+-- | ('rows' &&& 'cols')
+size :: Matrix t -> (Int, Int)
+size m = (rows m, cols m)
+
+-- | trans . inv
+mt :: Matrix Double -> Matrix Double
+mt = trans . inv
+
+--------------------------------------------------------------------------------
+
+{- | Compute mean vector and covariance matrix of the rows of a matrix.
+
+>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3])
+(fromList [4.010341078059521,5.0197204699640405],
+(2><2)
+ [     1.9862461923890056, -1.0127225830525157e-2
+ , -1.0127225830525157e-2,     3.0373954915729318 ])
+
+-}
+meanCov :: Matrix Double -> (Vector Double, Matrix Double)
+meanCov x = (med,cov) where
+    r    = rows x
+    k    = 1 / fromIntegral r
+    med  = konst k r `vXm` x
+    meds = konst 1 r `outer` med
+    xc   = x `sub` meds
+    cov  = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc)
+
+--------------------------------------------------------------------------------
+
+-- | Matrix of pairwise squared distances of row vectors
+-- (using the matrix product trick in blog.smola.org)
+pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double
+pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y
+               | otherwise = error $ "pairwiseD2 with different number of columns: "
+                                   ++ show (size x) ++ ", " ++ show (size y)
+  where
+    ox = one (rows x)
+    oy = one (rows y)
+    oc = one (cols x)
+    one k = constant 1 k
+    x2 = x * x <> oc
+    y2 = y * y <> oc
+    ok = cols x == cols y
+
+--------------------------------------------------------------------------------
+
+-- | outer products of rows
+rowOuters :: Matrix Double -> Matrix Double -> Matrix Double
+rowOuters a b = a' * b'
+  where
+    a' = kronecker a (ones 1 (cols b))
+    b' = kronecker (ones 1 (cols a)) b
+
+--------------------------------------------------------------------------------
+
+-- | solution of overconstrained homogeneous linear system
+null1 :: Matrix Double -> Vector Double
+null1 = last . toColumns . snd . rightSV
+
+-- | solution of overconstrained homogeneous symmetric linear system
+null1sym :: Matrix Double -> Vector Double
+null1sym = last . toColumns . snd . eigSH'
+
+--------------------------------------------------------------------------------
+
+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 = vjoin . 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
+
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs
new file mode 100644
index 0000000..d04c46b
--- /dev/null
+++ b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs
@@ -0,0 +1,115 @@
+{-# LANGUAGE FlexibleContexts #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  Numeric.LinearAlgebra.Util.Convolution
+Copyright   :  (c) Alberto Ruiz 2012
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+
+-}
+-----------------------------------------------------------------------------
+{-# OPTIONS_HADDOCK hide #-}
+
+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 = vjoin [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 (udot 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, Container Vector 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
+