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/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 +- 7 files changed, 5 insertions(+), 561 deletions(-) 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/hmatrix') 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