From 1838c4248679b7476bb8716a76171712dc3cd335 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 16 May 2014 13:35:35 +0200 Subject: linear algebra moved --- packages/base/hmatrix-base.cabal | 13 +- packages/base/src/Data/Packed/IO.hs | 141 ++++++++++ packages/base/src/Numeric/LinearAlgebra.hs | 30 +++ packages/base/src/Numeric/LinearAlgebra/LAPACK.hs | 2 + packages/base/src/Numeric/LinearAlgebra/Util.hs | 291 +++++++++++++++++++++ .../src/Numeric/LinearAlgebra/Util/Convolution.hs | 115 ++++++++ 6 files changed, 588 insertions(+), 4 deletions(-) create mode 100644 packages/base/src/Data/Packed/IO.hs create mode 100644 packages/base/src/Numeric/LinearAlgebra.hs create mode 100644 packages/base/src/Numeric/LinearAlgebra/Util.hs create mode 100644 packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs (limited to 'packages/base') 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 (). -- ----------------------------------------------------------------------------- +{-# 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 + -- cgit v1.2.3