hmatrix 0.8.1.1 → 0.8.2.0
raw patch · 17 files changed
+282/−86 lines, 17 filesdep ~arraydep ~haskell98dep ~process
Dependency ranges changed: array, haskell98, process
Files
- CHANGES +14/−0
- README +2/−0
- examples/lie.hs +3/−4
- examples/parallel.hs +3/−0
- hmatrix.cabal +6/−6
- lib/Data/Packed/Internal/Common.hs +7/−1
- lib/Data/Packed/Internal/Matrix.hs +8/−3
- lib/Data/Packed/Matrix.hs +114/−16
- lib/Data/Packed/Random.hs +2/−1
- lib/Data/Packed/Vector.hs +21/−0
- lib/Numeric/GSL/Minimization.hs +56/−37
- lib/Numeric/LinearAlgebra/Algorithms.hs +12/−7
- lib/Numeric/LinearAlgebra/Instances.hs +4/−4
- lib/Numeric/LinearAlgebra/Interface.hs +4/−0
- lib/Numeric/LinearAlgebra/Linear.hs +10/−3
- lib/Numeric/LinearAlgebra/Tests.hs +15/−3
- lib/Numeric/LinearAlgebra/Tests/Instances.hs +1/−1
CHANGES view
@@ -1,3 +1,17 @@+0.8.2.0+=======++- fromRows/fromColumns now automatically expand vectors of dim 1+ to match the common dimension.+ fromBlocks also replicates single row/column matrices.+ Previously all dimensions had to be exactly the same.++- display utilities: dispf, disps, vecdisp++- scalar++- minimizeV, minimizeVD, using Vector instead of lists.+ 0.8.1.0 =======
README view
@@ -106,3 +106,5 @@ - Erik de Castro Lopo added buildVector and buildMatrix, which take a size parameter(s) and a function that maps vector/matrix indices to the values at that position.++- Jean-Francois Tremblay discovered an error in the tutorial.
examples/lie.hs view
@@ -1,9 +1,8 @@ -- The magic of Lie Algebra import Numeric.LinearAlgebra-import Text.Printf(printf) -disp = putStrLn . format " " (printf "%.5f")+disp = putStrLn . dispf 5 rot1 :: Double -> Matrix Double rot1 a = (3><3)@@ -50,14 +49,14 @@ a & b = a <> b - b <> a infixl 6 |+|-a |+| b = a + b + a & b */2 + (a-b) & (a & b) */12+a |+| b = a + b + a&b /2 + (a-b)&(a & b) /12 main = do let a = 45*deg b = 50*deg c = -30*deg exact = rot3 a <> rot1 b <> rot2 c- lie = a.*g3 |+| b.*g1 |+| c.*g2+ lie = scalar a * g3 |+| scalar b * g1 |+| scalar c * g2 putStrLn "position in the tangent space:" disp lie putStrLn "exponential map back to the group (2 terms):"
examples/parallel.hs view
@@ -12,10 +12,13 @@ parMul p x y = fromBlocks [ inParallel ( map (x <>) ys ) ] where ys = splitColumns p y +-- x <||> y = fromColumns . inParallel . map (x <>) . toColumns $ y+ main = do n <- (read . head) `fmap` getArgs let m = ident n :: Matrix Double time $ print $ vectorMax $ takeDiag $ m <> m+-- time $ print $ vectorMax $ takeDiag $ m <||> m time $ print $ vectorMax $ takeDiag $ parMul 2 m m time $ print $ vectorMax $ takeDiag $ parMul 4 m m time $ print $ vectorMax $ takeDiag $ parMul 8 m m
hmatrix.cabal view
@@ -1,12 +1,12 @@ Name: hmatrix-Version: 0.8.1.1+Version: 0.8.2.0 License: GPL License-file: LICENSE Author: Alberto Ruiz Maintainer: Alberto Ruiz <aruiz@um.es> Stability: provisional-Homepage: http://www.hmatrix.googlepages.com-Synopsis: Linear algebra and numerical computations+Homepage: http://code.haskell.org/hmatrix+Synopsis: Linear algebra and numerical computation Description: Purely functional interface to basic linear algebra and other numerical computations, internally implemented using GSL, BLAS and LAPACK.@@ -64,14 +64,14 @@ library if flag(splitBase)- build-depends: base >= 3 && < 5, array < 0.3.0+ build-depends: base >= 3 && < 5, array else build-depends: base < 3 - Build-Depends: haskell98 < 1.0.1.1,+ Build-Depends: haskell98, QuickCheck, HUnit, storable-complex,- process < 1.0.1.2+ process Extensions: ForeignFunctionInterface, CPP
lib/Data/Packed/Internal/Common.hs view
@@ -18,7 +18,7 @@ Adapt, app1, app2, app3, app4, (//), check,- partit, common,+ partit, common, compatdim, fi ) where @@ -40,6 +40,12 @@ commonval [] = Nothing commonval [a] = Just a commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing++-- | common value with \"adaptable\" 1+compatdim :: [Int] -> Maybe Int+compatdim [] = Nothing+compatdim [a] = Just a+compatdim (a:b:xs) = if a==b || a==1 || b==1 then compatdim (max a b:xs) else Nothing -- | postfix function application (@flip ($)@) (//) :: x -> (x -> y) -> y
lib/Data/Packed/Internal/Matrix.hs view
@@ -135,11 +135,16 @@ toLists :: (Element t) => Matrix t -> [[t]] toLists m = partit (cols m) . toList . flatten $ m --- | creates a Matrix from a list of vectors+-- | Create a matrix from a list of vectors.+-- All vectors must have the same dimension,+-- or dimension 1, which is are automatically expanded. fromRows :: Element t => [Vector t] -> Matrix t-fromRows vs = case common dim vs of+fromRows vs = case compatdim (map dim vs) of Nothing -> error "fromRows applied to [] or to vectors with different sizes"- Just c -> reshape c (join vs)+ Just c -> reshape c . join . map (adapt c) $ vs+ where+ adapt c v | dim v == c = v+ | otherwise = constantD (v@>0) c -- | extracts the rows of a matrix as a list of vectors toRows :: Element t => Matrix t -> [Vector t]
lib/Data/Packed/Matrix.hs view
@@ -1,4 +1,3 @@-{-# OPTIONS_GHC -fglasgow-exts #-} ----------------------------------------------------------------------------- -- | -- Module : Data.Packed.Matrix@@ -29,7 +28,7 @@ extractRows, ident, diag, diagRect, takeDiag, liftMatrix, liftMatrix2,- format,+ format, dispf, disps, loadMatrix, saveMatrix, fromFile, fileDimensions, readMatrix, fromArray2D ) where@@ -40,33 +39,67 @@ import Data.List(transpose,intersperse) import Data.Array import System.Process(readProcess)+import Text.Printf(printf) -- | creates a matrix from a vertical list of matrices joinVert :: Element t => [Matrix t] -> Matrix t joinVert ms = case common cols ms of- Nothing -> error "joinVert on matrices with different number of columns"+ Nothing -> error "(impossible) joinVert on matrices with different number of columns" Just c -> reshape c $ join (map flatten ms) -- | creates a matrix from a horizontal list of matrices joinHoriz :: Element t => [Matrix t] -> Matrix t joinHoriz ms = trans. joinVert . map trans $ ms -{- | Creates a matrix from blocks given as a list of lists of matrices:+{- | Creates a matrix from blocks given as a list of lists of matrices. -@\> let a = 'diag' $ 'fromList' [5,7,2]-\> let b = 'reshape' 4 $ 'constant' (-1) 12-\> fromBlocks [[a,b],[b,a]]-(6><7)- [ 5.0, 0.0, 0.0, -1.0, -1.0, -1.0, -1.0- , 0.0, 7.0, 0.0, -1.0, -1.0, -1.0, -1.0- , 0.0, 0.0, 2.0, -1.0, -1.0, -1.0, -1.0- , -1.0, -1.0, -1.0, -1.0, 5.0, 0.0, 0.0- , -1.0, -1.0, -1.0, -1.0, 0.0, 7.0, 0.0- , -1.0, -1.0, -1.0, -1.0, 0.0, 0.0, 2.0 ]@+Single row/column components are automatically expanded to match the+corresponding common row and column:++@\> let disp = putStr . dispf 2+\> let vector xs = fromList xs :: Vector Double+\> let diagl = diag . vector+\> let rowm = asRow . vector++\> disp $ fromBlocks [[ident 5, 7, rowm[10,20]], [3, diagl[1,2,3], 0]]++8x10+1 0 0 0 0 7 7 7 10 20+0 1 0 0 0 7 7 7 10 20+0 0 1 0 0 7 7 7 10 20+0 0 0 1 0 7 7 7 10 20+0 0 0 0 1 7 7 7 10 20+3 3 3 3 3 1 0 0 0 0+3 3 3 3 3 0 2 0 0 0+3 3 3 3 3 0 0 3 0 0@ -} fromBlocks :: Element t => [[Matrix t]] -> Matrix t-fromBlocks = joinVert . map joinHoriz +fromBlocks = fromBlocksRaw . adaptBlocks +fromBlocksRaw mms = joinVert . map joinHoriz $ mms++adaptBlocks ms = ms' where+ bc = case common length ms of+ Just c -> c+ Nothing -> error "fromBlocks requires rectangular [[Matrix]]"+ rs = map (compatdim . map rows) ms+ cs = map (compatdim . map cols) (transpose ms)+ szs = sequence [rs,cs]+ ms' = partit bc $ zipWith g szs (concat ms)++ g [Just nr,Just nc] m+ | nr == r && nc == c = m+ | r == 1 && c == 1 = reshape nc (constant x (nr*nc))+ | r == 1 = fromRows (replicate nr (flatten m))+ | otherwise = fromColumns (replicate nc (flatten m))+ where+ r = rows m+ c = cols m+ x = m@@>(0,0)+ g _ _ = error "inconsistent dimensions in fromBlocks"++-----------------------------------------------------------+ -- | Reverse rows flipud :: Element t => Matrix t -> Matrix t flipud m = fromRows . reverse . toRows $ m@@ -229,7 +262,7 @@ @import Text.Printf(printf)@ -@disp = putStrLn . format \" \" (printf \"%.2f\")@+@disp = putStr . format \" \" (printf \"%.2f\")@ -} format :: (Element t) => String -> (t -> String) -> Matrix t -> String@@ -244,6 +277,71 @@ dispC :: Int -> Matrix (Complex Double) -> IO () dispC d m = disp m (shfc d) -}++-------------------------------------------------------------------+-- display utilities++{- | Show a matrix with \"autoscaling\" and a given number of decimal places.++@disp = putStr . disps 2++\> disp $ 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.++@disp = putStr . dispf 3++\> disp (1/3 + ident 4)+4x4+1.333 0.333 0.333 0.333+0.333 1.333 0.333 0.333+0.333 0.333 1.333 0.333+0.333 0.333 0.333 1.333+@+-}+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 where+ lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx+ where shx = show x++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 = 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.++@disp = putStr . vecdisp (dispf 2)++\> disp (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+++-------------------------------------------------------------------- -- | reads a matrix from a string containing a table of numbers. readMatrix :: String -> Matrix Double
lib/Data/Packed/Random.hs view
@@ -21,6 +21,7 @@ import Numeric.GSL.Vector import Data.Packed+import Numeric.LinearAlgebra.Linear import Numeric.LinearAlgebra.Algorithms import Numeric.LinearAlgebra.Instances() import Numeric.LinearAlgebra.Interface@@ -54,7 +55,7 @@ d = dim a dat = toRows $ reshape n $ randomVector seed Uniform (n*d) am = constant 1 n `outer` a- m = fromColumns (zipWith (.*) cs dat) + am+ m = fromColumns (zipWith scale cs dat) + am ------------ utilities -------------------------------
lib/Data/Packed/Vector.hs view
@@ -18,6 +18,7 @@ dim, (@>), subVector, join, constant, linspace,+ vecdisp, vectorMax, vectorMin, vectorMaxIndex, vectorMinIndex, mapVector, zipVector, fscanfVector, fprintfVector, freadVector, fwriteVector,@@ -32,6 +33,10 @@ @\> linspace 5 (-3,7) 5 |> [-3.0,-0.5,2.0,4.5,7.0]@++Logarithmic spacing can be defined as follows:++@logspace n (a,b) = 10 ** linspace n (a,b)@ -} linspace :: Int -> (Double, Double) -> Vector Double linspace n (a,b) = add a $ scale s $ fromList [0 .. fromIntegral n-1]@@ -70,3 +75,19 @@ buildVector :: Element a => Int -> (Int -> a) -> Vector a buildVector len f = fromList $ map f [0 .. (len - 1)]+++{- | Show a vector using a function for showing matrices.++@disp = putStr . vecdisp ('dispf' 2)++\> disp ('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
lib/Numeric/GSL/Minimization.hs view
@@ -51,8 +51,9 @@ ----------------------------------------------------------------------------- module Numeric.GSL.Minimization (- minimize, MinimizeMethod(..),- minimizeD, MinimizeMethodD(..),+ minimize, minimizeV, MinimizeMethod(..),+ minimizeD, minimizeVD, MinimizeMethodD(..),+ minimizeNMSimplex, minimizeConjugateGradient, minimizeVectorBFGS2@@ -82,7 +83,7 @@ | NMSimplex2 deriving (Enum,Eq,Show,Bounded) --- | Minimization without derivatives.+-- | Minimization without derivatives minimize :: MinimizeMethod -> Double -- ^ desired precision of the solution (size test) -> Int -- ^ maximum number of iterations allowed@@ -91,8 +92,40 @@ -> [Double] -- ^ starting point -> ([Double], Matrix Double) -- ^ solution vector and optimization path -minimize method = minimizeGen (fi (fromEnum method))+-- | Minimization without derivatives (vector version)+minimizeV :: MinimizeMethod+ -> Double -- ^ desired precision of the solution (size test)+ -> Int -- ^ maximum number of iterations allowed+ -> Vector Double -- ^ sizes of the initial search box+ -> (Vector Double -> Double) -- ^ function to minimize+ -> Vector Double -- ^ starting point+ -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path +minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi)+ where v2l (v,m) = (toList v, m)++ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2++minimizeV method eps maxit szv f xiv = unsafePerformIO $ do+ let n = dim xiv+ fp <- mkVecfun (iv f)+ rawpath <- ww2 withVector xiv withVector szv $ \xiv' szv' ->+ createMIO maxit (n+3)+ (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')+ "minimize"+ let it = round (rawpath @@> (maxit-1,0))+ path = takeRows it rawpath+ sol = cdat $ dropColumns 3 $ dropRows (it-1) path+ freeHaskellFunPtr fp+ return (sol, path)+++foreign import ccall "gsl-aux.h minimize"+ c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TVVM++----------------------------------------------------------------------------------++ data MinimizeMethodD = ConjugateFR | ConjugatePR | VectorBFGS@@ -111,49 +144,35 @@ -> [Double] -- ^ starting point -> ([Double], Matrix Double) -- ^ solution vector and optimization path -minimizeD method = minimizeDGen (fi (fromEnum method))-----------------------------------------------------------------------------ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2--minimizeGen method eps maxit sz f xi = unsafePerformIO $ do- let xiv = fromList xi- szv = fromList sz- n = dim xiv- fp <- mkVecfun (iv (f.toList))- rawpath <- ww2 withVector xiv withVector szv $ \xiv' szv' ->- createMIO maxit (n+3)- (c_minimize method fp eps (fi maxit) // xiv' // szv')- "minimize"- let it = round (rawpath @@> (maxit-1,0))- path = takeRows it rawpath- [sol] = toLists $ dropRows (it-1) path- freeHaskellFunPtr fp- return (drop 3 sol, path)---foreign import ccall "gsl-aux.h minimize"- c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TVVM------------------------------------------------------------------------------------+-- | Minimization with derivatives (vector version)+minimizeVD :: MinimizeMethodD+ -> Double -- ^ desired precision of the solution (gradient test)+ -> Int -- ^ maximum number of iterations allowed+ -> Double -- ^ size of the first trial step+ -> Double -- ^ tol (precise meaning depends on method)+ -> (Vector Double -> Double) -- ^ function to minimize+ -> (Vector Double -> Vector Double) -- ^ gradient+ -> Vector Double -- ^ starting point+ -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path +minimizeD method eps maxit istep tol f df xi = v2l $ minimizeVD+ method eps maxit istep tol (f.toList) (fromList.df.toList) (fromList xi)+ where v2l (v,m) = (toList v, m) -minimizeDGen method eps maxit istep tol f df xi = unsafePerformIO $ do- let xiv = fromList xi- n = dim xiv- f' = f . toList- df' = (checkdim1 n .fromList . df . toList)+minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do+ let n = dim xiv+ f' = f+ df' = (checkdim1 n . df) fp <- mkVecfun (iv f') dfp <- mkVecVecfun (aux_vTov df') rawpath <- withVector xiv $ \xiv' -> createMIO maxit (n+2)- (c_minimizeD method fp dfp istep tol eps (fi maxit) // xiv')+ (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv') "minimizeD" let it = round (rawpath @@> (maxit-1,0)) path = takeRows it rawpath- sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path+ sol = cdat $ dropColumns 2 $ dropRows (it-1) path freeHaskellFunPtr fp freeHaskellFunPtr dfp return (sol,path)
lib/Numeric/LinearAlgebra/Algorithms.hs view
@@ -147,6 +147,14 @@ -------------------------------------------------------------- +square m = rows m == cols m++vertical m = rows m >= cols m++exactHermitian m = m `equal` ctrans m++--------------------------------------------------------------+ -- | Full singular value decomposition. svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t) svd = svd'@@ -179,7 +187,6 @@ u' = takeColumns d u v' = takeColumns d v -vertical m = rows m >= cols m -- | Singular values and all right singular vectors. rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)@@ -255,12 +262,12 @@ -- -- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@ eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)-eigSH m | m `equal` ctrans m = eigSH' m+eigSH m | exactHermitian m = eigSH' m | otherwise = error "eigSH requires complex hermitian or real symmetric matrix" -- | Eigenvalues of a complex hermitian or real symmetric matrix. eigenvaluesSH :: Field t => Matrix t -> Vector Double-eigenvaluesSH m | m `equal` ctrans m = eigenvaluesSH' m+eigenvaluesSH m | exactHermitian m = eigenvaluesSH' m | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix" --------------------------------------------------------------@@ -317,14 +324,12 @@ -- -- If @c = chol m@ then @m == ctrans c \<> c@. chol :: Field t => Matrix t -> Matrix t-chol m | m `equal` ctrans m = cholSH m+chol m | exactHermitian m = cholSH m | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" -square m = rows m == cols m- -- | Determinant of a square matrix. det :: Field t => Matrix t -> t det m | square m = s * (product $ toList $ takeDiag $ lup)@@ -569,7 +574,7 @@ then Just (l,v) else Nothing where n = rows m- (l,v) = if m `equal` ctrans m+ (l,v) = if exactHermitian m then let (l',v') = eigSH m in (real l', v') else eig m
lib/Numeric/LinearAlgebra/Instances.hs view
@@ -80,8 +80,8 @@ || rows m2 == 1 && cols m2 == 1 || rows m1 == rows m2 && cols m1 == cols m2 -instance (Eq a, Element a) => Eq (Vector a) where- a == b = dim a == dim b && toList a == toList b+instance Linear Vector a => Eq (Vector a) where+ (==) = equal instance Num (Vector Double) where (+) = adaptScalar addConstant add (flip addConstant)@@ -99,8 +99,8 @@ abs = vectorMapC Abs fromInteger = fromList . return . fromInteger -instance (Eq a, Element a) => Eq (Matrix a) where- a == b = cols a == cols b && flatten a == flatten b+instance Linear Matrix a => Eq (Matrix a) where+ (==) = equal instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where (+) = liftMatrix2' (+)
lib/Numeric/LinearAlgebra/Interface.hs view
@@ -49,12 +49,16 @@ ---------------------------------------------------- +{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}+ -- | @x .* a = scale x a@ (.*) :: (Linear c a) => a -> c a -> c a infixl 7 .* a .* x = scale a x ----------------------------------------------------++{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-} -- | @a *\/ x = scale (recip x) a@ (*/) :: (Linear c a) => c a -> a -> c a
lib/Numeric/LinearAlgebra/Linear.hs view
@@ -23,7 +23,13 @@ -- | A generic interface for vectors and matrices to a few element-by-element functions in Numeric.GSL.Vector. class (Container c e) => Linear c e where+ -- | create a structure with a single element+ scalar :: e -> c e scale :: e -> c e -> c e+ -- | scale the element by element reciprocal of the object:+ --+ -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@+ scaleRecip :: e -> c e -> c e addConstant :: e -> c e -> c e add :: c e -> c e -> c e sub :: c e -> c e -> c e@@ -31,11 +37,9 @@ mul :: c e -> c e -> c e -- | element by element division divide :: c e -> c e -> c e- -- | scale the element by element reciprocal of the object: @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@- scaleRecip :: e -> c e -> c e equal :: c e -> c e -> Bool--- numequal :: Double -> c e -> c e -> Bool + instance Linear Vector Double where scale = vectorMapValR Scale scaleRecip = vectorMapValR Recip@@ -45,6 +49,7 @@ mul = vectorZipR Mul divide = vectorZipR Div equal u v = dim u == dim v && vectorMax (vectorMapR Abs (sub u v)) == 0.0+ scalar x = fromList [x] instance Linear Vector (Complex Double) where scale = vectorMapValC Scale@@ -55,6 +60,7 @@ mul = vectorZipC Mul divide = vectorZipC Div equal u v = dim u == dim v && vectorMax (mapVector magnitude (sub u v)) == 0.0+ scalar x = fromList [x] instance (Linear Vector a, Container Matrix a) => (Linear Matrix a) where scale x = liftMatrix (scale x)@@ -65,3 +71,4 @@ mul = liftMatrix2 mul divide = liftMatrix2 divide equal a b = cols a == cols b && flatten a `equal` flatten b+ scalar x = (1><1) [x]
lib/Numeric/LinearAlgebra/Tests.hs view
@@ -262,8 +262,8 @@ , utest "expm1" (expmTest1) , utest "expm2" (expmTest2) , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)- , utest "arith2" $ (((1+i) .* ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( (140*i-51).*1 :: CM)- , utest "arith3" $ exp (i.*ones(10,10)*pi) + 1 |~| 0+ , utest "arith2" $ ((scalar (1+i) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*i-51) :: CM)+ , utest "arith3" $ exp (scalar i * ones(10,10)*pi) + 1 |~| 0 , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3] , utest "gamma" (gamma 5 == 24.0) , besselTest@@ -279,10 +279,11 @@ && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0) , utest "rank" $ rank ((2><3)[1,0,0,1,6*eps,0]) == 1 && rank ((2><3)[1,0,0,1,7*eps,0]) == 2+ , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM) ] return () -makeUnitary v | realPart n > 1 = v */ n+makeUnitary v | realPart n > 1 = v / scalar n | otherwise = v where n = sqrt (conj v <.> v) @@ -321,6 +322,8 @@ where c = cos a s = sin a +multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Double))+ -------------------------------- multBench = do@@ -328,6 +331,15 @@ let b = ident 2000 :: Matrix Double a `seq` b `seq` putStrLn "" time "product of 1M different 3x3 matrices" (manymult (10^6))+ putStrLn ""+ time "product of 1M constant 1x1 matrices" (multb 1)+ time "product of 1M constant 3x3 matrices" (multb 3)+ --time "product of 1M constant 5x5 matrices" (multb 5)+ time "product of 1M const. 10x10 matrices" (multb 10)+ --time "product of 1M const. 15x15 matrices" (multb 15)+ time "product of 1M const. 20x20 matrices" (multb 20)+ --time "product of 1M const. 25x25 matrices" (multb 25)+ putStrLn "" time "product (1000 x 1000)<>(1000 x 1000)" (a<>a) time "product (2000 x 2000)<>(2000 x 2000)" (b<>b)
lib/Numeric/LinearAlgebra/Tests/Instances.hs view
@@ -180,7 +180,7 @@ l <- replicateM n (choose (0,100)) let s = diag (fromList l) p = v <> real s <> ctrans v- return $ PosDef (0.5 .* p + 0.5 .* ctrans p)+ return $ PosDef (0.5 * p + 0.5 * ctrans p) #if MIN_VERSION_QuickCheck(2,0,0) #else