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 ++++++++ packages/hmatrix/hmatrix.cabal | 5 +- packages/hmatrix/src/Numeric/HMatrix.hs | 4 +- packages/hmatrix/src/Numeric/IO.hs | 120 +-------- packages/hmatrix/src/Numeric/LinearAlgebra.hs | 30 --- packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs | 290 -------------------- .../src/Numeric/LinearAlgebra/Util/Convolution.hs | 114 -------- packages/hmatrix/src/Numeric/Random.hs | 3 +- 13 files changed, 593 insertions(+), 565 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 delete mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra.hs delete mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs delete mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs (limited to 'packages') 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 + diff --git a/packages/hmatrix/hmatrix.cabal b/packages/hmatrix/hmatrix.cabal index 9719c96..c79da8a 100644 --- a/packages/hmatrix/hmatrix.cabal +++ b/packages/hmatrix/hmatrix.cabal @@ -90,8 +90,6 @@ library Numeric.GSL.Fitting, Numeric.GSL.ODE, Numeric.GSL, - Numeric.LinearAlgebra, - Numeric.LinearAlgebra.Util, Graphics.Plot, Numeric.HMatrix, Numeric.HMatrix.Data, @@ -99,8 +97,7 @@ library other-modules: Numeric.Random, Numeric.GSL.Internal, Numeric.GSL.Vector, - Numeric.IO, - Numeric.LinearAlgebra.Util.Convolution + Numeric.IO C-sources: src/Numeric/GSL/gsl-aux.c diff --git a/packages/hmatrix/src/Numeric/HMatrix.hs b/packages/hmatrix/src/Numeric/HMatrix.hs index 36fcf70..fcd3e02 100644 --- a/packages/hmatrix/src/Numeric/HMatrix.hs +++ b/packages/hmatrix/src/Numeric/HMatrix.hs @@ -129,8 +129,8 @@ module Numeric.HMatrix ( import Numeric.HMatrix.Data -import Numeric.Matrix() -import Numeric.Vector() +--import Numeric.Matrix() +--import Numeric.Vector() import Numeric.Container import Numeric.LinearAlgebra.Algorithms import Numeric.LinearAlgebra.Util diff --git a/packages/hmatrix/src/Numeric/IO.hs b/packages/hmatrix/src/Numeric/IO.hs index 58fa2b4..08d812b 100644 --- a/packages/hmatrix/src/Numeric/IO.hs +++ b/packages/hmatrix/src/Numeric/IO.hs @@ -20,128 +20,10 @@ module Numeric.IO ( ) where import Data.Packed -import Data.Packed.Development +import Data.Packed.IO import System.Process(readProcess) -import Text.Printf(printf) -import Data.List(intersperse) -import Data.Complex import Numeric.GSL.Vector -{- | 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 {- | obtains the number of rows and columns in an ASCII data file (provisionally using unix's wc). diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra.hs b/packages/hmatrix/src/Numeric/LinearAlgebra.hs deleted file mode 100644 index 1db860c..0000000 --- a/packages/hmatrix/src/Numeric/LinearAlgebra.hs +++ /dev/null @@ -1,30 +0,0 @@ ------------------------------------------------------------------------------ -{- | -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/hmatrix/src/Numeric/LinearAlgebra/Util.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs deleted file mode 100644 index e4f21b0..0000000 --- a/packages/hmatrix/src/Numeric/LinearAlgebra/Util.hs +++ /dev/null @@ -1,290 +0,0 @@ -{-# 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 Numeric.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/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs deleted file mode 100644 index 82de476..0000000 --- a/packages/hmatrix/src/Numeric/LinearAlgebra/Util/Convolution.hs +++ /dev/null @@ -1,114 +0,0 @@ -{-# 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 = 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 - diff --git a/packages/hmatrix/src/Numeric/Random.hs b/packages/hmatrix/src/Numeric/Random.hs index 7bf9e8b..0e722fd 100644 --- a/packages/hmatrix/src/Numeric/Random.hs +++ b/packages/hmatrix/src/Numeric/Random.hs @@ -21,8 +21,7 @@ module Numeric.Random ( ) where import Numeric.GSL.Vector -import Data.Packed -import Data.Packed.Numeric +import Numeric.Container import Numeric.LinearAlgebra.Algorithms import System.Random(randomIO) -- cgit v1.2.3