From 2876998f04380c9e835c6177b440447368dfe623 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 5 Jun 2015 16:35:18 +0200 Subject: move util --- packages/base/src/Numeric/LinearAlgebra/Util.hs | 560 ------------------------ 1 file changed, 560 deletions(-) delete mode 100644 packages/base/src/Numeric/LinearAlgebra/Util.hs (limited to 'packages/base/src/Numeric/LinearAlgebra/Util.hs') diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs deleted file mode 100644 index 37cdc0f..0000000 --- a/packages/base/src/Numeric/LinearAlgebra/Util.hs +++ /dev/null @@ -1,560 +0,0 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FunctionalDependencies #-} -{-# LANGUAGE ViewPatterns #-} - - ------------------------------------------------------------------------------ -{- | -Module : Numeric.LinearAlgebra.Util -Copyright : (c) Alberto Ruiz 2013 -License : BSD3 -Maintainer : Alberto Ruiz -Stability : provisional - --} ------------------------------------------------------------------------------ - -module Numeric.LinearAlgebra.Util( - - -- * Convenience functions - vector, matrix, - disp, - formatSparse, - approxInt, - dispDots, - dispBlanks, - formatShort, - dispShort, - zeros, ones, - diagl, - row, - col, - (&), (¦), (|||), (——), (===), (#), - (?), (¿), - Indexable(..), size, - Numeric, - rand, randn, - cross, - norm, - ℕ,ℤ,ℝ,ℂ,iC, - Normed(..), norm_Frob, norm_nuclear, - unitary, - mt, - (~!~), - pairwiseD2, - rowOuters, - null1, - null1sym, - -- * Convolution - -- ** 1D - corr, conv, corrMin, - -- ** 2D - corr2, conv2, separable, - gaussElim -) where - -import Data.Packed.Numeric -import Numeric.LinearAlgebra.Algorithms hiding (i,Normed) ---import qualified Numeric.LinearAlgebra.Algorithms as A -import Numeric.Matrix() -import Numeric.Vector() -import Numeric.LinearAlgebra.Random -import Numeric.LinearAlgebra.Util.Convolution -import Control.Monad(when) -import Text.Printf -import Data.List.Split(splitOn) -import Data.List(intercalate,sortBy) -import Data.Function(on) -import Control.Arrow((&&&)) - -type ℝ = Double -type ℕ = Int -type ℤ = Int -type ℂ = Complex Double - --- | imaginary unit -iC :: ℂ -iC = 0:+1 - -{- | Create a real vector. - ->>> vector [1..5] -fromList [1.0,2.0,3.0,4.0,5.0] - --} -vector :: [ℝ] -> Vector ℝ -vector = fromList - -{- | Create a real matrix. - ->>> matrix 5 [1..15] -(3><5) - [ 1.0, 2.0, 3.0, 4.0, 5.0 - , 6.0, 7.0, 8.0, 9.0, 10.0 - , 11.0, 12.0, 13.0, 14.0, 15.0 ] - --} -matrix - :: Int -- ^ number of columns - -> [ℝ] -- ^ elements in row order - -> Matrix ℝ -matrix c = reshape c . fromList - - -{- | 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 = putStr . 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 - ->>> 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]] - --- | a synonym for ('|||') (unicode 0x00a6, broken bar) -infixl 3 ¦ -(¦) :: Matrix Double -> Matrix Double -> Matrix Double -(¦) = (|||) - - --- | vertical concatenation of real matrices --- -(===) :: Matrix Double -> Matrix Double -> Matrix Double -infixl 2 === -a === b = fromBlocks [[a],[b]] - --- | a synonym for ('===') (unicode 0x2014, em dash) -(——) :: Matrix Double -> Matrix Double -> Matrix Double -infixl 2 —— -(——) = (===) - - -(#) :: Matrix Double -> Matrix Double -> Matrix Double -infixl 2 # -a # b = fromBlocks [[a],[b]] - --- | create a single row real matrix from a list --- --- >>> row [2,3,1,8] --- (1><4) --- [ 2.0, 3.0, 1.0, 8.0 ] --- -row :: [Double] -> Matrix Double -row = asRow . fromList - --- | create a single column real matrix from a list --- --- >>> col [7,-2,4] --- (3><1) --- [ 7.0 --- , -2.0 --- , 4.0 ] --- -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 :: Product t => Vector t -> Vector t -> Vector t --- ^ cross product (for three-element vectors) -cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3] - | otherwise = error $ "the cross product requires 3-element vectors (sizes given: " - ++show (dim x)++" and "++show (dim 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 - -{-# SPECIALIZE cross :: Vector Double -> Vector Double -> Vector Double #-} -{-# SPECIALIZE cross :: Vector (Complex Double) -> Vector (Complex Double) -> Vector (Complex Double) #-} - -norm :: Vector Double -> Double --- ^ 2-norm of real vector -norm = pnorm PNorm2 - -class Normed a - where - norm_0 :: a -> ℝ - norm_1 :: a -> ℝ - norm_2 :: a -> ℝ - norm_Inf :: a -> ℝ - - -instance Normed (Vector ℝ) - where - norm_0 v = sumElements (step (abs v - scalar (eps*normInf v))) - norm_1 = pnorm PNorm1 - norm_2 = pnorm PNorm2 - norm_Inf = pnorm Infinity - -instance Normed (Vector ℂ) - where - norm_0 v = sumElements (step (fst (fromComplex (abs v)) - scalar (eps*normInf v))) - norm_1 = pnorm PNorm1 - norm_2 = pnorm PNorm2 - norm_Inf = pnorm Infinity - -instance Normed (Matrix ℝ) - where - norm_0 = norm_0 . flatten - norm_1 = pnorm PNorm1 - norm_2 = pnorm PNorm2 - norm_Inf = pnorm Infinity - -instance Normed (Matrix ℂ) - where - norm_0 = norm_0 . flatten - norm_1 = pnorm PNorm1 - norm_2 = pnorm PNorm2 - norm_Inf = pnorm Infinity - -instance Normed (Vector I) - where - norm_0 = fromIntegral . sumElements . step . abs - norm_1 = fromIntegral . norm1 - norm_2 v = sqrt . fromIntegral $ dot v v - norm_Inf = fromIntegral . normInf - - - -norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> ℝ -norm_Frob = norm_2 . flatten - -norm_nuclear :: Field t => Matrix t -> ℝ -norm_nuclear = sumElements . singularValues - - --- | Obtains a vector in the same direction with 2-norm=1 -unitary :: Vector Double -> Vector Double -unitary v = v / scalar (norm v) - - --- | trans . inv -mt :: Matrix Double -> Matrix Double -mt = trans . inv - --------------------------------------------------------------------------------- -{- | - ->>> size $ vector [1..10] -10 ->>> size $ (2><5)[1..10::Double] -(2,5) - --} -size :: Container c t => c t -> IndexOf c -size = size' - -{- | Alternative indexing function. - ->>> vector [1..10] ! 3 -4.0 - -On a matrix it gets the k-th row as a vector: - ->>> matrix 5 [1..15] ! 1 -fromList [6.0,7.0,8.0,9.0,10.0] - ->>> matrix 5 [1..15] ! 1 ! 3 -9.0 - --} -class Indexable c t | c -> t , t -> c - where - infixl 9 ! - (!) :: c -> Int -> t - -instance Indexable (Vector Double) Double - where - (!) = (@>) - -instance Indexable (Vector Float) Float - where - (!) = (@>) - -instance Indexable (Vector I) I - where - (!) = (@>) - -instance Indexable (Vector (Complex Double)) (Complex Double) - where - (!) = (@>) - -instance Indexable (Vector (Complex Float)) (Complex Float) - where - (!) = (@>) - -instance Element t => Indexable (Matrix t) (Vector t) - where - m!j = subVector (j*c) c (flatten m) - where - c = cols m - --------------------------------------------------------------------------------- - --- | 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 = konst 1 k - x2 = x * x <> oc - y2 = y * y <> oc - ok = cols x == cols y - --------------------------------------------------------------------------------- - -{- | outer products of rows - ->>> a -(3><2) - [ 1.0, 2.0 - , 10.0, 20.0 - , 100.0, 200.0 ] ->>> b -(3><3) - [ 1.0, 2.0, 3.0 - , 4.0, 5.0, 6.0 - , 7.0, 8.0, 9.0 ] - ->>> rowOuters a (b ||| 1) -(3><8) - [ 1.0, 2.0, 3.0, 1.0, 2.0, 4.0, 6.0, 2.0 - , 40.0, 50.0, 60.0, 10.0, 80.0, 100.0, 120.0, 20.0 - , 700.0, 800.0, 900.0, 100.0, 1400.0, 1600.0, 1800.0, 200.0 ] - --} -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' - --------------------------------------------------------------------------------- - -infixl 0 ~!~ -c ~!~ msg = when c (error msg) - --------------------------------------------------------------------------------- - -formatSparse :: String -> String -> String -> Int -> Matrix Double -> String - -formatSparse zeroI _zeroF sep _ (approxInt -> Just m) = format sep f m - where - f 0 = zeroI - f x = printf "%.0f" x - -formatSparse zeroI zeroF sep n m = format sep f m - where - f x | abs (x::Double) < 2*peps = zeroI++zeroF - | abs (fromIntegral (round x::Int) - x) / abs x < 2*peps - = printf ("%.0f."++replicate n ' ') x - | otherwise = printf ("%."++show n++"f") x - -approxInt m - | norm_Inf (v - vi) < 2*peps * norm_Inf v = Just (reshape (cols m) vi) - | otherwise = Nothing - where - v = flatten m - vi = roundVector v - -dispDots n = putStr . formatSparse "." (replicate n ' ') " " n - -dispBlanks n = putStr . formatSparse "" "" " " n - -formatShort sep fmt maxr maxc m = auxm4 - where - (rm,cm) = size m - (r1,r2,r3) - | rm <= maxr = (rm,0,0) - | otherwise = (maxr-3,rm-maxr+1,2) - (c1,c2,c3) - | cm <= maxc = (cm,0,0) - | otherwise = (maxc-3,cm-maxc+1,2) - [ [a,_,b] - ,[_,_,_] - ,[c,_,d]] = toBlocks [r1,r2,r3] - [c1,c2,c3] m - auxm = fromBlocks [[a,b],[c,d]] - auxm2 - | cm > maxc = format "|" fmt auxm - | otherwise = format sep fmt auxm - auxm3 - | cm > maxc = map (f . splitOn "|") (lines auxm2) - | otherwise = (lines auxm2) - f items = intercalate sep (take (maxc-3) items) ++ " .. " ++ - intercalate sep (drop (maxc-3) items) - auxm4 - | rm > maxr = unlines (take (maxr-3) auxm3 ++ vsep : drop (maxr-3) auxm3) - | otherwise = unlines auxm3 - vsep = map g (head auxm3) - g '.' = ':' - g _ = ' ' - - -dispShort :: Int -> Int -> Int -> Matrix Double -> IO () -dispShort maxr maxc dec m = - printf "%dx%d\n%s" (rows m) (cols m) (formatShort " " fmt maxr maxc m) - where - fmt = printf ("%."++show dec ++"f") - --------------------------------------------------------------------------------- - --- | generic reference implementation of gaussian elimination --- --- @a <> gauss a b = b@ --- -gaussElim - :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t) - => Matrix t -> Matrix t -> Matrix t - -gaussElim x y = dropColumns (rows x) (flipud $ fromRows s2) - where - rs = toRows $ fromBlocks [[x , y]] - s1 = pivotDown (rows x) 0 rs - s2 = pivotUp (rows x-1) (reverse s1) - -pivotDown t n xs - | t == n = [] - | otherwise = y : pivotDown t (n+1) ys - where - y:ys = redu (pivot n xs) - - pivot k = (const k &&& id) - . reverse . sortBy (compare `on` (abs. (!k))) -- FIXME - - redu (k,x:zs) - | p == 0 = error "gauss: singular!" -- FIXME - | otherwise = u : map f zs - where - p = x!k - u = scale (recip (x!k)) x - f z = z - scale (z!k) u - redu (_,[]) = [] - - -pivotUp n xs - | n == -1 = [] - | otherwise = y : pivotUp (n-1) ys - where - y:ys = redu' (n,xs) - - redu' (k,x:zs) = u : map f zs - where - u = x - f z = z - scale (z!k) u - redu' (_,[]) = [] - --------------------------------------------------------------------------------- - -instance Testable (Matrix I) where - checkT _ = test - -test :: (Bool, IO()) -test = (and ok, return ()) - where - m = (3><4) [1..12] :: Matrix I - r = (2><3) [1,2,3,4,3,2] - c = (3><2) [0,4,4,1,2,3] - p = (9><10) [0..89] :: Matrix I - ep = (2><3) [10,24,32,44,31,23] - md = fromInt m :: Matrix Double - ok = [ tr m <> m == toInt (tr md <> md) - , m <> tr m == toInt (md <> tr md) - , m ?? (Take 2, Take 3) == remap (asColumn (range 2)) (asRow (range 3)) m - , remap r (tr c) p == ep - , tr p ?? (PosCyc (idxs[-5,13]), Pos (idxs[3,7,1])) == (2><3) [35,75,15,33,73,13] - ] - -- cgit v1.2.3