diff --git a/CHANGES b/CHANGES --- a/CHANGES +++ b/CHANGES @@ -1,3 +1,21 @@ +0.10.0.0 +======== + +- Module reorganization + +- Support for Float and Complex Float elements (excluding LAPACK computations) + +- Binary instances for Vector and Matrix + +- optimiseMult + +- mapVectorM, mapVectorWithIndexM, unzipVectorWith, and related functions. + +- diagRect admits diagonal vectors of any length without producing an error, + and takes an additional argument for the off-diagonal elements. + +- different signatures in some functions + 0.9.3.0 ======= diff --git a/Setup.lhs b/Setup.lhs --- a/Setup.lhs +++ b/Setup.lhs @@ -1,7 +1,7 @@ #! /usr/bin/env runhaskell > import Distribution.Simple -> import System(system) +> import System.Process(system) > main = defaultMainWithHooks autoconfUserHooks {runTests = t} diff --git a/THANKS b/THANKS --- a/THANKS +++ b/THANKS @@ -1,6 +1,11 @@ I thank Don Stewart, Henning Thielemann, Bulat Ziganshin, Heinrich Apfelmus, and all the people in the Haskell mailing lists for their help. +I am particularly grateful to Vivian McPhail for his excellent +contributions: improved configure.hs, Binary instances for +Vector and Matrix, support for Float and Complex Float elements, +module reorganization, monadic mapVectorM, and many other improvements. + - Nico Mahlo discovered a bug in the eigendecomposition wrapper. - Frederik Eaton discovered a bug in the design of the wrappers. @@ -71,4 +76,6 @@ - Tim Sears reported the zgesdd problem also in intel mac. - Max Suica simplified the installation on Windows and improved the instructions. + +- John Billings reported an incompatibility with QuickCheck>=2.1.1 diff --git a/configure.hs b/configure.hs --- a/configure.hs +++ b/configure.hs @@ -17,7 +17,10 @@ -} -import System +import System.Process +import System.Exit +import System.Environment +import System.Directory(createDirectoryIfMissing) import Data.List(isPrefixOf, intercalate) import Distribution.Simple.LocalBuildInfo import Distribution.Simple.Configure @@ -37,13 +40,14 @@ ] -- compile a simple program with symbols from GSL and LAPACK with the given libs -testprog buildInfo libs fmks = +testprog bInfo buildInfo libs fmks = "echo \"#include \nint main(){zgesvd_(); gsl_sf_gamma(5);}\"" - ++" > /tmp/dummy.c; gcc " + ++" > " ++ (buildDir bInfo) ++ "/dummy.c; gcc " ++ (join $ ccOptions buildInfo) ++ " " ++ (join $ cppOptions buildInfo) ++ " " - ++ (join $ map ("-I"++) $ includeDirs buildInfo) - ++" /tmp/dummy.c -o /tmp/dummy " + ++ (join $ map ("-I"++) $ includeDirs buildInfo) ++ " " + ++ (buildDir bInfo) ++ "/dummy.c -o " + ++ (buildDir bInfo) ++ "/dummy " ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " " ++ (prepend "-l" $ libs) ++ " " ++ (prepend "-framework " fmks) ++ " > /dev/null 2> /dev/null" @@ -51,26 +55,28 @@ join = intercalate " " prepend x = unwords . map (x++) . words -check buildInfo libs fmks = (ExitSuccess ==) `fmap` system (testprog buildInfo libs fmks) +check bInfo buildInfo libs fmks = (ExitSuccess ==) `fmap` system (testprog bInfo buildInfo libs fmks) -- simple test for GSL -gsl buildInfo = "echo \"#include \nint main(){gsl_sf_gamma(5);}\"" - ++" > /tmp/dummy.c; gcc " +gsl bInfo buildInfo = "echo \"#include \nint main(){gsl_sf_gamma(5);}\"" + ++" > " ++ (buildDir bInfo) ++ "/dummy.c; gcc " ++ (join $ ccOptions buildInfo) ++ " " ++ (join $ cppOptions buildInfo) ++ " " - ++ (join $ map ("-I"++) $ includeDirs buildInfo) - ++ " /tmp/dummy.c -o /tmp/dummy " + ++ (join $ map ("-I"++) $ includeDirs buildInfo) ++ " " + ++ (buildDir bInfo) ++ "/dummy.c -o " + ++ (buildDir bInfo) ++ "/dummy " ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas" ++ " > /dev/null 2> /dev/null" -- test for gsl >= 1.12 -gsl112 buildInfo = +gsl112 bInfo buildInfo = "echo \"#include \nint main(){gsl_sf_exprel_n_CF_e(1,1,0);}\"" - ++" > /tmp/dummy.c; gcc /tmp/dummy.c " + ++" > " ++ (buildDir bInfo) ++ "/dummy.c; gcc " + ++ (buildDir bInfo) ++ "/dummy.c " ++ (join $ ccOptions buildInfo) ++ " " ++ (join $ cppOptions buildInfo) ++ " " ++ (join $ map ("-I"++) $ includeDirs buildInfo) - ++" -o /tmp/dummy " + ++" -o " ++ (buildDir bInfo) ++ "/dummy " ++ (join $ map ("-L"++) $ extraLibDirs buildInfo) ++ " -lgsl -lgslcblas" ++ " > /dev/null 2> /dev/null" @@ -78,11 +84,11 @@ checkCommand c = (ExitSuccess ==) `fmap` system c -- test different configurations until the first one works -try _ _ _ [] = return Nothing -try i b f (opt:rest) = do - ok <- check i (b ++ " " ++ opt) f +try _ _ _ _ [] = return Nothing +try l i b f (opt:rest) = do + ok <- check l i (b ++ " " ++ opt) f if ok then return (Just opt) - else try i b f rest + else try l i b f rest -- read --configure-option=link:lib1,lib2,lib3,etc linkop = "link:" @@ -110,11 +116,14 @@ let pref = if null (words (base ++ " " ++ auxpref)) then "gsl lapack" else auxpref fullOpts = map ((pref++" ")++) opts - r <- try buildInfo base fwks fullOpts + -- create the build directory (used for tmp files) if necessary + createDirectoryIfMissing True $ buildDir bInfo + + r <- try bInfo buildInfo base fwks fullOpts case r of Nothing -> do putStrLn " FAIL" - g <- checkCommand $ gsl buildInfo + g <- checkCommand $ gsl bInfo buildInfo if g then putStrLn " *** Sorry, I can't link LAPACK." else putStrLn " *** Sorry, I can't link GSL." @@ -124,7 +133,7 @@ writeFile "hmatrix.buildinfo" ("buildable: False\n") Just ops -> do putStrLn " OK" - g <- checkCommand $ gsl112 buildInfo + g <- checkCommand $ gsl112 bInfo buildInfo writeFile "hmatrix.buildinfo" $ "extra-libraries: " ++ ops ++ "\n" ++ if g diff --git a/examples/Real.hs b/examples/Real.hs --- a/examples/Real.hs +++ b/examples/Real.hs @@ -67,13 +67,13 @@ zeros :: Int -- ^ rows -> Int -- ^ columns -> Matrix Double -zeros r c = reshape c (constant 0 (r*c)) +zeros r c = konst 0 (r,c) -- | Create a matrix or ones. ones :: Int -- ^ rows -> Int -- ^ columns -> Matrix Double -ones r c = reshape c (constant 1 (r*c)) +ones r c = konst 1 (r,c) -- | Concatenation of real vectors. infixl 9 # diff --git a/examples/monadic.hs b/examples/monadic.hs new file mode 100644 --- /dev/null +++ b/examples/monadic.hs @@ -0,0 +1,118 @@ +-- monadic computations +-- (contributed by Vivian McPhail) + +import Numeric.LinearAlgebra +import Control.Monad.State.Strict +import Control.Monad.Maybe +import Foreign.Storable(Storable) +import System.Random(randomIO) + +------------------------------------------- + +-- an instance of MonadIO, a monad transformer +type VectorMonadT = StateT Int IO + +test1 :: Vector Int -> IO (Vector Int) +test1 = mapVectorM $ \x -> do + putStr $ (show x) ++ " " + return (x + 1) + +-- we can have an arbitrary monad AND do IO +addInitialM :: Vector Int -> VectorMonadT () +addInitialM = mapVectorM_ $ \x -> do + i <- get + liftIO $ putStr $ (show $ x + i) ++ " " + put $ x + i + +-- sum the values of the even indiced elements +sumEvens :: Vector Int -> Int +sumEvens = foldVectorWithIndex (\x a b -> if x `mod` 2 == 0 then a + b else b) 0 + +-- sum and print running total of evens +sumEvensAndPrint :: Vector Int -> VectorMonadT () +sumEvensAndPrint = mapVectorWithIndexM_ $ \ i x -> do + when (i `mod` 2 == 0) $ do + v <- get + put $ v + x + v' <- get + liftIO $ putStr $ (show v') ++ " " + + +indexPlusSum :: Vector Int -> VectorMonadT () +indexPlusSum v' = do + let f i x = do + s <- get + let inc = x+s + liftIO $ putStr $ show (i,inc) ++ " " + put inc + return inc + v <- mapVectorWithIndexM f v' + liftIO $ do + putStrLn "" + putStrLn $ show v + +------------------------------------------- + +-- short circuit +monoStep :: Double -> MaybeT (State Double) () +monoStep d = do + dp <- get + when (d < dp) (fail "negative difference") + put d +{-# INLINE monoStep #-} + +isMonotoneIncreasing :: Vector Double -> Bool +isMonotoneIncreasing v = + let res = evalState (runMaybeT $ (mapVectorM_ monoStep v)) (v @> 0) + in case res of + Nothing -> False + Just _ -> True + + +------------------------------------------- + +-- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs +successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool +successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (dim v - 1) v))) (v @> 0) + where step e = do + ep <- lift $ get + if t e ep + then lift $ put e + else (fail "successive_ test failed") + +-- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input +successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b +successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0) + where step e = do + ep <- get + put e + return $ f ep e + +------------------------------------------- + +v :: Vector Int +v = 10 |> [0..] + +w = fromList ([1..10]++[10,9..1]) :: Vector Double + + +main = do + v' <- test1 v + putStrLn "" + putStrLn $ show v' + evalStateT (addInitialM v) 0 + putStrLn "" + putStrLn $ show (sumEvens v) + evalStateT (sumEvensAndPrint v) 0 + putStrLn "" + evalStateT (indexPlusSum v) 0 + putStrLn "-----------------------" + mapVectorM_ print v + print =<< (mapVectorM (const randomIO) v :: IO (Vector Double)) + print =<< (mapVectorM (\a -> fmap (+a) randomIO) (5|>[0,100..1000]) :: IO (Vector Double)) + putStrLn "-----------------------" + print $ isMonotoneIncreasing w + print $ isMonotoneIncreasing (subVector 0 7 w) + print $ successive_ (>) v + print $ successive_ (>) w + print $ successive (+) v diff --git a/examples/parallel.hs b/examples/parallel.hs --- a/examples/parallel.hs +++ b/examples/parallel.hs @@ -1,6 +1,6 @@ -- $ runhaskell parallel.hs 2000 -import System(getArgs) +import System.Environment(getArgs) import Numeric.LinearAlgebra import Control.Parallel.Strategies import System.Time @@ -15,10 +15,10 @@ main = do n <- (read . head) `fmap` getArgs let m = ident n :: Matrix Double - 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 + time $ print $ maxElement $ takeDiag $ m <> m + time $ print $ maxElement $ takeDiag $ parMul 2 m m + time $ print $ maxElement $ takeDiag $ parMul 4 m m + time $ print $ maxElement $ takeDiag $ parMul 8 m m time act = do t0 <- getClockTime diff --git a/examples/pca1.hs b/examples/pca1.hs --- a/examples/pca1.hs +++ b/examples/pca1.hs @@ -2,7 +2,7 @@ import Numeric.LinearAlgebra import System.Directory(doesFileExist) -import System(system) +import System.Process(system) import Control.Monad(when) type Vec = Vector Double diff --git a/examples/pca2.hs b/examples/pca2.hs --- a/examples/pca2.hs +++ b/examples/pca2.hs @@ -3,7 +3,7 @@ import Numeric.LinearAlgebra import Graphics.Plot import System.Directory(doesFileExist) -import System(system) +import System.Process(system) import Control.Monad(when) type Vec = Vector Double diff --git a/examples/vector.hs b/examples/vector.hs --- a/examples/vector.hs +++ b/examples/vector.hs @@ -14,7 +14,7 @@ fromVector v = unsafeFromForeignPtr p i n where (p,i,n) = V.unsafeToForeignPtr v -toVector :: H.Vector t -> V.Vector t +toVector :: Storable t => H.Vector t -> V.Vector t toVector v = V.unsafeFromForeignPtr p i n where (p,i,n) = unsafeToForeignPtr v @@ -22,11 +22,10 @@ v = V.slice 5 10 (V.fromList [1 .. 10::Double] V.++ V.replicate 10 7) -w = subVector 2 3 (linspace 10 (0,2)) +w = subVector 2 3 (linspace 5 (0,1)) :: Vector Double main = do print v print $ fromVector v print w print $ toVector w - diff --git a/hmatrix.cabal b/hmatrix.cabal --- a/hmatrix.cabal +++ b/hmatrix.cabal @@ -1,5 +1,5 @@ Name: hmatrix -Version: 0.9.3.0 +Version: 0.10.0.0 License: GPL License-file: LICENSE Author: Alberto Ruiz @@ -11,11 +11,19 @@ and other numerical computations, internally implemented using GSL, BLAS and LAPACK. . - See also hmatrix-special and hmatrix-glpk. + The Linear Algebra API is organized as follows: + . + - "Data.Packed": structure manipulation + . + - "Numeric.Container": simple numeric functions + . + - "Numeric.LinearAlgebra.Algorithms": matrix computations + . + - "Numeric.LinearAlgebra": everything + instances of standard Haskell numeric classes Category: Math tested-with: GHC ==6.10.4, GHC ==6.12.1 -cabal-version: >=1.2 +cabal-version: >=1.6 build-type: Custom @@ -45,6 +53,7 @@ examples/devel/ej2/functions.c examples/Real.hs examples/vector.hs + examples/monadic.hs extra-source-files: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h, lib/Numeric/LinearAlgebra/LAPACK/clapack.h @@ -74,7 +83,8 @@ Build-Depends: base >= 4 && < 5, array, storable-complex, - process + process, + binary Extensions: ForeignFunctionInterface, CPP @@ -91,35 +101,38 @@ Numeric.GSL.Root, Numeric.GSL.Fitting, Numeric.GSL.ODE, - Numeric.GSL.Vector, Numeric.GSL, + Numeric.Container, Numeric.LinearAlgebra, Numeric.LinearAlgebra.LAPACK, - Numeric.LinearAlgebra.Interface, Numeric.LinearAlgebra.Algorithms, Graphics.Plot, - -- Data.Packed.Convert, Data.Packed.ST, - Data.Packed.Development, - Data.Packed.Random + Data.Packed.Development other-modules: Data.Packed.Internal, Data.Packed.Internal.Common, Data.Packed.Internal.Signatures, Data.Packed.Internal.Vector, Data.Packed.Internal.Matrix, - Numeric.LinearAlgebra.Linear, - Numeric.LinearAlgebra.Instances, - Numeric.GSL.Internal + Data.Packed.Random, + Numeric.GSL.Internal, + Numeric.GSL.Vector, + Numeric.Conversion, + Numeric.ContainerBoot, + Numeric.IO, + Numeric.Chain, + Numeric.Vector, + Numeric.Matrix C-sources: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c, lib/Numeric/GSL/gsl-aux.c if flag(vector) - Build-Depends: vector + Build-Depends: vector >= 0.7 cpp-options: -DVECTOR if flag(tests) - Build-Depends: QuickCheck, HUnit + Build-Depends: QuickCheck, HUnit, random exposed-modules: Numeric.LinearAlgebra.Tests other-modules: Numeric.LinearAlgebra.Tests.Instances, Numeric.LinearAlgebra.Tests.Properties @@ -160,7 +173,6 @@ extra-libraries: extra-lib-dirs: - source-repository head - type: darcs - location: http://code.haskell.org/hmatrix - +source-repository head + type: darcs + location: http://code.haskell.org/hmatrix diff --git a/lib/Data/Packed.hs b/lib/Data/Packed.hs --- a/lib/Data/Packed.hs +++ b/lib/Data/Packed.hs @@ -1,14 +1,14 @@ ----------------------------------------------------------------------------- {- | Module : Data.Packed -Copyright : (c) Alberto Ruiz 2006-7 +Copyright : (c) Alberto Ruiz 2006-2010 License : GPL-style Maintainer : Alberto Ruiz (aruiz at um dot es) Stability : provisional Portability : uses ffi -The Vector and Matrix types and some utilities. +Types for dense 'Vector' and 'Matrix' of 'Storable' elements. -} ----------------------------------------------------------------------------- @@ -16,12 +16,13 @@ module Data.Packed ( module Data.Packed.Vector, module Data.Packed.Matrix, - module Data.Packed.Random, - module Data.Complex +-- module Numeric.Conversion, +-- module Data.Packed.Random, +-- module Data.Complex ) where import Data.Packed.Vector import Data.Packed.Matrix -import Data.Packed.Random -import Data.Complex - +--import Data.Packed.Random +--import Data.Complex +--import Numeric.Conversion diff --git a/lib/Data/Packed/Development.hs b/lib/Data/Packed/Development.hs --- a/lib/Data/Packed/Development.hs +++ b/lib/Data/Packed/Development.hs @@ -17,7 +17,6 @@ module Data.Packed.Development ( createVector, createMatrix, - Adapt, vec, mat, app1, app2, app3, app4, app5, app6, app7, app8, app9, app10, diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs --- a/lib/Data/Packed/Internal/Common.hs +++ b/lib/Data/Packed/Internal/Common.hs @@ -82,41 +82,27 @@ type Adapt f t r = t -> ((f -> r) -> IO()) -> IO() -app1 :: f - -> Adapt f t (IO CInt) - -> t - -> String - -> IO() - -app2 :: f - -> Adapt f t1 r - -> t1 - -> Adapt r t2 (IO CInt) - -> t2 - -> String - -> IO() - -app3 :: f - -> Adapt f t1 r1 - -> t1 - -> Adapt r1 t2 r2 - -> t2 - -> Adapt r2 t3 (IO CInt) - -> t3 - -> String - -> IO() +type Adapt1 f t1 = Adapt f t1 (IO CInt) -> t1 -> String -> IO() +type Adapt2 f t1 r1 t2 = Adapt f t1 r1 -> t1 -> Adapt1 r1 t2 +type Adapt3 f t1 r1 t2 r2 t3 = Adapt f t1 r1 -> t1 -> Adapt2 r1 t2 r2 t3 +type Adapt4 f t1 r1 t2 r2 t3 r3 t4 = Adapt f t1 r1 -> t1 -> Adapt3 r1 t2 r2 t3 r3 t4 +type Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5 = Adapt f t1 r1 -> t1 -> Adapt4 r1 t2 r2 t3 r3 t4 r4 t5 +type Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 = Adapt f t1 r1 -> t1 -> Adapt5 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 +type Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 = Adapt f t1 r1 -> t1 -> Adapt6 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 +type Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 = Adapt f t1 r1 -> t1 -> Adapt7 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 +type Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 = Adapt f t1 r1 -> t1 -> Adapt8 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 +type Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 = Adapt f t1 r1 -> t1 -> Adapt9 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 -app4 :: f - -> Adapt f t1 r1 - -> t1 - -> Adapt r1 t2 r2 - -> t2 - -> Adapt r2 t3 r3 - -> t3 - -> Adapt r3 t4 (IO CInt) - -> t4 - -> String - -> IO() +app1 :: f -> Adapt1 f t1 +app2 :: f -> Adapt2 f t1 r1 t2 +app3 :: f -> Adapt3 f t1 r1 t2 r2 t3 +app4 :: f -> Adapt4 f t1 r1 t2 r2 t3 r3 t4 +app5 :: f -> Adapt5 f t1 r1 t2 r2 t3 r3 t4 r4 t5 +app6 :: f -> Adapt6 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 +app7 :: f -> Adapt7 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 +app8 :: f -> Adapt8 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 +app9 :: f -> Adapt9 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 +app10 :: f -> Adapt10 f t1 r1 t2 r2 t3 r3 t4 r4 t5 r5 t6 r6 t7 r7 t8 r8 t9 r9 t10 app1 f w1 o1 s = w1 o1 $ \a1 -> f // a1 // check s app2 f w1 o1 w2 o2 s = ww2 w1 o1 w2 o2 $ \a1 a2 -> f // a1 // a2 // check s diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs @@ -29,7 +29,6 @@ liftMatrix, liftMatrix2, (@@>), saveMatrix, - fromComplexV, toComplexV, conjV, singleton ) where @@ -76,7 +75,11 @@ data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq) --- | Matrix representation suitable for GSL and LAPACK computations. +{- | Matrix representation suitable for GSL and LAPACK computations. + +The elements are stored in a continuous memory array. + +-} data Matrix t = MC { irows :: {-# UNPACK #-} !Int , icols :: {-# UNPACK #-} !Int , cdat :: {-# UNPACK #-} !(Vector t) } @@ -222,13 +225,13 @@ , 9.0, 10.0, 11.0, 12.0 ]@ -} -reshape :: Element t => Int -> Vector t -> Matrix t +reshape :: Storable t => Int -> Vector t -> Matrix t reshape c v = matrixFromVector RowMajor c v singleton x = reshape 1 (fromList [x]) -- | application of a vector function on the flattened matrix elements -liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b +liftMatrix :: (Storable a, Storable b) => (Vector a -> Vector b) -> Matrix a -> Matrix b liftMatrix f MC { icols = c, cdat = d } = matrixFromVector RowMajor c (f d) liftMatrix f MF { icols = c, fdat = d } = matrixFromVector ColumnMajor c (f d) @@ -246,21 +249,37 @@ ------------------------------------------------------------------ --- | Auxiliary class. -class (Storable a, Floating a) => Element a where +{- | Supported matrix elements. + + This class provides optimized internal + operations for selected element types. + It provides unoptimised defaults for any 'Storable' type, + so you can create instances simply as: + @instance Element Foo@. +-} +class (Storable a) => Element a where subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix -> Matrix a -> Matrix a subMatrixD = subMatrix' transdata :: Int -> Vector a -> Int -> Vector a - transdata = transdata' + transdata = transdataP -- transdata' constantD :: a -> Int -> Vector a - constantD = constant' + constantD = constantP -- constant' + +instance Element Float where + transdata = transdataAux ctransF + constantD = constantAux cconstantF + instance Element Double where transdata = transdataAux ctransR constantD = constantAux cconstantR +instance Element (Complex Float) where + transdata = transdataAux ctransQ + constantD = constantAux cconstantQ + instance Element (Complex Double) where transdata = transdataAux ctransC constantD = constantAux cconstantC @@ -308,8 +327,27 @@ r2 = dim d `div` c2 noneed = r1 == 1 || c1 == 1 +transdataP :: Storable a => Int -> Vector a -> Int -> Vector a +transdataP c1 d c2 = + if noneed + then d + else unsafePerformIO $ do + v <- createVector (dim d) + unsafeWith d $ \pd -> + unsafeWith v $ \pv -> + ctransP (fi r1) (fi c1) (castPtr pd) (fi sz) (fi r2) (fi c2) (castPtr pv) (fi sz) // check "transdataP" + return v + where r1 = dim d `div` c1 + r2 = dim d `div` c2 + sz = sizeOf (d @> 0) + noneed = r1 == 1 || c1 == 1 + +foreign import ccall "transF" ctransF :: TFMFM foreign import ccall "transR" ctransR :: TMM +foreign import ccall "transQ" ctransQ :: TQMQM foreign import ccall "transC" ctransC :: TCMCM +foreign import ccall "transP" ctransP :: CInt -> CInt -> Ptr () -> CInt -> CInt -> CInt -> Ptr () -> CInt -> IO CInt + ---------------------------------------------------------------------- constant' v n = unsafePerformIO $ do @@ -329,13 +367,33 @@ free px return v +constantF :: Float -> Int -> Vector Float +constantF = constantAux cconstantF +foreign import ccall "constantF" cconstantF :: Ptr Float -> TF + constantR :: Double -> Int -> Vector Double constantR = constantAux cconstantR foreign import ccall "constantR" cconstantR :: Ptr Double -> TV +constantQ :: Complex Float -> Int -> Vector (Complex Float) +constantQ = constantAux cconstantQ +foreign import ccall "constantQ" cconstantQ :: Ptr (Complex Float) -> TQV + constantC :: Complex Double -> Int -> Vector (Complex Double) constantC = constantAux cconstantC foreign import ccall "constantC" cconstantC :: Ptr (Complex Double) -> TCV + +constantP :: Storable a => a -> Int -> Vector a +constantP a n = unsafePerformIO $ do + let sz = sizeOf a + v <- createVector n + unsafeWith v $ \p -> do + alloca $ \k -> do + poke k a + cconstantP (castPtr k) (fi n) (castPtr p) (fi sz) // check "constantP" + return v +foreign import ccall "constantP" cconstantP :: Ptr () -> CInt -> Ptr () -> CInt -> IO CInt + ---------------------------------------------------------------------- -- | Extracts a submatrix from a matrix. @@ -364,21 +422,6 @@ subMatrix' (r0,c0) (rt,ct) (MC _r c v) = MC rt ct $ subMatrix'' (r0,c0) (rt,ct) c v subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m) - --------------------------------------------------------------------------- - --- | obtains the complex conjugate of a complex vector -conjV :: Vector (Complex Double) -> Vector (Complex Double) -conjV = mapVector conjugate - --- | creates a complex vector from vectors with real and imaginary parts -toComplexV :: (Vector Double, Vector Double) -> Vector (Complex Double) -toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i] - --- | the inverse of 'toComplex' -fromComplexV :: Vector (Complex Double) -> (Vector Double, Vector Double) -fromComplexV z = (r,i) where - [r,i] = toColumns $ reshape 2 $ asReal z -------------------------------------------------------------------------- diff --git a/lib/Data/Packed/Internal/Signatures.hs b/lib/Data/Packed/Internal/Signatures.hs --- a/lib/Data/Packed/Internal/Signatures.hs +++ b/lib/Data/Packed/Internal/Signatures.hs @@ -18,11 +18,21 @@ import Data.Complex import Foreign.C.Types +type PF = Ptr Float -- type PD = Ptr Double -- +type PQ = Ptr (Complex Float) -- type PC = Ptr (Complex Double) -- +type TF = CInt -> PF -> IO CInt -- +type TFF = CInt -> PF -> TF -- +type TFV = CInt -> PF -> TV -- +type TVF = CInt -> PD -> TF -- +type TFFF = CInt -> PF -> TFF -- type TV = CInt -> PD -> IO CInt -- type TVV = CInt -> PD -> TV -- type TVVV = CInt -> PD -> TVV -- +type TFM = CInt -> CInt -> PF -> IO CInt -- +type TFMFM = CInt -> CInt -> PF -> TFM -- +type TFMFMFM = CInt -> CInt -> PF -> TFMFM -- type TM = CInt -> CInt -> PD -> IO CInt -- type TMM = CInt -> CInt -> PD -> TM -- type TVMM = CInt -> PD -> TMM -- @@ -47,6 +57,14 @@ type TCV = CInt -> PC -> IO CInt -- type TCVCV = CInt -> PC -> TCV -- type TCVCVCV = CInt -> PC -> TCVCV -- +type TCVV = CInt -> PC -> TV -- +type TQV = CInt -> PQ -> IO CInt -- +type TQVQV = CInt -> PQ -> TQV -- +type TQVQVQV = CInt -> PQ -> TQVQV -- +type TQVF = CInt -> PQ -> TF -- +type TQM = CInt -> CInt -> PQ -> IO CInt -- +type TQMQM = CInt -> CInt -> PQ -> TQM -- +type TQMQMQM = CInt -> CInt -> PQ -> TQMQM -- type TCMCV = CInt -> CInt -> PC -> TCV -- type TVCV = CInt -> PD -> TCV -- type TCVM = CInt -> PC -> TM -- diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs --- a/lib/Data/Packed/Internal/Vector.hs +++ b/lib/Data/Packed/Internal/Vector.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns #-} +{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-} ----------------------------------------------------------------------------- -- | -- Module : Data.Packed.Internal.Vector @@ -17,10 +17,11 @@ Vector, dim, fromList, toList, (|>), join, (@>), safe, at, at', subVector, takesV, - mapVector, zipVector, - foldVector, foldVectorG, foldLoop, + mapVector, zipVectorWith, unzipVectorWith, + mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_, + foldVector, foldVectorG, foldLoop, foldVectorWithIndex, createVector, vec, - asComplex, asReal, + asComplex, asReal, float2DoubleV, double2FloatV, fwriteVector, freadVector, fprintfVector, fscanfVector, cloneVector, unsafeToForeignPtr, @@ -44,7 +45,7 @@ import GHC.Base #if __GLASGOW_HASKELL__ < 612 -import GHC.IOBase +import GHC.IOBase hiding (liftIO) #endif #ifdef VECTOR @@ -70,11 +71,11 @@ , fptr :: {-# UNPACK #-} !(ForeignPtr t) -- ^ foreign pointer to the memory block } -unsafeToForeignPtr :: Vector a -> (ForeignPtr a, Int, Int) +unsafeToForeignPtr :: Storable a => Vector a -> (ForeignPtr a, Int, Int) unsafeToForeignPtr v = (fptr v, ioff v, idim v) -- | Same convention as in Roman Leshchinskiy's vector package. -unsafeFromForeignPtr :: ForeignPtr a -> Int -> Int -> Vector a +unsafeFromForeignPtr :: Storable a => ForeignPtr a -> Int -> Int -> Vector a unsafeFromForeignPtr fp i n | n > 0 = V {ioff = i, idim = n, fptr = fp} | otherwise = error "unsafeFromForeignPtr with dim < 1" @@ -264,17 +265,33 @@ --------------------------------------------------------------- -- | transforms a complex vector into a real vector with alternating real and imaginary parts -asReal :: Vector (Complex Double) -> Vector Double ---asReal v = V { ioff = 2*ioff v, idim = 2*dim v, fptr = castForeignPtr (fptr v) } +asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n) where (fp,i,n) = unsafeToForeignPtr v -- | transforms a real vector into a complex vector with alternating real and imaginary parts -asComplex :: Vector Double -> Vector (Complex Double) ---asComplex v = V { ioff = ioff v `div` 2, idim = dim v `div` 2, fptr = castForeignPtr (fptr v) } +asComplex :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a) asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2) where (fp,i,n) = unsafeToForeignPtr v +--------------------------------------------------------------- + +float2DoubleV :: Vector Float -> Vector Double +float2DoubleV v = unsafePerformIO $ do + r <- createVector (dim v) + app2 c_float2double vec v vec r "float2double" + return r + +double2FloatV :: Vector Double -> Vector Float +double2FloatV v = unsafePerformIO $ do + r <- createVector (dim v) + app2 c_double2float vec v vec r "double2float2" + return r + + +foreign import ccall "float2double" c_float2double:: TFV +foreign import ccall "double2float" c_double2float:: TVF + ---------------------------------------------------------------- cloneVector :: Storable t => Vector t -> IO (Vector t) @@ -302,8 +319,8 @@ {-# INLINE mapVector #-} -- | zipWith for Vectors -zipVector :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c -zipVector f u v = unsafePerformIO $ do +zipVectorWith :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c +zipVectorWith f u v = unsafePerformIO $ do let n = min (dim u) (dim v) w <- createVector n unsafeWith u $ \pu -> @@ -316,16 +333,47 @@ go (k-1) go (n -1) return w -{-# INLINE zipVector #-} +{-# INLINE zipVectorWith #-} +-- | unzipWith for Vectors +unzipVectorWith :: (Storable (a,b), Storable c, Storable d) + => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d) +unzipVectorWith f u = unsafePerformIO $ do + let n = dim u + v <- createVector n + w <- createVector n + unsafeWith u $ \pu -> + unsafeWith v $ \pv -> + unsafeWith w $ \pw -> do + let go (-1) = return () + go !k = do z <- peekElemOff pu k + let (x,y) = f z + pokeElemOff pv k x + pokeElemOff pw k y + go (k-1) + go (n-1) + return (v,w) +{-# INLINE unzipVectorWith #-} + +foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b foldVector f x v = unsafePerformIO $ - unsafeWith (v::Vector Double) $ \p -> do + unsafeWith v $ \p -> do let go (-1) s = return s go !k !s = do y <- peekElemOff p k go (k-1::Int) (f y s) go (dim v -1) x {-# INLINE foldVector #-} +-- the zero-indexed index is passed to the folding function +foldVectorWithIndex :: Storable a => (Int -> a -> b -> b) -> b -> Vector a -> b +foldVectorWithIndex f x v = unsafePerformIO $ + unsafeWith v $ \p -> do + let go (-1) s = return s + go !k !s = do y <- peekElemOff p k + go (k-1::Int) (f k y s) + go (dim v -1) x +{-# INLINE foldVectorWithIndex #-} + foldLoop f s0 d = go (d - 1) s0 where go 0 s = f (0::Int) s @@ -338,6 +386,75 @@ ------------------------------------------------------------------- +-- | monadic map over Vectors +-- the monad @m@ must be strict +mapVectorM :: (Storable a, Storable b, Monad m) => (a -> m b) -> Vector a -> m (Vector b) +mapVectorM f v = do + w <- return $! unsafePerformIO $! createVector (dim v) + mapVectorM' w 0 (dim v -1) + return w + where mapVectorM' w' !k !t + | k == t = do + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + y <- f x + return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y + | otherwise = do + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + y <- f x + _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y + mapVectorM' w' (k+1) t +{-# INLINE mapVectorM #-} + +-- | monadic map over Vectors +mapVectorM_ :: (Storable a, Monad m) => (a -> m ()) -> Vector a -> m () +mapVectorM_ f v = do + mapVectorM' 0 (dim v -1) + where mapVectorM' !k !t + | k == t = do + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + f x + | otherwise = do + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + _ <- f x + mapVectorM' (k+1) t +{-# INLINE mapVectorM_ #-} + +-- | monadic map over Vectors with the zero-indexed index passed to the mapping function +-- the monad @m@ must be strict +mapVectorWithIndexM :: (Storable a, Storable b, Monad m) => (Int -> a -> m b) -> Vector a -> m (Vector b) +mapVectorWithIndexM f v = do + w <- return $! unsafePerformIO $! createVector (dim v) + mapVectorM' w 0 (dim v -1) + return w + where mapVectorM' w' !k !t + | k == t = do + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + y <- f k x + return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y + | otherwise = do + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + y <- f k x + _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y + mapVectorM' w' (k+1) t +{-# INLINE mapVectorWithIndexM #-} + +-- | monadic map over Vectors with the zero-indexed index passed to the mapping function +mapVectorWithIndexM_ :: (Storable a, Monad m) => (Int -> a -> m ()) -> Vector a -> m () +mapVectorWithIndexM_ f v = do + mapVectorM' 0 (dim v -1) + where mapVectorM' !k !t + | k == t = do + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + f k x + | otherwise = do + x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k + _ <- f k x + mapVectorM' (k+1) t +{-# INLINE mapVectorWithIndexM_ #-} + +------------------------------------------------------------------- + + -- | Loads a vector from an ASCII file (the number of elements must be known in advance). fscanfVector :: FilePath -> Int -> IO (Vector Double) fscanfVector filename n = do @@ -379,3 +496,4 @@ free charname foreign import ccall "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV + diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs --- a/lib/Data/Packed/Matrix.hs +++ b/lib/Data/Packed/Matrix.hs @@ -1,8 +1,13 @@ -{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} + ----------------------------------------------------------------------------- -- | -- Module : Data.Packed.Matrix --- Copyright : (c) Alberto Ruiz 2007 +-- Copyright : (c) Alberto Ruiz 2007-10 -- License : GPL-style -- -- Maintainer : Alberto Ruiz @@ -11,11 +16,14 @@ -- -- A Matrix representation suitable for numerical computations using LAPACK and GSL. -- +-- This module provides basic functions for manipulation of structure. + ----------------------------------------------------------------------------- module Data.Packed.Matrix ( - Element, Container(..), - Matrix,rows,cols, + Matrix, + Element, + rows,cols, (><), trans, reshape, flatten, @@ -28,22 +36,67 @@ flipud, fliprl, subMatrix, takeRows, dropRows, takeColumns, dropColumns, extractRows, - ident, diag, diagRect, takeDiag, - liftMatrix, liftMatrix2, liftMatrix2Auto, - dispf, disps, dispcf, latexFormat, format, - loadMatrix, saveMatrix, fromFile, fileDimensions, - readMatrix, fromArray2D + diagRect, takeDiag, + liftMatrix, liftMatrix2, liftMatrix2Auto,fromArray2D ) where import Data.Packed.Internal import qualified Data.Packed.ST as ST -import Data.Packed.Vector -import Data.Array -import System.Process(readProcess) -import Text.Printf(printf) import Data.List(transpose,intersperse) -import Data.Complex +import Data.Array + +import Data.Binary +import Foreign.Storable +import Control.Monad(replicateM) +--import Control.Arrow((***)) +--import GHC.Float(double2Float,float2Double) + + +------------------------------------------------------------------- + +instance (Binary a, Element a, Storable a) => Binary (Matrix a) where + put m = do + let r = rows m + let c = cols m + put r + put c + mapM_ (\i -> mapM_ (\j -> put $ m @@> (i,j)) [0..(c-1)]) [0..(r-1)] + get = do + r <- get + c <- get + xs <- replicateM r $ replicateM c get + return $ fromLists xs + +------------------------------------------------------------------- + +instance (Show a, Element a) => (Show (Matrix a)) where + show m = (sizes++) . dsp . map (map show) . toLists $ m + where sizes = "("++show (rows m)++"><"++show (cols m)++")\n" + +dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp + where + mt = transpose as + longs = map (maximum . map length) mt + mtp = zipWith (\a b -> map (pad a) b) longs mt + pad n str = replicate (n - length str) ' ' ++ str + unwords' = concat . intersperse ", " + +------------------------------------------------------------------ + +instance (Element a, Read a) => Read (Matrix a) where + readsPrec _ s = [((rs>' $ dims + + +breakAt c l = (a++[c],tail b) where + (a,b) = break (==c) l + +------------------------------------------------------------------ + -- | creates a matrix from a vertical list of matrices joinVert :: Element t => [Matrix t] -> Matrix t joinVert ms = case common cols ms of @@ -92,7 +145,7 @@ g [Just nr,Just nc] m | nr == r && nc == c = m - | r == 1 && c == 1 = reshape nc (constant x (nr*nc)) + | r == 1 && c == 1 = reshape nc (constantD x (nr*nc)) | r == 1 = fromRows (replicate nr (flatten m)) | otherwise = fromColumns (replicate nc (flatten m)) where @@ -113,28 +166,19 @@ ------------------------------------------------------------ --- | Creates a square matrix with a given diagonal. -diag :: Element a => Vector a -> Matrix a -diag v = ST.runSTMatrix $ do - let d = dim v - m <- ST.newMatrix 0 d d - mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1] - return m - -{- | creates a rectangular diagonal matrix +{- | creates a rectangular diagonal matrix: -@> diagRect (constant 5 3) 3 4 :: Matrix Double -(3><4) - [ 5.0, 0.0, 0.0, 0.0 - , 0.0, 5.0, 0.0, 0.0 - , 0.0, 0.0, 5.0, 0.0 ]@ +@> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double +(4><5) + [ 10.0, 7.0, 7.0, 7.0, 7.0 + , 7.0, 20.0, 7.0, 7.0, 7.0 + , 7.0, 7.0, 30.0, 7.0, 7.0 + , 7.0, 7.0, 7.0, 7.0, 7.0 ]@ -} -diagRect :: (Element t, Num t) => Vector t -> Int -> Int -> Matrix t -diagRect v r c - | dim v < min r c = error "diagRect called with dim v < min r c" - | otherwise = ST.runSTMatrix $ do - m <- ST.newMatrix 0 r c - let d = min r c +diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t +diagRect z v r c = ST.runSTMatrix $ do + m <- ST.newMatrix z r c + let d = min r c `min` (dim v) mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1] return m @@ -142,10 +186,6 @@ takeDiag :: (Element t) => Matrix t -> Vector t takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (cols m) -1]] --- | creates the identity matrix of given dimension -ident :: Element a => Int -> Matrix a -ident n = diag (constant 1 n) - ------------------------------------------------------------ {- | An easy way to create a matrix: @@ -169,7 +209,7 @@ , 4.0, 5.0, 6.0 ]@ -} -(><) :: (Element a) => Int -> Int -> [a] -> Matrix a +(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a r >< c = f where f l | dim v == r*c = matrixFromVector RowMajor c v | otherwise = error $ "inconsistent list size = " @@ -205,22 +245,28 @@ fromLists = fromRows . map fromList -- | creates a 1-row matrix from a vector -asRow :: Element a => Vector a -> Matrix a +asRow :: Storable a => Vector a -> Matrix a asRow v = reshape (dim v) v -- | creates a 1-column matrix from a vector -asColumn :: Element a => Vector a -> Matrix a +asColumn :: Storable a => Vector a -> Matrix a asColumn v = reshape 1 v + {- | creates a Matrix of the specified size using the supplied function to - to map the row/column position to the value at that row/column position. + to map the row\/column position to the value at that row\/column position. @> buildMatrix 3 4 (\ (r,c) -> fromIntegral r * fromIntegral c) (3><4) [ 0.0, 0.0, 0.0, 0.0, 0.0 , 0.0, 1.0, 2.0, 3.0, 4.0 , 0.0, 2.0, 4.0, 6.0, 8.0]@ + +Hilbert matrix of order N: + +@hilb n = buildMatrix n n (\(i,j)->1/(fromIntegral i + fromIntegral j +1))@ + -} buildMatrix :: Element a => Int -> Int -> ((Int, Int) -> a) -> Matrix a buildMatrix rc cc f = @@ -229,135 +275,13 @@ ----------------------------------------------------- -fromArray2D :: (Element e) => Array (Int, Int) e -> Matrix e +fromArray2D :: (Storable e) => Array (Int, Int) e -> Matrix e fromArray2D m = (r> 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. - -@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 - -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" - --- | Tool to display matrices with latex syntax. -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). --} -fileDimensions :: FilePath -> IO (Int,Int) -fileDimensions fname = do - wcres <- readProcess "wc" ["-w",fname] "" - contents <- readFile fname - let tot = read . head . words $ wcres - c = length . head . dropWhile null . map words . lines $ contents - if tot > 0 - then return (tot `div` c, c) - else return (0,0) - --- | Loads a matrix from an ASCII file formatted as a 2D table. -loadMatrix :: FilePath -> IO (Matrix Double) -loadMatrix file = fromFile file =<< fileDimensions file - --- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance). -fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) -fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c) - - -- | rearranges the rows of a matrix according to the order given in a list of integers. extractRows :: Element t => [Int] -> Matrix t -> Matrix t extractRows l m = fromRows $ extract (toRows $ m) l @@ -424,47 +348,3 @@ cs = replicate qc c ++ if rc > 0 then [rc] else [] ------------------------------------------------------------------- - --- | conversion utilities -class (Element e) => Container c e where - toComplex :: RealFloat e => (c e, c e) -> c (Complex e) - fromComplex :: RealFloat e => c (Complex e) -> (c e, c e) - comp :: RealFloat e => c e -> c (Complex e) - conj :: RealFloat e => c (Complex e) -> c (Complex e) - real :: c Double -> c e - complex :: c e -> c (Complex Double) - -instance Container Vector Double where - toComplex = toComplexV - fromComplex = fromComplexV - comp v = toComplex (v,constant 0 (dim v)) - conj = conjV - real = id - complex = comp - -instance Container Vector (Complex Double) where - toComplex = undefined -- can't match - fromComplex = undefined - comp = undefined - conj = undefined - real = comp - complex = id - -instance Container Matrix Double where - toComplex = uncurry $ liftMatrix2 $ curry toComplex - fromComplex z = (reshape c r, reshape c i) - where (r,i) = fromComplex (flatten z) - c = cols z - comp = liftMatrix comp - conj = liftMatrix conj - real = id - complex = comp - -instance Container Matrix (Complex Double) where - toComplex = undefined - fromComplex = undefined - comp = undefined - conj = undefined - real = comp - complex = id - diff --git a/lib/Data/Packed/Random.hs b/lib/Data/Packed/Random.hs --- a/lib/Data/Packed/Random.hs +++ b/lib/Data/Packed/Random.hs @@ -20,11 +20,11 @@ ) where import Numeric.GSL.Vector -import Data.Packed.Matrix -import Data.Packed.Vector +import Data.Packed +import Numeric.ContainerBoot import Numeric.LinearAlgebra.Algorithms -import Numeric.LinearAlgebra.Interface + -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate -- Gaussian distribution. gaussianSample :: Int -- ^ seed @@ -34,9 +34,9 @@ -> Matrix Double -- ^ result gaussianSample seed n med cov = m where c = dim med - meds = constant 1 n `outer` med + meds = konst 1 n `outer` med rs = reshape c $ randomVector seed Gaussian (c * n) - m = rs <> cholSH cov + meds + m = rs `mXm` cholSH cov `add` meds -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate -- uniform distribution. @@ -50,8 +50,8 @@ cs = zipWith subtract as bs d = dim a dat = toRows $ reshape n $ randomVector seed Uniform (n*d) - am = constant 1 n `outer` a - m = fromColumns (zipWith scale cs dat) + am + am = konst 1 n `outer` a + m = fromColumns (zipWith scale cs dat) `add` am ------------ utilities ------------------------------- @@ -60,7 +60,7 @@ meanCov x = (med,cov) where r = rows x k = 1 / fromIntegral r - med = constant k r <> x - meds = constant 1 r `outer` med - xc = x - meds - cov = (trans xc <> xc) / fromIntegral (r-1) + med = konst k r `vXm` x + meds = konst 1 r `outer` med + xc = x `sub` meds + cov = flip scale (trans xc `mXm` xc) (recip (fromIntegral (r-1))) diff --git a/lib/Data/Packed/ST.hs b/lib/Data/Packed/ST.hs --- a/lib/Data/Packed/ST.hs +++ b/lib/Data/Packed/ST.hs @@ -90,11 +90,11 @@ writeVector = safeIndexV unsafeWriteVector {-# NOINLINE newUndefinedVector #-} -newUndefinedVector :: Element t => Int -> ST s (STVector s t) +newUndefinedVector :: Storable t => Int -> ST s (STVector s t) newUndefinedVector = unsafeIOToST . fmap STVector . createVector {-# INLINE newVector #-} -newVector :: Element t => t -> Int -> ST s (STVector s t) +newVector :: Storable t => t -> Int -> ST s (STVector s t) newVector x n = do v <- newUndefinedVector n let go (-1) = return v @@ -164,9 +164,9 @@ writeMatrix = safeIndexM unsafeWriteMatrix {-# NOINLINE newUndefinedMatrix #-} -newUndefinedMatrix :: Element t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t) +newUndefinedMatrix :: Storable t => MatrixOrder -> Int -> Int -> ST s (STMatrix s t) newUndefinedMatrix order r c = unsafeIOToST $ fmap STMatrix $ createMatrix order r c {-# NOINLINE newMatrix #-} -newMatrix :: Element t => t -> Int -> Int -> ST s (STMatrix s t) +newMatrix :: Storable t => t -> Int -> Int -> ST s (STMatrix s t) newMatrix v r c = unsafeThawMatrix $ reshape c $ runSTVector $ newVector v (r*c) diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs --- a/lib/Data/Packed/Vector.hs +++ b/lib/Data/Packed/Vector.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE FlexibleContexts #-} ----------------------------------------------------------------------------- -- | -- Module : Data.Packed.Vector @@ -10,6 +11,8 @@ -- -- 1D arrays suitable for numeric computations using external libraries. -- +-- This module provides basic functions for manipulation of structure. +-- ----------------------------------------------------------------------------- module Data.Packed.Vector ( @@ -17,53 +20,47 @@ fromList, (|>), toList, buildVector, dim, (@>), subVector, takesV, join, - constant, linspace, - vecdisp, - vectorMax, vectorMin, vectorMaxIndex, vectorMinIndex, - mapVector, zipVector, - fscanfVector, fprintfVector, freadVector, fwriteVector, - foldLoop, foldVector, foldVectorG + mapVector, zipVector, zipVectorWith, unzipVector, unzipVectorWith, + mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_, + foldLoop, foldVector, foldVectorG, foldVectorWithIndex ) where -import Data.Packed.Internal -import Numeric.GSL.Vector --- import Data.Packed.ST - -{- | Creates a real vector containing a range of values: - -@\> linspace 5 (-3,7) -5 |> [-3.0,-0.5,2.0,4.5,7.0]@ - -Logarithmic spacing can be defined as follows: +import Data.Packed.Internal.Vector +import Data.Binary +import Foreign.Storable +import Control.Monad(replicateM) -@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] - where scale = vectorMapValR Scale - add = vectorMapValR AddConstant - s = (b-a)/fromIntegral (n-1) +------------------------------------------------------------------- -vectorMax :: Vector Double -> Double -vectorMax = toScalarR Max +-- a 64K cache, with a Double taking 13 bytes in Bytestring, +-- implies a chunk size of 5041 +chunk :: Int +chunk = 5000 -vectorMin :: Vector Double -> Double -vectorMin = toScalarR Min +chunks :: Int -> [Int] +chunks d = let c = d `div` chunk + m = d `mod` chunk + in if m /= 0 then reverse (m:(replicate c chunk)) else (replicate c chunk) -vectorMaxIndex :: Vector Double -> Int -vectorMaxIndex = round . toScalarR MaxIdx +putVector v = do + let d = dim v + mapM_ (\i -> put $ v @> i) [0..(d-1)] -vectorMinIndex :: Vector Double -> Int -vectorMinIndex = round . toScalarR MinIdx +getVector d = do + xs <- replicateM d get + return $! fromList xs -{- | creates a vector with a given number of equal components: +instance (Binary a, Storable a) => Binary (Vector a) where + put v = do + let d = dim v + put d + mapM_ putVector $! takesV (chunks d) v + get = do + d <- get + vs <- mapM getVector $ chunks d + return $! join vs -@> constant 2 7 -7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ --} -constant :: Element a => a -> Int -> Vector a --- constant x n = runSTVector (newVector x n) -constant = constantD -- about 2x faster +------------------------------------------------------------------- {- | creates a Vector of the specified length using the supplied function to to map the index to the value at that index. @@ -72,22 +69,17 @@ 4 |> [0.0,1.0,2.0,3.0]@ -} -buildVector :: Element a => Int -> (Int -> a) -> Vector a +buildVector :: Storable a => Int -> (Int -> a) -> Vector a buildVector len f = fromList $ map f [0 .. (len - 1)] -{- | Show a vector using a function for showing matrices. +-- | zip for Vectors +zipVector :: (Storable a, Storable b, Storable (a,b)) => Vector a -> Vector b -> Vector (a,b) +zipVector = zipVectorWith (,) -@disp = putStr . vecdisp ('dispf' 2) +-- | unzip for Vectors +unzipVector :: (Storable a, Storable b, Storable (a,b)) => Vector (a,b) -> (Vector a,Vector b) +unzipVector = unzipVectorWith id -\> 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 +------------------------------------------------------------------- diff --git a/lib/Graphics/Plot.hs b/lib/Graphics/Plot.hs --- a/lib/Graphics/Plot.hs +++ b/lib/Graphics/Plot.hs @@ -3,15 +3,13 @@ -- Module : Graphics.Plot -- Copyright : (c) Alberto Ruiz 2005-8 -- License : GPL-style --- +-- -- Maintainer : Alberto Ruiz (aruiz at um dot es) -- Stability : provisional -- Portability : uses gnuplot and ImageMagick -- --- Very basic (and provisional) drawing tools using gnuplot and imageMagick. --- --- This module is deprecated. It will be replaced by improved drawing tools based --- on the Gnuplot package by Henning Thielemann. +-- This module is deprecated. It can be replaced by improved drawing tools +-- available in the plot\\plot-gtk packages by Vivian McPhail or Gnuplot by Henning Thielemann. ----------------------------------------------------------------------------- module Graphics.Plot( @@ -20,7 +18,7 @@ plot, parametricPlot, - splot, mesh, mesh', meshdom, + splot, mesh, meshdom, matrixToPGM, imshow, @@ -28,40 +26,18 @@ ) where -import Data.Packed -import Numeric.LinearAlgebra(outer) +import Numeric.Container import Data.List(intersperse) import System.Process (system) -size = dim - --- | Loads a real matrix from a formatted ASCII text file ---fromFile :: FilePath -> IO Matrix ---fromFile filename = readFile filename >>= return . readMatrix read - --- | Saves a real matrix to a formatted ascii text file -toFile' :: FilePath -> Matrix Double -> IO () -toFile' filename matrix = writeFile filename (unlines . map unwords. map (map show) . toLists $ matrix) - ------------------------------------------------------------------------- - - -- | From vectors x and y, it generates a pair of matrices to be used as x and y arguments for matrix functions. meshdom :: Vector Double -> Vector Double -> (Matrix Double , Matrix Double) -meshdom r1 r2 = (outer r1 (constant 1 (size r2)), outer (constant 1 (size r1)) r2) - -gnuplotX :: String -> IO () -gnuplotX command = do { _ <- system cmdstr; return()} where - cmdstr = "echo \""++command++"\" | gnuplot -persist" - -datafollows = "\\\"-\\\"" - -prep = (++"e\n\n") . unlines . map (unwords . (map show)) +meshdom r1 r2 = (outer r1 (constant 1 (dim r2)), outer (constant 1 (dim r1)) r2) {- | Draws a 3D surface representation of a real matrix. -> > mesh (hilb 20) +> > mesh $ build (10,10) (\\i j -> i + (j-5)^2) In certain versions you can interactively rotate the graphic using the mouse. @@ -71,15 +47,6 @@ command = "splot "++datafollows++" matrix with lines\n" dat = prep $ toLists $ m -mesh' :: Matrix Double -> IO () -mesh' m = do - writeFile "splot-gnu-command" "splot \"splot-tmp.txt\" matrix with lines; pause -1"; - toFile' "splot-tmp.txt" m - putStr "Press [Return] to close the graphic and continue... " - _ <- system "gnuplot -persist splot-gnu-command" - _ <- system "rm splot-tmp.txt splot-gnu-command" - return () - {- | Draws the surface represented by the function f in the desired ranges and number of points, internally using 'mesh'. > > let f x y = cos (x + y) @@ -87,11 +54,15 @@ -} splot :: (Matrix Double->Matrix Double->Matrix Double) -> (Double,Double) -> (Double,Double) -> Int -> IO () -splot f rx ry n = mesh' z where +splot f rx ry n = mesh z where (x,y) = meshdom (linspace n rx) (linspace n ry) z = f x y -{- | plots several vectors against the first one -} +{- | plots several vectors against the first one + +> > let t = linspace 100 (-3,3) in mplot [t, sin t, exp (-t^2)] + +-} mplot :: [Vector Double] -> IO () mplot m = gnuplotX (commands++dats) where commands = if length m == 1 then command1 else commandmore @@ -103,26 +74,6 @@ dats = concat (replicate (length m-1) dat) -{- -mplot' m = do - writeFile "plot-gnu-command" (commands++endcmd) - toFile "plot-tmp.txt" (fromColumns m) - putStr "Press [Return] to close the graphic and continue... " - system "gnuplot plot-gnu-command" - system "rm plot-tmp.txt plot-gnu-command" - return () - where - commands = if length m == 1 then command1 else commandmore - command1 = "plot \"plot-tmp.txt\" with lines\n" - commandmore = "plot " ++ plots ++ "\n" - plots = concat $ intersperse ", " (map cmd [2 .. length m]) - cmd k = "\"plot-tmp.txt\" using 1:"++show k++" with lines" - endcmd = "pause -1" --} - --- apply several functions to one object -mapf fs x = map ($ x) fs - {- | Draws a list of functions over a desired range and with a desired number of points > > plot [sin, cos, sin.(3*)] (0,2*pi) 1000 @@ -130,7 +81,8 @@ -} plot :: [Vector Double->Vector Double] -> (Double,Double) -> Int -> IO () plot fs rx n = mplot (x: mapf fs x) - where x = linspace n rx + where x = linspace n rx + mapf gs y = map ($ y) gs {- | Draws a parametric curve. For instance, to draw a spiral we can do something like: @@ -150,12 +102,12 @@ r = rows m header = "P2 "++show c++" "++show r++" "++show (round maxgray :: Int)++"\n" maxgray = 255.0 - maxval = vectorMax $ flatten $ m - minval = vectorMin $ flatten $ m - scale = if (maxval == minval) + maxval = maxElement m + minval = minElement m + scale' = if (maxval == minval) then 0.0 else maxgray / (maxval - minval) - f x = show ( round ( scale *(x - minval) ) :: Int ) + f x = show ( round ( scale' *(x - minval) ) :: Int ) ll = map (map f) (toLists m) -- | imshow shows a representation of a matrix as a gray level image using ImageMagick's display. @@ -165,6 +117,15 @@ return () ---------------------------------------------------- + +gnuplotX :: String -> IO () +gnuplotX command = do { _ <- system cmdstr; return()} where + cmdstr = "echo \""++command++"\" | gnuplot -persist" + +datafollows = "\\\"-\\\"" + +prep = (++"e\n\n") . unlines . map (unwords . (map show)) + gnuplotpdf :: String -> String -> [([[Double]], String)] -> IO () gnuplotpdf title command ds = gnuplot (prelude ++ command ++" "++ draw) >> postproc where diff --git a/lib/Numeric/Chain.hs b/lib/Numeric/Chain.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Chain.hs @@ -0,0 +1,140 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Chain +-- Copyright : (c) Vivian McPhail 2010 +-- License : GPL-style +-- +-- Maintainer : Vivian McPhail gmail.com> +-- Stability : provisional +-- Portability : portable +-- +-- optimisation of association order for chains of matrix multiplication +-- +----------------------------------------------------------------------------- + +module Numeric.Chain ( + optimiseMult, + ) where + +import Data.Maybe + +import Data.Packed.Matrix +import Numeric.ContainerBoot + +import qualified Data.Array.IArray as A + +----------------------------------------------------------------------------- +{- | + Provide optimal association order for a chain of matrix multiplications + and apply the multiplications. + + The algorithm is the well-known O(n\^3) dynamic programming algorithm + that builds a pyramid of optimal associations. + +> m1, m2, m3, m4 :: Matrix Double +> m1 = (10><15) [1..] +> m2 = (15><20) [1..] +> m3 = (20><5) [1..] +> m4 = (5><10) [1..] + +> >>> optimiseMult [m1,m2,m3,m4] + +will perform @((m1 `multiply` (m2 `multiply` m3)) `multiply` m4)@ + +The naive left-to-right multiplication would take @4500@ scalar multiplications +whereas the optimised version performs @2750@ scalar multiplications. The complexity +in this case is 32 (= 4^3/2) * (2 comparisons, 3 scalar multiplications, 3 scalar additions, +5 lookups, 2 updates) + a constant (= three table allocations) +-} +optimiseMult :: Product t => [Matrix t] -> Matrix t +optimiseMult = chain + +----------------------------------------------------------------------------- + +type Matrices a = A.Array Int (Matrix a) +type Sizes = A.Array Int (Int,Int) +type Cost = A.Array Int (A.Array Int (Maybe Int)) +type Indexes = A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int)))) + +update :: A.Array Int (A.Array Int a) -> (Int,Int) -> a -> A.Array Int (A.Array Int a) +update a (r,c) e = a A.// [(r,(a A.! r) A.// [(c,e)])] + +newWorkSpaceCost :: Int -> A.Array Int (A.Array Int (Maybe Int)) +newWorkSpaceCost n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n] + where subArray i = A.listArray (1,i) (repeat Nothing) + +newWorkSpaceIndexes :: Int -> A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int)))) +newWorkSpaceIndexes n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n] + where subArray i = A.listArray (1,i) (repeat Nothing) + +matricesToSizes :: [Matrix a] -> Sizes +matricesToSizes ms = A.listArray (1,length ms) $ map (\m -> (rows m,cols m)) ms + +chain :: Product a => [Matrix a] -> Matrix a +chain [] = error "chain: zero matrices to multiply" +chain [m] = m +chain [ml,mr] = ml `multiply` mr +chain ms = let ln = length ms + ma = A.listArray (1,ln) ms + mz = matricesToSizes ms + i = chain_cost mz + in chain_paren (ln,ln) i ma + +chain_cost :: Sizes -> Indexes +chain_cost mz = let (_,u) = A.bounds mz + cost = newWorkSpaceCost u + ixes = newWorkSpaceIndexes u + (_,_,i) = foldl chain_cost' (mz,cost,ixes) (order u) + in i + +chain_cost' :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes) +chain_cost' sci@(mz,cost,ixes) (r,c) + | c == 1 = let cost' = update cost (r,c) (Just 0) + ixes' = update ixes (r,c) (Just ((r,c),(r,c))) + in (mz,cost',ixes') + | otherwise = minimum_cost sci (r,c) + +minimum_cost :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes) +minimum_cost sci fu = foldl (smaller_cost fu) sci (fulcrum_order fu) + +smaller_cost :: (Int,Int) -> (Sizes,Cost,Indexes) -> ((Int,Int),(Int,Int)) -> (Sizes,Cost,Indexes) +smaller_cost (r,c) (mz,cost,ixes) ix@((lr,lc),(rr,rc)) = let op_cost = (fromJust ((cost A.! lr) A.! lc)) + + (fromJust ((cost A.! rr) A.! rc)) + + ((fst $ mz A.! (lr-lc+1)) + *(snd $ mz A.! lc) + *(snd $ mz A.! rr)) + cost' = (cost A.! r) A.! c + in case cost' of + Nothing -> let cost'' = update cost (r,c) (Just op_cost) + ixes'' = update ixes (r,c) (Just ix) + in (mz,cost'',ixes'') + Just ct -> if op_cost < ct then + let cost'' = update cost (r,c) (Just op_cost) + ixes'' = update ixes (r,c) (Just ix) + in (mz,cost'',ixes'') + else (mz,cost,ixes) + + +fulcrum_order (r,c) = let fs' = zip (repeat r) [1..(c-1)] + in map (partner (r,c)) fs' + +partner (r,c) (a,b) = (((r-b),(c-b)),(a,b)) + +order 0 = [] +order n = (order (n-1)) ++ (zip (repeat n) [1..n]) + +chain_paren :: Product a => (Int,Int) -> Indexes -> Matrices a -> Matrix a +chain_paren (r,c) ixes ma = let ((lr,lc),(rr,rc)) = fromJust $ (ixes A.! r) A.! c + in if lr == rr && lc == rc then (ma A.! lr) + else (chain_paren (lr,lc) ixes ma) `multiply` (chain_paren (rr,rc) ixes ma) + +-------------------------------------------------------------------------- + +{- TESTS -} + +-- optimal association is ((m1*(m2*m3))*m4) +m1, m2, m3, m4 :: Matrix Double +m1 = (10><15) [1..] +m2 = (15><20) [1..] +m3 = (20><5) [1..] +m4 = (5><10) [1..] diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Container.hs @@ -0,0 +1,132 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Container +-- Copyright : (c) Alberto Ruiz 2010 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines. +-- +-- The 'Container' class is used to define optimized generic functions which work +-- on 'Vector' and 'Matrix' with real or complex elements. +-- +-- Some of these functions are also available in the instances of the standard +-- numeric Haskell classes provided by "Numeric.LinearAlgebra". +-- +----------------------------------------------------------------------------- + +module Numeric.Container ( + -- * Basic functions + module Data.Packed, + constant, linspace, + diag, ident, + ctrans, + -- * Generic operations + Container(..), + -- * Matrix product + Product(..), + optimiseMult, + mXm,mXv,vXm,(<.>),(<>),(<\>), + outer, kronecker, + -- * Random numbers + RandDist(..), + randomVector, + gaussianSample, + uniformSample, + meanCov, + -- * Element conversion + Convert(..), + Complexable(), + RealElement(), + + RealOf, ComplexOf, SingleOf, DoubleOf, + + IndexOf, + module Data.Complex, + -- * Input / Output + dispf, disps, dispcf, vecdisp, latexFormat, format, + loadMatrix, saveMatrix, fromFile, fileDimensions, + readMatrix, + fscanfVector, fprintfVector, freadVector, fwriteVector, + -- * Experimental + build', konst', + -- * Deprecated + (.*),(*/),(<|>),(<->), + vectorMax,vectorMin, + vectorMaxIndex, vectorMinIndex +) where + +import Data.Packed +import Data.Packed.Internal(constantD) +import Numeric.ContainerBoot +import Numeric.Chain +import Numeric.IO +import Data.Complex +import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD) +import Data.Packed.Random + +------------------------------------------------------------------ + +{- | creates a vector with a given number of equal components: + +@> constant 2 7 +7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ +-} +constant :: Element a => a -> Int -> Vector a +-- constant x n = runSTVector (newVector x n) +constant = constantD-- about 2x faster + +{- | Creates a real vector containing a range of values: + +@\> 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 :: (Enum e, Container Vector e) => Int -> (e, e) -> Vector e +linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1] + where s = (b-a)/fromIntegral (n-1) + +-- | Dot product: @u \<.\> v = dot u v@ +(<.>) :: Product t => Vector t -> Vector t -> t +infixl 7 <.> +(<.>) = dot + + + +-------------------------------------------------------- + +class Mul a b c | a b -> c where + infixl 7 <> + -- | Matrix-matrix, matrix-vector, and vector-matrix products. + (<>) :: Product t => a t -> b t -> c t + +instance Mul Matrix Matrix Matrix where + (<>) = mXm + +instance Mul Matrix Vector Vector where + (<>) m v = flatten $ m <> (asColumn v) + +instance Mul Vector Matrix Vector where + (<>) v m = flatten $ (asRow v) <> m + +-------------------------------------------------------- + +-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD). +(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a +infixl 7 <\> +m <\> v = flatten (linearSolveSVD m (reshape 1 v)) + +-------------------------------------------------------- diff --git a/lib/Numeric/ContainerBoot.hs b/lib/Numeric/ContainerBoot.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/ContainerBoot.hs @@ -0,0 +1,583 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.ContainerBoot +-- Copyright : (c) Alberto Ruiz 2010 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Module to avoid cyclyc dependencies. +-- +----------------------------------------------------------------------------- + +module Numeric.ContainerBoot ( + -- * Basic functions + ident, diag, ctrans, + -- * Generic operations + Container(..), + -- * Matrix product and related functions + Product(..), + mXm,mXv,vXm, + outer, kronecker, + -- * Element conversion + Convert(..), + Complexable(), + RealElement(), + + RealOf, ComplexOf, SingleOf, DoubleOf, + + IndexOf, + module Data.Complex, + -- * Experimental + build', konst', + -- * Deprecated + (.*),(*/),(<|>),(<->), + vectorMax,vectorMin, + vectorMaxIndex, vectorMinIndex +) where + +import Data.Packed +import Numeric.Conversion +import Data.Packed.Internal +import Numeric.GSL.Vector + +import Data.Complex +import Control.Monad(ap) + +import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) + +import System.IO.Unsafe + +------------------------------------------------------------------- + +type family IndexOf c + +type instance IndexOf Vector = Int +type instance IndexOf Matrix = (Int,Int) + +type family ArgOf c a + +type instance ArgOf Vector a = a -> a +type instance ArgOf Matrix a = a -> a -> a + +------------------------------------------------------------------- + +-- | Basic element-by-element functions for numeric containers +class (Complexable c, Fractional e, Element e) => Container c e where + -- | create a structure with a single element + scalar :: e -> c e + -- | complex conjugate + conj :: c 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 + -- | element by element multiplication + mul :: c e -> c e -> c e + -- | element by element division + divide :: c e -> c e -> c e + equal :: c e -> c e -> Bool + -- + -- element by element inverse tangent + arctan2 :: c e -> c e -> c e + -- + -- | cannot implement instance Functor because of Element class constraint + cmap :: (Element a, Element b) => (a -> b) -> c a -> c b + -- | constant structure of given size + konst :: e -> IndexOf c -> c e + -- | create a structure using a function + -- + -- Hilbert matrix of order N: + -- + -- @hilb n = build (n,n) (\\i j -> 1/(i+j+1))@ + build :: IndexOf c -> (ArgOf c e) -> c e + --build :: BoundsOf f -> f -> (ContainerOf f) e + -- + -- | indexing function + atIndex :: c e -> IndexOf c -> e + -- | index of min element + minIndex :: c e -> IndexOf c + -- | index of max element + maxIndex :: c e -> IndexOf c + -- | value of min element + minElement :: c e -> e + -- | value of max element + maxElement :: c e -> e + -- the C functions sumX/prodX are twice as fast as using foldVector + -- | the sum of elements (faster than using @fold@) + sumElements :: c e -> e + -- | the product of elements (faster than using @fold@) + prodElements :: c e -> e + +-------------------------------------------------------------------------- + +instance Container Vector Float where + scale = vectorMapValF Scale + scaleRecip = vectorMapValF Recip + addConstant = vectorMapValF AddConstant + add = vectorZipF Add + sub = vectorZipF Sub + mul = vectorZipF Mul + divide = vectorZipF Div + equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 + arctan2 = vectorZipF ATan2 + scalar x = fromList [x] + konst = constantD + build = buildV + conj = id + cmap = mapVector + atIndex = (@>) + minIndex = round . toScalarF MinIdx + maxIndex = round . toScalarF MaxIdx + minElement = toScalarF Min + maxElement = toScalarF Max + sumElements = sumF + prodElements = prodF + +instance Container Vector Double where + scale = vectorMapValR Scale + scaleRecip = vectorMapValR Recip + addConstant = vectorMapValR AddConstant + add = vectorZipR Add + sub = vectorZipR Sub + mul = vectorZipR Mul + divide = vectorZipR Div + equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 + arctan2 = vectorZipR ATan2 + scalar x = fromList [x] + konst = constantD + build = buildV + conj = id + cmap = mapVector + atIndex = (@>) + minIndex = round . toScalarR MinIdx + maxIndex = round . toScalarR MaxIdx + minElement = toScalarR Min + maxElement = toScalarR Max + sumElements = sumR + prodElements = prodR + +instance Container Vector (Complex Double) where + scale = vectorMapValC Scale + scaleRecip = vectorMapValC Recip + addConstant = vectorMapValC AddConstant + add = vectorZipC Add + sub = vectorZipC Sub + mul = vectorZipC Mul + divide = vectorZipC Div + equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 + arctan2 = vectorZipC ATan2 + scalar x = fromList [x] + konst = constantD + build = buildV + conj = conjugateC + cmap = mapVector + atIndex = (@>) + minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) + maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) + minElement = ap (@>) minIndex + maxElement = ap (@>) maxIndex + sumElements = sumC + prodElements = prodC + +instance Container Vector (Complex Float) where + scale = vectorMapValQ Scale + scaleRecip = vectorMapValQ Recip + addConstant = vectorMapValQ AddConstant + add = vectorZipQ Add + sub = vectorZipQ Sub + mul = vectorZipQ Mul + divide = vectorZipQ Div + equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 + arctan2 = vectorZipQ ATan2 + scalar x = fromList [x] + konst = constantD + build = buildV + conj = conjugateQ + cmap = mapVector + atIndex = (@>) + minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) + maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) + minElement = ap (@>) minIndex + maxElement = ap (@>) maxIndex + sumElements = sumQ + prodElements = prodQ + +--------------------------------------------------------------- + +instance (Container Vector a) => Container Matrix a where + scale x = liftMatrix (scale x) + scaleRecip x = liftMatrix (scaleRecip x) + addConstant x = liftMatrix (addConstant x) + add = liftMatrix2 add + sub = liftMatrix2 sub + mul = liftMatrix2 mul + divide = liftMatrix2 divide + equal a b = cols a == cols b && flatten a `equal` flatten b + arctan2 = liftMatrix2 arctan2 + scalar x = (1><1) [x] + konst v (r,c) = reshape c (konst v (r*c)) + build = buildM + conj = liftMatrix conj + cmap f = liftMatrix (mapVector f) + atIndex = (@@>) + minIndex m = let (r,c) = (rows m,cols m) + i = (minIndex $ flatten m) + in (i `div` c,i `mod` c) + maxIndex m = let (r,c) = (rows m,cols m) + i = (maxIndex $ flatten m) + in (i `div` c,i `mod` c) + minElement = ap (@@>) minIndex + maxElement = ap (@@>) maxIndex + sumElements = sumElements . flatten + prodElements = prodElements . flatten + +---------------------------------------------------- + +-- | Matrix product and related functions +class Element e => Product e where + -- | matrix product + multiply :: Matrix e -> Matrix e -> Matrix e + -- | dot (inner) product + dot :: Vector e -> Vector e -> e + -- | sum of absolute value of elements (differs in complex case from @norm1@) + absSum :: Vector e -> RealOf e + -- | sum of absolute value of elements + norm1 :: Vector e -> RealOf e + -- | euclidean norm + norm2 :: Vector e -> RealOf e + -- | element of maximum magnitude + normInf :: Vector e -> RealOf e + +instance Product Float where + norm2 = toScalarF Norm2 + absSum = toScalarF AbsSum + dot = dotF + norm1 = toScalarF AbsSum + normInf = maxElement . vectorMapF Abs + multiply = multiplyF + +instance Product Double where + norm2 = toScalarR Norm2 + absSum = toScalarR AbsSum + dot = dotR + norm1 = toScalarR AbsSum + normInf = maxElement . vectorMapR Abs + multiply = multiplyR + +instance Product (Complex Float) where + norm2 = toScalarQ Norm2 + absSum = toScalarQ AbsSum + dot = dotQ + norm1 = sumElements . fst . fromComplex . vectorMapQ Abs + normInf = maxElement . fst . fromComplex . vectorMapQ Abs + multiply = multiplyQ + +instance Product (Complex Double) where + norm2 = toScalarC Norm2 + absSum = toScalarC AbsSum + dot = dotC + norm1 = sumElements . fst . fromComplex . vectorMapC Abs + normInf = maxElement . fst . fromComplex . vectorMapC Abs + multiply = multiplyC + +---------------------------------------------------------- + +-- synonym for matrix product +mXm :: Product t => Matrix t -> Matrix t -> Matrix t +mXm = multiply + +-- matrix - vector product +mXv :: Product t => Matrix t -> Vector t -> Vector t +mXv m v = flatten $ m `mXm` (asColumn v) + +-- vector - matrix product +vXm :: Product t => Vector t -> Matrix t -> Vector t +vXm v m = flatten $ (asRow v) `mXm` m + +{- | Outer product of two vectors. + +@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3] +(3><3) + [ 5.0, 2.0, 3.0 + , 10.0, 4.0, 6.0 + , 15.0, 6.0, 9.0 ]@ +-} +outer :: (Product t) => Vector t -> Vector t -> Matrix t +outer u v = asColumn u `multiply` asRow v + +{- | Kronecker product of two matrices. + +@m1=(2><3) + [ 1.0, 2.0, 0.0 + , 0.0, -1.0, 3.0 ] +m2=(4><3) + [ 1.0, 2.0, 3.0 + , 4.0, 5.0, 6.0 + , 7.0, 8.0, 9.0 + , 10.0, 11.0, 12.0 ]@ + +@\> kronecker m1 m2 +(8><9) + [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0 + , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0 + , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0 + , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0 + , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0 + , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0 + , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 + , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@ +-} +kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t +kronecker a b = fromBlocks + . splitEvery (cols a) + . map (reshape (cols b)) + . toRows + $ flatten a `outer` flatten b + +------------------------------------------------------------------- + + +class Convert t where + real :: Container c t => c (RealOf t) -> c t + complex :: Container c t => c t -> c (ComplexOf t) + single :: Container c t => c t -> c (SingleOf t) + double :: Container c t => c t -> c (DoubleOf t) + toComplex :: (Container c t, RealElement t) => (c t, c t) -> c (Complex t) + fromComplex :: (Container c t, RealElement t) => c (Complex t) -> (c t, c t) + + +instance Convert Double where + real = id + complex = comp' + single = single' + double = id + toComplex = toComplex' + fromComplex = fromComplex' + +instance Convert Float where + real = id + complex = comp' + single = id + double = double' + toComplex = toComplex' + fromComplex = fromComplex' + +instance Convert (Complex Double) where + real = comp' + complex = id + single = single' + double = id + toComplex = toComplex' + fromComplex = fromComplex' + +instance Convert (Complex Float) where + real = comp' + complex = id + single = id + double = double' + toComplex = toComplex' + fromComplex = fromComplex' + +------------------------------------------------------------------- + +type family RealOf x + +type instance RealOf Double = Double +type instance RealOf (Complex Double) = Double + +type instance RealOf Float = Float +type instance RealOf (Complex Float) = Float + +type family ComplexOf x + +type instance ComplexOf Double = Complex Double +type instance ComplexOf (Complex Double) = Complex Double + +type instance ComplexOf Float = Complex Float +type instance ComplexOf (Complex Float) = Complex Float + +type family SingleOf x + +type instance SingleOf Double = Float +type instance SingleOf Float = Float + +type instance SingleOf (Complex a) = Complex (SingleOf a) + +type family DoubleOf x + +type instance DoubleOf Double = Double +type instance DoubleOf Float = Double + +type instance DoubleOf (Complex a) = Complex (DoubleOf a) + +type family ElementOf c + +type instance ElementOf (Vector a) = a +type instance ElementOf (Matrix a) = a + +------------------------------------------------------------ + +conjugateAux fun x = unsafePerformIO $ do + v <- createVector (dim x) + app2 fun vec x vec v "conjugateAux" + return v + +conjugateQ :: Vector (Complex Float) -> Vector (Complex Float) +conjugateQ = conjugateAux c_conjugateQ +foreign import ccall "conjugateQ" c_conjugateQ :: TQVQV + +conjugateC :: Vector (Complex Double) -> Vector (Complex Double) +conjugateC = conjugateAux c_conjugateC +foreign import ccall "conjugateC" c_conjugateC :: TCVCV + +---------------------------------------------------- + +{-# 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 +infixl 7 */ +v */ x = scale (recip x) v + + +------------------------------------------------ + +{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-} +{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-} + +class Joinable a b where + joinH :: Element t => a t -> b t -> Matrix t + joinV :: Element t => a t -> b t -> Matrix t + +instance Joinable Matrix Matrix where + joinH m1 m2 = fromBlocks [[m1,m2]] + joinV m1 m2 = fromBlocks [[m1],[m2]] + +instance Joinable Matrix Vector where + joinH m v = joinH m (asColumn v) + joinV m v = joinV m (asRow v) + +instance Joinable Vector Matrix where + joinH v m = joinH (asColumn v) m + joinV v m = joinV (asRow v) m + +infixl 4 <|> +infixl 3 <-> + +{-- - | Horizontal concatenation of matrices and vectors: + +@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0] +(6><4) + [ 1.0, 0.0, 0.0, 1.0 + , 0.0, 1.0, 0.0, 2.0 + , 0.0, 0.0, 1.0, 3.0 + , 3.0, 0.0, 0.0, 4.0 + , 0.0, 3.0, 0.0, 5.0 + , 0.0, 0.0, 3.0, 6.0 ]@ +-} +-- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t +a <|> b = joinH a b + +-- -- | Vertical concatenation of matrices and vectors. +-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t +a <-> b = joinV a b + +------------------------------------------------------------------- + +{-# DEPRECATED vectorMin "use minElement" #-} +vectorMin :: (Container Vector t, Element t) => Vector t -> t +vectorMin = minElement + +{-# DEPRECATED vectorMax "use maxElement" #-} +vectorMax :: (Container Vector t, Element t) => Vector t -> t +vectorMax = maxElement + + +{-# DEPRECATED vectorMaxIndex "use minIndex" #-} +vectorMaxIndex :: Vector Double -> Int +vectorMaxIndex = round . toScalarR MaxIdx + +{-# DEPRECATED vectorMinIndex "use maxIndex" #-} +vectorMinIndex :: Vector Double -> Int +vectorMinIndex = round . toScalarR MinIdx + +----------------------------------------------------- + +class Build f where + build' :: BoundsOf f -> f -> ContainerOf f + +type family BoundsOf x + +type instance BoundsOf (a->a) = Int +type instance BoundsOf (a->a->a) = (Int,Int) + +type family ContainerOf x + +type instance ContainerOf (a->a) = Vector a +type instance ContainerOf (a->a->a) = Matrix a + +instance (Element a, Num a) => Build (a->a) where + build' = buildV + +instance (Element a, Num a) => Build (a->a->a) where + build' = buildM + +buildM (rc,cc) f = fromLists [ [f r c | c <- cs] | r <- rs ] + where rs = map fromIntegral [0 .. (rc-1)] + cs = map fromIntegral [0 .. (cc-1)] + +buildV n f = fromList [f k | k <- ks] + where ks = map fromIntegral [0 .. (n-1)] + +---------------------------------------------------- +-- experimental + +class Konst s where + konst' :: Element e => e -> s -> ContainerOf' s e + +type family ContainerOf' x y + +type instance ContainerOf' Int a = Vector a +type instance ContainerOf' (Int,Int) a = Matrix a + +instance Konst Int where + konst' = constantD + +instance Konst (Int,Int) where + konst' k (r,c) = reshape c $ konst' k (r*c) + +-------------------------------------------------------- +-- | conjugate transpose +ctrans :: (Container Vector e, Element e) => Matrix e -> Matrix e +ctrans = liftMatrix conj . trans + +-- | Creates a square matrix with a given diagonal. +diag :: (Num a, Element a) => Vector a -> Matrix a +diag v = diagRect 0 v n n where n = dim v + +-- | creates the identity matrix of given dimension +ident :: (Num a, Element a) => Int -> Matrix a +ident n = diag (constantD 1 n) diff --git a/lib/Numeric/Conversion.hs b/lib/Numeric/Conversion.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Conversion.hs @@ -0,0 +1,91 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Conversion +-- Copyright : (c) Alberto Ruiz 2010 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Conversion routines +-- +----------------------------------------------------------------------------- + +module Numeric.Conversion ( + Complexable(..), RealElement, + module Data.Complex +) where + +import Data.Packed.Internal.Vector +import Data.Packed.Internal.Matrix +import Data.Complex +import Control.Arrow((***)) + +------------------------------------------------------------------- + +-- | Supported single-double precision type pairs +class (Element s, Element d) => Precision s d | s -> d, d -> s where + double2FloatG :: Vector d -> Vector s + float2DoubleG :: Vector s -> Vector d + +instance Precision Float Double where + double2FloatG = double2FloatV + float2DoubleG = float2DoubleV + +instance Precision (Complex Float) (Complex Double) where + double2FloatG = asComplex . double2FloatV . asReal + float2DoubleG = asComplex . float2DoubleV . asReal + +-- | Supported real types +class (Element t, Element (Complex t), RealFloat t +-- , RealOf t ~ t, RealOf (Complex t) ~ t + ) + => RealElement t + +instance RealElement Double +instance RealElement Float + + +-- | Structures that may contain complex numbers +class Complexable c where + toComplex' :: (RealElement e) => (c e, c e) -> c (Complex e) + fromComplex' :: (RealElement e) => c (Complex e) -> (c e, c e) + comp' :: (RealElement e) => c e -> c (Complex e) + single' :: Precision a b => c b -> c a + double' :: Precision a b => c a -> c b + + +instance Complexable Vector where + toComplex' = toComplexV + fromComplex' = fromComplexV + comp' v = toComplex' (v,constantD 0 (dim v)) + single' = double2FloatG + double' = float2DoubleG + + +-- | creates a complex vector from vectors with real and imaginary parts +toComplexV :: (RealElement a) => (Vector a, Vector a) -> Vector (Complex a) +toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i] + +-- | the inverse of 'toComplex' +fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a) +fromComplexV z = (r,i) where + [r,i] = toColumns $ reshape 2 $ asReal z + + +instance Complexable Matrix where + toComplex' = uncurry $ liftMatrix2 $ curry toComplex' + fromComplex' z = (reshape c *** reshape c) . fromComplex' . flatten $ z + where c = cols z + comp' = liftMatrix comp' + single' = liftMatrix single' + double' = liftMatrix double' + diff --git a/lib/Numeric/GSL/Fitting.hs b/lib/Numeric/GSL/Fitting.hs --- a/lib/Numeric/GSL/Fitting.hs +++ b/lib/Numeric/GSL/Fitting.hs @@ -109,7 +109,7 @@ sol = toList vsol c = max 1 (chi/sqrt (fromIntegral dof)) dof = length dat - (rows cov) - chi = pnorm PNorm2 (fromList $ cost (resMs model) dat sol) + chi = norm2 (fromList $ cost (resMs model) dat sol) js = fromLists $ jacobian (resDs deriv) dat sol cov = inv $ trans js <> js errs = toList $ scalar c * sqrt (takeDiag cov) diff --git a/lib/Numeric/GSL/Vector.hs b/lib/Numeric/GSL/Vector.hs --- a/lib/Numeric/GSL/Vector.hs +++ b/lib/Numeric/GSL/Vector.hs @@ -14,10 +14,13 @@ ----------------------------------------------------------------------------- module Numeric.GSL.Vector ( - FunCodeS(..), toScalarR, - FunCodeV(..), vectorMapR, vectorMapC, - FunCodeSV(..), vectorMapValR, vectorMapValC, - FunCodeVV(..), vectorZipR, vectorZipC, + sumF, sumR, sumQ, sumC, + prodF, prodR, prodQ, prodC, + dotF, dotR, dotQ, dotC, + FunCodeS(..), toScalarR, toScalarF, toScalarC, toScalarQ, + FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ, + FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ, + FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ, RandDist(..), randomVector ) where @@ -76,6 +79,107 @@ ------------------------------------------------------------------ +-- | sum of elements +sumF :: Vector Float -> Float +sumF x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumF vec x vec r "sumF" + return $ r @> 0 + +-- | sum of elements +sumR :: Vector Double -> Double +sumR x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumR vec x vec r "sumR" + return $ r @> 0 + +-- | sum of elements +sumQ :: Vector (Complex Float) -> Complex Float +sumQ x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumQ vec x vec r "sumQ" + return $ r @> 0 + +-- | sum of elements +sumC :: Vector (Complex Double) -> Complex Double +sumC x = unsafePerformIO $ do + r <- createVector 1 + app2 c_sumC vec x vec r "sumC" + return $ r @> 0 + +foreign import ccall safe "gsl-aux.h sumF" c_sumF :: TFF +foreign import ccall safe "gsl-aux.h sumR" c_sumR :: TVV +foreign import ccall safe "gsl-aux.h sumQ" c_sumQ :: TQVQV +foreign import ccall safe "gsl-aux.h sumC" c_sumC :: TCVCV + +-- | product of elements +prodF :: Vector Float -> Float +prodF x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodF vec x vec r "prodF" + return $ r @> 0 + +-- | product of elements +prodR :: Vector Double -> Double +prodR x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodR vec x vec r "prodR" + return $ r @> 0 + +-- | product of elements +prodQ :: Vector (Complex Float) -> Complex Float +prodQ x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodQ vec x vec r "prodQ" + return $ r @> 0 + +-- | product of elements +prodC :: Vector (Complex Double) -> Complex Double +prodC x = unsafePerformIO $ do + r <- createVector 1 + app2 c_prodC vec x vec r "prodC" + return $ r @> 0 + +foreign import ccall safe "gsl-aux.h prodF" c_prodF :: TFF +foreign import ccall safe "gsl-aux.h prodR" c_prodR :: TVV +foreign import ccall safe "gsl-aux.h prodQ" c_prodQ :: TQVQV +foreign import ccall safe "gsl-aux.h prodC" c_prodC :: TCVCV + +-- | dot product +dotF :: Vector Float -> Vector Float -> Float +dotF x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotF vec x vec y vec r "dotF" + return $ r @> 0 + +-- | dot product +dotR :: Vector Double -> Vector Double -> Double +dotR x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotR vec x vec y vec r "dotR" + return $ r @> 0 + +-- | dot product +dotQ :: Vector (Complex Float) -> Vector (Complex Float) -> Complex Float +dotQ x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotQ vec x vec y vec r "dotQ" + return $ r @> 0 + +-- | dot product +dotC :: Vector (Complex Double) -> Vector (Complex Double) -> Complex Double +dotC x y = unsafePerformIO $ do + r <- createVector 1 + app3 c_dotC vec x vec y vec r "dotC" + return $ r @> 0 + +foreign import ccall safe "gsl-aux.h dotF" c_dotF :: TFFF +foreign import ccall safe "gsl-aux.h dotR" c_dotR :: TVVV +foreign import ccall safe "gsl-aux.h dotQ" c_dotQ :: TQVQVQV +foreign import ccall safe "gsl-aux.h dotC" c_dotC :: TCVCVCV + +------------------------------------------------------------------ + toScalarAux fun code v = unsafePerformIO $ do r <- createVector 1 app2 (fun (fromei code)) vec v vec r "toScalarAux" @@ -106,6 +210,24 @@ foreign import ccall safe "gsl-aux.h toScalarR" c_toScalarR :: CInt -> TVV +-- | obtains different functions of a vector: norm1, norm2, max, min, posmax, posmin, etc. +toScalarF :: FunCodeS -> Vector Float -> Float +toScalarF oper = toScalarAux c_toScalarF (fromei oper) + +foreign import ccall safe "gsl-aux.h toScalarF" c_toScalarF :: CInt -> TFF + +-- | obtains different functions of a vector: only norm1, norm2 +toScalarC :: FunCodeS -> Vector (Complex Double) -> Double +toScalarC oper = toScalarAux c_toScalarC (fromei oper) + +foreign import ccall safe "gsl-aux.h toScalarC" c_toScalarC :: CInt -> TCVV + +-- | obtains different functions of a vector: only norm1, norm2 +toScalarQ :: FunCodeS -> Vector (Complex Float) -> Float +toScalarQ oper = toScalarAux c_toScalarQ (fromei oper) + +foreign import ccall safe "gsl-aux.h toScalarQ" c_toScalarQ :: CInt -> TQVF + ------------------------------------------------------------------ -- | map of real vectors with given function @@ -120,6 +242,18 @@ foreign import ccall safe "gsl-aux.h mapC" c_vectorMapC :: CInt -> TCVCV +-- | map of real vectors with given function +vectorMapF :: FunCodeV -> Vector Float -> Vector Float +vectorMapF = vectorMapAux c_vectorMapF + +foreign import ccall safe "gsl-aux.h mapF" c_vectorMapF :: CInt -> TFF + +-- | map of real vectors with given function +vectorMapQ :: FunCodeV -> Vector (Complex Float) -> Vector (Complex Float) +vectorMapQ = vectorMapAux c_vectorMapQ + +foreign import ccall safe "gsl-aux.h mapQ" c_vectorMapQ :: CInt -> TQVQV + ------------------------------------------------------------------- -- | map of real vectors with given function @@ -134,6 +268,18 @@ foreign import ccall safe "gsl-aux.h mapValC" c_vectorMapValC :: CInt -> Ptr (Complex Double) -> TCVCV +-- | map of real vectors with given function +vectorMapValF :: FunCodeSV -> Float -> Vector Float -> Vector Float +vectorMapValF oper = vectorMapValAux c_vectorMapValF (fromei oper) + +foreign import ccall safe "gsl-aux.h mapValF" c_vectorMapValF :: CInt -> Ptr Float -> TFF + +-- | map of complex vectors with given function +vectorMapValQ :: FunCodeSV -> Complex Float -> Vector (Complex Float) -> Vector (Complex Float) +vectorMapValQ oper = vectorMapValAux c_vectorMapValQ (fromei oper) + +foreign import ccall safe "gsl-aux.h mapValQ" c_vectorMapValQ :: CInt -> Ptr (Complex Float) -> TQVQV + ------------------------------------------------------------------- -- | elementwise operation on real vectors @@ -147,6 +293,18 @@ vectorZipC = vectorZipAux c_vectorZipC foreign import ccall safe "gsl-aux.h zipC" c_vectorZipC :: CInt -> TCVCVCV + +-- | elementwise operation on real vectors +vectorZipF :: FunCodeVV -> Vector Float -> Vector Float -> Vector Float +vectorZipF = vectorZipAux c_vectorZipF + +foreign import ccall safe "gsl-aux.h zipF" c_vectorZipF :: CInt -> TFFF + +-- | elementwise operation on complex vectors +vectorZipQ :: FunCodeVV -> Vector (Complex Float) -> Vector (Complex Float) -> Vector (Complex Float) +vectorZipQ = vectorZipAux c_vectorZipQ + +foreign import ccall safe "gsl-aux.h zipQ" c_vectorZipQ :: CInt -> TQVQVQV ----------------------------------------------------------------------- diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c @@ -10,6 +10,16 @@ #define KCVEC(A) int A##n, const gsl_complex*A##p #define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p +#define FVEC(A) int A##n, float*A##p +#define FMAT(A) int A##r, int A##c, float* A##p +#define KFVEC(A) int A##n, const float*A##p +#define KFMAT(A) int A##r, int A##c, const float* A##p + +#define QVEC(A) int A##n, gsl_complex_float*A##p +#define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p +#define KQVEC(A) int A##n, const gsl_complex_float*A##p +#define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p + #include #include #include @@ -64,12 +74,24 @@ #define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n) #define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c) +#define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n) +#define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c) +#define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n) +#define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c) +#define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n) +#define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c) +#define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n) +#define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c) + #define V(a) (&a.vector) #define M(a) (&a.matrix) #define GCVEC(A) int A##n, gsl_complex*A##p #define KGCVEC(A) int A##n, const gsl_complex*A##p +#define GQVEC(A) int A##n, gsl_complex_float*A##p +#define KGQVEC(A) int A##n, const gsl_complex_float*A##p + #define BAD_SIZE 2000 #define BAD_CODE 2001 #define MEM 2002 @@ -81,6 +103,154 @@ } +int sumF(KFVEC(x),FVEC(r)) { + DEBUGMSG("sumF"); + REQUIRES(rn==1,BAD_SIZE); + int i; + float res = 0; + for (i = 0; i < xn; i++) res += xp[i]; + rp[0] = res; + OK +} + +int sumR(KRVEC(x),RVEC(r)) { + DEBUGMSG("sumR"); + REQUIRES(rn==1,BAD_SIZE); + int i; + double res = 0; + for (i = 0; i < xn; i++) res += xp[i]; + rp[0] = res; + OK +} + +int sumQ(KQVEC(x),QVEC(r)) { + DEBUGMSG("sumQ"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex_float res; + res.dat[0] = 0; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + res.dat[0] += xp[i].dat[0]; + res.dat[1] += xp[i].dat[1]; + } + rp[0] = res; + OK +} + +int sumC(KCVEC(x),CVEC(r)) { + DEBUGMSG("sumC"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex res; + res.dat[0] = 0; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + res.dat[0] += xp[i].dat[0]; + res.dat[1] += xp[i].dat[1]; + } + rp[0] = res; + OK +} + +int prodF(KFVEC(x),FVEC(r)) { + DEBUGMSG("prodF"); + REQUIRES(rn==1,BAD_SIZE); + int i; + float res = 1; + for (i = 0; i < xn; i++) res *= xp[i]; + rp[0] = res; + OK +} + +int prodR(KRVEC(x),RVEC(r)) { + DEBUGMSG("prodR"); + REQUIRES(rn==1,BAD_SIZE); + int i; + double res = 1; + for (i = 0; i < xn; i++) res *= xp[i]; + rp[0] = res; + OK +} + +int prodQ(KQVEC(x),QVEC(r)) { + DEBUGMSG("prodQ"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex_float res; + float temp; + res.dat[0] = 1; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; + res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; + res.dat[0] = temp; + } + rp[0] = res; + OK +} + +int prodC(KCVEC(x),CVEC(r)) { + DEBUGMSG("prodC"); + REQUIRES(rn==1,BAD_SIZE); + int i; + gsl_complex res; + double temp; + res.dat[0] = 1; + res.dat[1] = 0; + for (i = 0; i < xn; i++) { + temp = res.dat[0] * xp[i].dat[0] - res.dat[1] * xp[i].dat[1]; + res.dat[1] = res.dat[0] * xp[i].dat[1] + res.dat[1] * xp[i].dat[0]; + res.dat[0] = temp; + } + rp[0] = res; + OK +} + +int dotF(KFVEC(x), KFVEC(y), FVEC(r)) { + DEBUGMSG("dotF"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotF"); + KFVVIEW(x); + KFVVIEW(y); + gsl_blas_sdot(V(x),V(y),rp); + OK +} + +int dotR(KRVEC(x), KRVEC(y), RVEC(r)) { + DEBUGMSG("dotR"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotR"); + KDVVIEW(x); + KDVVIEW(y); + gsl_blas_ddot(V(x),V(y),rp); + OK +} + +int dotQ(KQVEC(x), KQVEC(y), QVEC(r)) { + DEBUGMSG("dotQ"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotQ"); + KQVVIEW(x); + KQVVIEW(y); + gsl_blas_cdotu(V(x),V(y),rp); + OK +} + +int dotC(KCVEC(x), KCVEC(y), CVEC(r)) { + DEBUGMSG("dotC"); + REQUIRES(xn==yn,BAD_SIZE); + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("dotC"); + KCVVIEW(x); + KCVVIEW(y); + gsl_blas_zdotu(V(x),V(y),rp); + OK +} + int toScalarR(int code, KRVEC(x), RVEC(r)) { REQUIRES(rn==1,BAD_SIZE); DEBUGMSG("toScalarR"); @@ -99,7 +269,54 @@ OK } +int toScalarF(int code, KFVEC(x), FVEC(r)) { + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("toScalarF"); + KFVVIEW(x); + float res; + switch(code) { + case 0: { res = gsl_blas_snrm2(V(x)); break; } + case 1: { res = gsl_blas_sasum(V(x)); break; } + case 2: { res = gsl_vector_float_max_index(V(x)); break; } + case 3: { res = gsl_vector_float_max(V(x)); break; } + case 4: { res = gsl_vector_float_min_index(V(x)); break; } + case 5: { res = gsl_vector_float_min(V(x)); break; } + default: ERROR(BAD_CODE); + } + rp[0] = res; + OK +} + +int toScalarC(int code, KCVEC(x), RVEC(r)) { + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("toScalarC"); + KCVVIEW(x); + double res; + switch(code) { + case 0: { res = gsl_blas_dznrm2(V(x)); break; } + case 1: { res = gsl_blas_dzasum(V(x)); break; } + default: ERROR(BAD_CODE); + } + rp[0] = res; + OK +} + +int toScalarQ(int code, KQVEC(x), FVEC(r)) { + REQUIRES(rn==1,BAD_SIZE); + DEBUGMSG("toScalarQ"); + KQVVIEW(x); + float res; + switch(code) { + case 0: { res = gsl_blas_scnrm2(V(x)); break; } + case 1: { res = gsl_blas_scasum(V(x)); break; } + default: ERROR(BAD_CODE); + } + rp[0] = res; + OK +} + + inline double sign(double x) { if(x>0) { return +1.0; @@ -110,6 +327,16 @@ } } +inline float float_sign(float x) { + if(x>0) { + return +1.0; + } else if (x<0) { + return -1.0; + } else { + return 0.0; + } +} + inline gsl_complex complex_abs(gsl_complex z) { gsl_complex r; r.dat[0] = gsl_complex_abs(z); @@ -159,7 +386,33 @@ } } +int mapF(int code, KFVEC(x), FVEC(r)) { + int k; + REQUIRES(xn == rn,BAD_SIZE); + DEBUGMSG("mapF"); + switch (code) { + OP(0,sin) + OP(1,cos) + OP(2,tan) + OP(3,fabs) + OP(4,asin) + OP(5,acos) + OP(6,atan) /* atan2 mediante vectorZip */ + OP(7,sinh) + OP(8,cosh) + OP(9,tanh) + OP(10,gsl_asinh) + OP(11,gsl_acosh) + OP(12,gsl_atanh) + OP(13,exp) + OP(14,log) + OP(15,sign) + OP(16,sqrt) + default: ERROR(BAD_CODE); + } +} + int mapCAux(int code, KGCVEC(x), GCVEC(r)) { int k; REQUIRES(xn == rn,BAD_SIZE); @@ -194,6 +447,83 @@ } +gsl_complex_float complex_float_math_fun(gsl_complex (*cf)(gsl_complex), gsl_complex_float a) +{ + gsl_complex c; + gsl_complex r; + + gsl_complex_float float_r; + + c.dat[0] = a.dat[0]; + c.dat[1] = a.dat[1]; + + r = (*cf)(c); + + float_r.dat[0] = r.dat[0]; + float_r.dat[1] = r.dat[1]; + + return float_r; +} + +gsl_complex_float complex_float_math_op(gsl_complex (*cf)(gsl_complex,gsl_complex), + gsl_complex_float a,gsl_complex_float b) +{ + gsl_complex c1; + gsl_complex c2; + gsl_complex r; + + gsl_complex_float float_r; + + c1.dat[0] = a.dat[0]; + c1.dat[1] = a.dat[1]; + + c2.dat[0] = b.dat[0]; + c2.dat[1] = b.dat[1]; + + r = (*cf)(c1,c2); + + float_r.dat[0] = r.dat[0]; + float_r.dat[1] = r.dat[1]; + + return float_r; +} + +#define OPC(C,F) case C: { for(k=0;k +-- Stability : provisional +-- Portability : portable +-- +-- Display, formatting and IO functions for numeric 'Vector' and 'Matrix' +-- +----------------------------------------------------------------------------- + +module Numeric.IO ( + dispf, disps, dispcf, vecdisp, latexFormat, format, + loadMatrix, saveMatrix, fromFile, fileDimensions, + readMatrix, fromArray2D, + fscanfVector, fprintfVector, freadVector, fwriteVector +) where + +import Data.Packed +import Data.Packed.Internal +import System.Process(readProcess) +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. + +@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 + +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 + +-- | Tool to display matrices with latex syntax. +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). +-} +fileDimensions :: FilePath -> IO (Int,Int) +fileDimensions fname = do + wcres <- readProcess "wc" ["-w",fname] "" + contents <- readFile fname + let tot = read . head . words $ wcres + c = length . head . dropWhile null . map words . lines $ contents + if tot > 0 + then return (tot `div` c, c) + else return (0,0) + +-- | Loads a matrix from an ASCII file formatted as a 2D table. +loadMatrix :: FilePath -> IO (Matrix Double) +loadMatrix file = fromFile file =<< fileDimensions file + +-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance). +fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double) +fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c) + diff --git a/lib/Numeric/LinearAlgebra.hs b/lib/Numeric/LinearAlgebra.hs --- a/lib/Numeric/LinearAlgebra.hs +++ b/lib/Numeric/LinearAlgebra.hs @@ -1,7 +1,7 @@ ----------------------------------------------------------------------------- {- | Module : Numeric.LinearAlgebra -Copyright : (c) Alberto Ruiz 2006-9 +Copyright : (c) Alberto Ruiz 2006-10 License : GPL-style Maintainer : Alberto Ruiz (aruiz at um dot es) @@ -10,14 +10,19 @@ 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. + -} ----------------------------------------------------------------------------- module Numeric.LinearAlgebra ( - module Data.Packed, - module Numeric.LinearAlgebra.Algorithms, - module Numeric.LinearAlgebra.Interface + module Numeric.Container, + module Numeric.LinearAlgebra.Algorithms ) where -import Data.Packed +import Numeric.Container import Numeric.LinearAlgebra.Algorithms -import Numeric.LinearAlgebra.Interface +import Numeric.Matrix() +import Numeric.Vector() diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs @@ -1,5 +1,8 @@ -{-# OPTIONS_GHC -XFlexibleContexts -XFlexibleInstances #-} +{-# LANGUAGE FlexibleContexts, FlexibleInstances #-} {-# LANGUAGE CPP #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TypeFamilies #-} ----------------------------------------------------------------------------- {- | Module : Numeric.LinearAlgebra.Algorithms @@ -10,7 +13,7 @@ Stability : provisional Portability : uses ffi -Generic interface for the most common functions. Using it we can write higher level algorithms and testing properties for both real and complex matrices. +High level generic interface to common matrix computations. Specific functions for particular base types can also be explicitly imported from "Numeric.LinearAlgebra.LAPACK". @@ -21,9 +24,6 @@ module Numeric.LinearAlgebra.Algorithms ( -- * Supported types Field(), --- * Products - multiply, dot, - outer, kronecker, -- * Linear Systems linearSolve, luSolve, @@ -63,10 +63,9 @@ nullspaceSVD, -- * Norms Normed(..), NormType(..), + relativeError, -- * Misc - ctrans, - eps, i, - Linear(..), + eps, peps, i, -- * Util haussholder, unpackQR, unpackHess, @@ -77,18 +76,23 @@ import Data.Packed.Internal hiding ((//)) -import Data.Packed.Vector import Data.Packed.Matrix -import Data.Complex -import Numeric.GSL.Vector import Numeric.LinearAlgebra.LAPACK as LAPACK -import Numeric.LinearAlgebra.Linear import Data.List(foldl1') import Data.Array +import Numeric.ContainerBoot hiding ((.*),(*/)) --- | Auxiliary typeclass used to define generic computations for both real and complex matrices. -class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where +{- | Class used to define generic linear algebra computations for both real and complex matrices. Only double precision is supported in this version (we can +transform single precision objects using 'single' and 'double'). + +-} +class (Product t, + Convert t, + Container Vector t, + Container Matrix t, + Normed Matrix t, + Normed Vector t) => Field t where svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) sv' :: Matrix t -> Vector Double @@ -107,8 +111,6 @@ qr' :: Matrix t -> (Matrix t, Matrix t) hess' :: Matrix t -> (Matrix t, Matrix t) schur' :: Matrix t -> (Matrix t, Matrix t) - ctrans' :: Matrix t -> Matrix t - multiply' :: Matrix t -> Matrix t -> Matrix t instance Field Double where @@ -121,7 +123,6 @@ cholSolve' = cholSolveR linearSolveLS' = linearSolveLSR linearSolveSVD' = linearSolveSVDR Nothing - ctrans' = trans eig' = eigR eigSH'' = eigS eigOnly = eigOnlyR @@ -131,7 +132,6 @@ qr' = unpackQR . qrR hess' = unpackHess hessR schur' = schurR - multiply' = multiplyR instance Field (Complex Double) where #ifdef NOZGESDD @@ -148,7 +148,6 @@ cholSolve' = cholSolveC linearSolveLS' = linearSolveLSC linearSolveSVD' = linearSolveSVDC Nothing - ctrans' = conj . trans eig' = eigC eigOnly = eigOnlyC eigSH'' = eigH @@ -158,7 +157,6 @@ qr' = unpackQR . qrC hess' = unpackHess hessC schur' = schurC - multiply' = multiplyC -------------------------------------------------------------- @@ -190,7 +188,7 @@ fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t) fullSVD m = (u,d,v) where (u,s,v) = svd m - d = diagRect s r c + d = diagRect 0 s r c r = rows m c = cols m @@ -217,7 +215,7 @@ {-# DEPRECATED full "use fullSVD instead" #-} full svdFun m = (u, d ,v) where (u,s,v) = svdFun m - d = diagRect s r c + d = diagRect 0 s r c r = rows m c = cols m @@ -326,14 +324,7 @@ schur :: Field t => Matrix t -> (Matrix t, Matrix t) schur = schur' --- | Generic conjugate transpose. -ctrans :: Field t => Matrix t -> Matrix t -ctrans = ctrans' --- | Matrix product. -multiply :: Field t => Matrix t -> Matrix t -> Matrix t -multiply = {-# SCC "multiply" #-} multiply' - -- | Similar to 'cholSH', but instead of an error (e.g., caused by a matrix not positive definite) it returns 'Nothing'. mbCholSH :: Field t => Matrix t -> Maybe (Matrix t) mbCholSH = {-# SCC "mbCholSH" #-} mbCholSH' @@ -398,74 +389,15 @@ eps :: Double eps = 2.22044604925031e-16 --- | The imaginary unit: @i = 0.0 :+ 1.0@ -i :: Complex Double -i = 0:+1 - --- matrix product -mXm :: (Num t, Field t) => Matrix t -> Matrix t -> Matrix t -mXm = multiply - --- matrix - vector product -mXv :: (Num t, Field t) => Matrix t -> Vector t -> Vector t -mXv m v = flatten $ m `mXm` (asColumn v) - --- vector - matrix product -vXm :: (Num t, Field t) => Vector t -> Matrix t -> Vector t -vXm v m = flatten $ (asRow v) `mXm` m +-- | 1 + 0.5*peps == 1, 1 + 0.6*peps /= 1 +peps :: RealFloat x => x +peps = x where x = 2.0**(fromIntegral $ 1-floatDigits x) ---------------------------------------------------------------------------- - -norm2 :: Vector Double -> Double -norm2 = toScalarR Norm2 - -norm1 :: Vector Double -> Double -norm1 = toScalarR AbsSum - -data NormType = Infinity | PNorm1 | PNorm2 -- PNorm Int - -pnormRV PNorm2 = norm2 -pnormRV PNorm1 = norm1 -pnormRV Infinity = vectorMax . vectorMapR Abs ---pnormRV _ = error "pnormRV not yet defined" - -pnormCV PNorm2 = norm2 . asReal -pnormCV PNorm1 = norm1 . mapVector magnitude -pnormCV Infinity = vectorMax . mapVector magnitude ---pnormCV _ = error "pnormCV not yet defined" - -pnormRM PNorm2 m = singularValues m @> 0 -pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m -pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m) ---pnormRM _ _ = error "p norm not yet defined" - -pnormCM PNorm2 m = singularValues m @> 0 -pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m -pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m) ---pnormCM _ _ = error "p norm not yet defined" - --- | Objects which have a p-norm. --- Using it you can define convenient shortcuts: --- --- @norm2 x = pnorm PNorm2 x@ --- --- @frobenius m = norm2 . flatten $ m@ -class Normed t where - pnorm :: NormType -> t -> Double - -instance Normed (Vector Double) where - pnorm = pnormRV - -instance Normed (Vector (Complex Double)) where - pnorm = pnormCV - -instance Normed (Matrix Double) where - pnorm = pnormRM - -instance Normed (Matrix (Complex Double)) where - pnorm = pnormCM +-- | The imaginary unit: @i = 0.0 :+ 1.0@ +i :: Complex Double +i = 0:+1 ----------------------------------------------------------------------- @@ -543,7 +475,7 @@ where xs = toList v zt 0 v = v -zt k v = join [subVector 0 (dim v - k) v, constant 0 k] +zt k v = join [subVector 0 (dim v - k) v, konst 0 k] unpackQR :: (Field t) => (Matrix t, Vector t) -> (Matrix t, Matrix t) @@ -602,8 +534,8 @@ -- -- @logm = matFunc log@ -- -matFunc :: Field t => (Complex Double -> Complex Double) -> Matrix t -> Matrix (Complex Double) -matFunc f m = case diagonalize (complex m) of +matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +matFunc f m = case diagonalize m of Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v Nothing -> error "Sorry, matFunc requires a diagonalizable matrix" @@ -660,11 +592,11 @@ [ 2.0, 2.25 , 0.0, 2.0 ]@ -} -sqrtm :: Field t => Matrix t -> Matrix t +sqrtm :: Field t => Matrix t -> Matrix t sqrtm = sqrtmInv sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x)) - where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < eps = a + where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < peps = a | otherwise = fixedPoint (b:rest) fixedPoint _ = error "fixedpoint with impossible inputs" f (y,z) = (0.5 .* (y |+| inv z), @@ -697,59 +629,72 @@ c = cols l_u tu = triang r c 0 1 tl = triang r c 0 0 - l = takeColumns r (l_u |*| tl) |+| diagRect (constant 1 r) r r + l = takeColumns r (l_u |*| tl) |+| diagRect 0 (konst 1 r) r r u = l_u |*| tu (p,s) = fixPerm r perm - l' = (l_u |*| tl) |+| diagRect (constant 1 c) r c + l' = (l_u |*| tl) |+| diagRect 0 (konst 1 c) r c u' = takeRows c (l_u |*| tu) (|+|) = add (|*|) = mul --------------------------------------------------- +--------------------------------------------------------------------------- --- | Euclidean inner product. -dot :: (Field t) => Vector t -> Vector t -> t -dot u v = multiply r c @@> (0,0) - where r = asRow u - c = asColumn v +data NormType = Infinity | PNorm1 | PNorm2 | Frobenius +class (RealFloat (RealOf t)) => Normed c t where + pnorm :: NormType -> c t -> RealOf t -{- | Outer product of two vectors. +instance Normed Vector Double where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = norm2 -@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3] -(3><3) - [ 5.0, 2.0, 3.0 - , 10.0, 4.0, 6.0 - , 15.0, 6.0, 9.0 ]@ --} -outer :: (Field t) => Vector t -> Vector t -> Matrix t -outer u v = asColumn u `multiply` asRow v +instance Normed Vector (Complex Double) where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = pnorm PNorm2 -{- | Kronecker product of two matrices. +instance Normed Vector Float where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = pnorm PNorm2 -@m1=(2><3) - [ 1.0, 2.0, 0.0 - , 0.0, -1.0, 3.0 ] -m2=(4><3) - [ 1.0, 2.0, 3.0 - , 4.0, 5.0, 6.0 - , 7.0, 8.0, 9.0 - , 10.0, 11.0, 12.0 ]@ +instance Normed Vector (Complex Float) where + pnorm PNorm1 = norm1 + pnorm PNorm2 = norm2 + pnorm Infinity = normInf + pnorm Frobenius = pnorm PNorm2 -@\> kronecker m1 m2 -(8><9) - [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0 - , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0 - , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0 - , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0 - , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0 - , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0 - , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 - , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@ --} -kronecker :: (Field t) => Matrix t -> Matrix t -> Matrix t -kronecker a b = fromBlocks - . splitEvery (cols a) - . map (reshape (cols b)) - . toRows - $ flatten a `outer` flatten b + +instance Normed Matrix Double where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = (@>0) . singularValues + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +instance Normed Matrix (Complex Double) where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = (@>0) . singularValues + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +instance Normed Matrix Float where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = realToFrac . (@>0) . singularValues . double + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +instance Normed Matrix (Complex Float) where + pnorm PNorm1 = maximum . map (pnorm PNorm1) . toColumns + pnorm PNorm2 = realToFrac . (@>0) . singularValues . double + pnorm Infinity = pnorm PNorm1 . trans + pnorm Frobenius = pnorm PNorm2 . flatten + +-- | Approximate number of common digits in the maximum element. +relativeError :: (Normed c t, Container c t) => c t -> c t -> Int +relativeError x y = dig (norm (x `sub` y) / norm x) + where norm = pnorm Infinity + dig r = round $ -logBase 10 (realToFrac r :: Double) diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs deleted file mode 100644 --- a/lib/Numeric/LinearAlgebra/Instances.hs +++ /dev/null @@ -1,218 +0,0 @@ -{-# LANGUAGE UndecidableInstances, FlexibleInstances #-} ------------------------------------------------------------------------------ -{- | -Module : Numeric.LinearAlgebra.Instances -Copyright : (c) Alberto Ruiz 2006 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) -Stability : provisional -Portability : portable - -This module exports Show, Read, Eq, Num, Fractional, and Floating instances for Vector and Matrix. - -In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand. - --} ------------------------------------------------------------------------------ - -module Numeric.LinearAlgebra.Instances( -) where - -import Numeric.LinearAlgebra.Linear -import Numeric.GSL.Vector -import Data.Packed.Matrix -import Data.Complex -import Data.List(transpose,intersperse) -import Data.Packed.Internal.Vector - -#ifndef VECTOR -import Foreign(Storable) -#endif - ------------------------------------------------------------------- - -instance (Show a, Element a) => (Show (Matrix a)) where - show m = (sizes++) . dsp . map (map show) . toLists $ m - where sizes = "("++show (rows m)++"><"++show (cols m)++")\n" - -dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp - where - mt = transpose as - longs = map (maximum . map length) mt - mtp = zipWith (\a b -> map (pad a) b) longs mt - pad n str = replicate (n - length str) ' ' ++ str - unwords' = concat . intersperse ", " - -#ifndef VECTOR - -instance (Show a, Storable a) => (Show (Vector a)) where - show v = (show (dim v))++" |> " ++ show (toList v) - -#endif - ------------------------------------------------------------------- - -instance (Element a, Read a) => Read (Matrix a) where - readsPrec _ s = [((rs>' $ dims - -#ifdef VECTOR - -instance (Element a, Read a) => Read (Vector a) where - readsPrec _ s = [(fromList . read $ listnums, rest)] - where (thing,trest) = breakAt ']' s - (dims,listnums) = breakAt ' ' (dropWhile (==' ') thing) - rest = drop 31 trest -#else - -instance (Element a, Read a) => Read (Vector a) where - readsPrec _ s = [((d |>) . read $ listnums, rest)] - where (thing,rest) = breakAt ']' s - (dims,listnums) = breakAt '>' thing - d = read . init . fst . breakAt '|' $ dims - -#endif - -breakAt c l = (a++[c],tail b) where - (a,b) = break (==c) l - ------------------------------------------------------------------- - -adaptScalar f1 f2 f3 x y - | dim x == 1 = f1 (x@>0) y - | dim y == 1 = f3 x (y@>0) - | otherwise = f2 x y - -#ifndef VECTOR - -instance Linear Vector a => Eq (Vector a) where - (==) = equal - -#endif - -instance Num (Vector Double) where - (+) = adaptScalar addConstant add (flip addConstant) - negate = scale (-1) - (*) = adaptScalar scale mul (flip scale) - signum = vectorMapR Sign - abs = vectorMapR Abs - fromInteger = fromList . return . fromInteger - -instance Num (Vector (Complex Double)) where - (+) = adaptScalar addConstant add (flip addConstant) - negate = scale (-1) - (*) = adaptScalar scale mul (flip scale) - signum = vectorMapC Sign - abs = vectorMapC Abs - fromInteger = fromList . return . fromInteger - -instance Linear Matrix a => Eq (Matrix a) where - (==) = equal - -instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where - (+) = liftMatrix2Auto (+) - (-) = liftMatrix2Auto (-) - negate = liftMatrix negate - (*) = liftMatrix2Auto (*) - signum = liftMatrix signum - abs = liftMatrix abs - fromInteger = (1><1) . return . fromInteger - ---------------------------------------------------- - -instance (Linear Vector a, Num (Vector a)) => Fractional (Vector a) where - fromRational n = fromList [fromRational n] - (/) = adaptScalar f divide g where - r `f` v = scaleRecip r v - v `g` r = scale (recip r) v - -------------------------------------------------------- - -instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where - fromRational n = (1><1) [fromRational n] - (/) = liftMatrix2Auto (/) - ---------------------------------------------------------- - -instance Floating (Vector Double) where - sin = vectorMapR Sin - cos = vectorMapR Cos - tan = vectorMapR Tan - asin = vectorMapR ASin - acos = vectorMapR ACos - atan = vectorMapR ATan - sinh = vectorMapR Sinh - cosh = vectorMapR Cosh - tanh = vectorMapR Tanh - asinh = vectorMapR ASinh - acosh = vectorMapR ACosh - atanh = vectorMapR ATanh - exp = vectorMapR Exp - log = vectorMapR Log - sqrt = vectorMapR Sqrt - (**) = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS)) - pi = fromList [pi] - -------------------------------------------------------------- - -instance Floating (Vector (Complex Double)) where - sin = vectorMapC Sin - cos = vectorMapC Cos - tan = vectorMapC Tan - asin = vectorMapC ASin - acos = vectorMapC ACos - atan = vectorMapC ATan - sinh = vectorMapC Sinh - cosh = vectorMapC Cosh - tanh = vectorMapC Tanh - asinh = vectorMapC ASinh - acosh = vectorMapC ACosh - atanh = vectorMapC ATanh - exp = vectorMapC Exp - log = vectorMapC Log - sqrt = vectorMapC Sqrt - (**) = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS)) - pi = fromList [pi] - ------------------------------------------------------------ - -instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where - sin = liftMatrix sin - cos = liftMatrix cos - tan = liftMatrix tan - asin = liftMatrix asin - acos = liftMatrix acos - atan = liftMatrix atan - sinh = liftMatrix sinh - cosh = liftMatrix cosh - tanh = liftMatrix tanh - asinh = liftMatrix asinh - acosh = liftMatrix acosh - atanh = liftMatrix atanh - exp = liftMatrix exp - log = liftMatrix log - (**) = liftMatrix2Auto (**) - sqrt = liftMatrix sqrt - pi = (1><1) [pi] - ---------------------------------------------------------------- - --- instance (Storable a, Num (Vector a)) => Monoid (Vector a) where --- mempty = 0 { idim = 0 } --- mappend a b = mconcat [a,b] --- mconcat = j . filter ((>0).dim) --- where j [] = mempty --- j l = join l - ---------------------------------------------------------------- - --- instance (NFData a, Storable a) => NFData (Vector a) where --- rnf = rnf . (@>0) --- --- instance (NFData a, Element a) => NFData (Matrix a) where --- rnf = rnf . flatten - diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs deleted file mode 100644 --- a/lib/Numeric/LinearAlgebra/Interface.hs +++ /dev/null @@ -1,117 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} ------------------------------------------------------------------------------ -{- | -Module : Numeric.LinearAlgebra.Interface -Copyright : (c) Alberto Ruiz 2007 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) -Stability : provisional -Portability : portable - -Some useful operators, and Show, Read, Eq, Num, Fractional, and Floating instances for Vector and Matrix. - -In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand. - - --} ------------------------------------------------------------------------------ - -module Numeric.LinearAlgebra.Interface( - (<>),(<.>), - (<\>), - (.*),(*/), - (<|>),(<->), -) where - -import Numeric.LinearAlgebra.Instances() -import Data.Packed.Vector -import Data.Packed.Matrix -import Numeric.LinearAlgebra.Algorithms - -class Mul a b c | a b -> c where - infixl 7 <> - -- | Matrix-matrix, matrix-vector, and vector-matrix products. - (<>) :: Field t => a t -> b t -> c t - -instance Mul Matrix Matrix Matrix where - (<>) = multiply - -instance Mul Matrix Vector Vector where - (<>) m v = flatten $ m <> (asColumn v) - -instance Mul Vector Matrix Vector where - (<>) v m = flatten $ (asRow v) <> m - ---------------------------------------------------- - --- | Dot product: @u \<.\> v = dot u v@ -(<.>) :: (Field t) => Vector t -> Vector t -> t -infixl 7 <.> -(<.>) = dot - ----------------------------------------------------- - -{-# 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 -infixl 7 */ -v */ x = scale (recip x) v - --- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD). -(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a -infixl 7 <\> -m <\> v = flatten (linearSolveSVD m (reshape 1 v)) - ------------------------------------------------- - -{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-} -{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-} - -class Joinable a b where - joinH :: Element t => a t -> b t -> Matrix t - joinV :: Element t => a t -> b t -> Matrix t - -instance Joinable Matrix Matrix where - joinH m1 m2 = fromBlocks [[m1,m2]] - joinV m1 m2 = fromBlocks [[m1],[m2]] - -instance Joinable Matrix Vector where - joinH m v = joinH m (asColumn v) - joinV m v = joinV m (asRow v) - -instance Joinable Vector Matrix where - joinH v m = joinH (asColumn v) m - joinV v m = joinV (asRow v) m - -infixl 4 <|> -infixl 3 <-> - -{-- - | Horizontal concatenation of matrices and vectors: - -@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0] -(6><4) - [ 1.0, 0.0, 0.0, 1.0 - , 0.0, 1.0, 0.0, 2.0 - , 0.0, 0.0, 1.0, 3.0 - , 3.0, 0.0, 0.0, 4.0 - , 0.0, 3.0, 0.0, 5.0 - , 0.0, 0.0, 3.0, 6.0 ]@ --} --- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t -a <|> b = joinH a b - --- -- | Vertical concatenation of matrices and vectors. --- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t -a <-> b = joinV a b - diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs --- a/lib/Numeric/LinearAlgebra/LAPACK.hs +++ b/lib/Numeric/LinearAlgebra/LAPACK.hs @@ -14,7 +14,7 @@ module Numeric.LinearAlgebra.LAPACK ( -- * Matrix product - multiplyR, multiplyC, + multiplyR, multiplyC, multiplyF, multiplyQ, -- * Linear systems linearSolveR, linearSolveC, lusR, lusC, @@ -43,7 +43,8 @@ import Data.Packed.Internal import Data.Packed.Matrix -import Data.Complex +--import Data.Complex +import Numeric.Conversion import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) import Foreign import Foreign.C.Types (CInt) @@ -51,8 +52,10 @@ ----------------------------------------------------------------------------------- -foreign import ccall "LAPACK/lapack-aux.h multiplyR" dgemmc :: CInt -> CInt -> TMMM -foreign import ccall "LAPACK/lapack-aux.h multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM +foreign import ccall "multiplyR" dgemmc :: CInt -> CInt -> TMMM +foreign import ccall "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM +foreign import ccall "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM +foreign import ccall "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM isT MF{} = 0 isT MC{} = 1 @@ -69,12 +72,20 @@ -- | Matrix product based on BLAS's /dgemm/. multiplyR :: Matrix Double -> Matrix Double -> Matrix Double -multiplyR a b = multiplyAux dgemmc "dgemmc" a b +multiplyR a b = {-# SCC "multiplyR" #-} multiplyAux dgemmc "dgemmc" a b -- | Matrix product based on BLAS's /zgemm/. multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) multiplyC a b = multiplyAux zgemmc "zgemmc" a b +-- | Matrix product based on BLAS's /sgemm/. +multiplyF :: Matrix Float -> Matrix Float -> Matrix Float +multiplyF a b = multiplyAux sgemmc "sgemmc" a b + +-- | Matrix product based on BLAS's /cgemm/. +multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float) +multiplyQ a b = multiplyAux cgemmc "cgemmc" a b + ----------------------------------------------------------------------------- foreign import ccall "svd_l_R" dgesvd :: TMMVM foreign import ccall "svd_l_C" zgesvd :: TCMCMVCM @@ -248,14 +259,14 @@ where r = rows m g ra ca pa = dgeev ra ca pa 0 0 nullPtr -fixeig1 s = toComplex (subVector 0 r (asReal s), subVector r r (asReal s)) +fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s)) where r = dim s fixeig [] _ = [] -fixeig [_] [v] = [comp v] +fixeig [_] [v] = [comp' v] fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) - | r1 == r2 && i1 == (-i2) = toComplex (v1,v2) : toComplex (v1,scale (-1) v2) : fixeig r vs - | otherwise = comp v1 : fixeig ((r2:+i2):r) (v2:vs) + | r1 == r2 && i1 == (-i2) = toComplex' (v1,v2) : toComplex' (v1,scale (-1) v2) : fixeig r vs + | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) where scale = vectorMapValR Scale fixeig _ _ = error "fixeig with impossible inputs" diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c --- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c +++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c @@ -11,15 +11,25 @@ #define MIN(A,B) ((A)<(B)?(A):(B)) #define MAX(A,B) ((A)>(B)?(A):(B)) - + +// #define DBGL + #ifdef DBGL -#define DEBUGMSG(M) printf("LAPACK Wrapper "M"\n: "); size_t t0 = time(NULL); -#define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;); +#define DEBUGMSG(M) printf("\nLAPACK "M"\n"); #else #define DEBUGMSG(M) -#define OK return 0; #endif +#define OK return 0; + +// #ifdef DBGL +// #define DEBUGMSG(M) printf("LAPACK Wrapper "M"\n: "); size_t t0 = time(NULL); +// #define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;); +// #else +// #define DEBUGMSG(M) +// #define OK return 0; +// #endif + #define TRACEMAT(M) {int q; printf(" %d x %d: ",M##r,M##c); \ for(q=0;q=1 && ar==ac && ar==br,BAD_SIZE); DEBUGMSG("linearSolveC_l"); - double*AC = (double*)malloc(2*n*n*sizeof(double)); - memcpy(AC,ap,2*n*n*sizeof(double)); - memcpy(xp,bp,2*n*nhrs*sizeof(double)); + doublecomplex*AC = (doublecomplex*)malloc(n*n*sizeof(doublecomplex)); + memcpy(AC,ap,n*n*sizeof(doublecomplex)); + memcpy(xp,bp,n*nhrs*sizeof(doublecomplex)); integer * ipiv = (integer*)malloc(n*sizeof(integer)); integer res; zgesv_ (&n,&nhrs, - (doublecomplex*)AC, &n, + AC, &n, ipiv, - (doublecomplex*)xp, &n, + xp, &n, &res); if(res>0) { return SINGULAR; @@ -517,12 +527,12 @@ integer nhrs = bc; REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); DEBUGMSG("cholSolveC_l"); - memcpy(xp,bp,2*n*nhrs*sizeof(double)); + memcpy(xp,bp,n*nhrs*sizeof(doublecomplex)); integer res; zpotrs_ ("U", &n,&nhrs, (doublecomplex*)ap, &n, - (doublecomplex*)xp, &n, + xp, &n, &res); CHECK(res,res); OK @@ -581,31 +591,30 @@ integer ldb = xr; REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); DEBUGMSG("linearSolveLSC_l"); - double*AC = (double*)malloc(2*m*n*sizeof(double)); - memcpy(AC,ap,2*m*n*sizeof(double)); - memcpy(AC,ap,2*m*n*sizeof(double)); + doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); + memcpy(AC,ap,m*n*sizeof(doublecomplex)); if (m>=n) { - memcpy(xp,bp,2*m*nrhs*sizeof(double)); + memcpy(xp,bp,m*nrhs*sizeof(doublecomplex)); } else { int k; for(k = 0; k0) { @@ -685,16 +694,16 @@ integer ldb = xr; REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); DEBUGMSG("linearSolveSVDC_l"); - double*AC = (double*)malloc(2*m*n*sizeof(double)); + doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); double*S = (double*)malloc(MIN(m,n)*sizeof(double)); double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double)); - memcpy(AC,ap,2*m*n*sizeof(double)); + memcpy(AC,ap,m*n*sizeof(doublecomplex)); if (m>=n) { - memcpy(xp,bp,2*m*nrhs*sizeof(double)); + memcpy(xp,bp,m*nrhs*sizeof(doublecomplex)); } else { int k; for(k = 0; k0,NODEFPOS); CHECK(res,res); doublecomplex zero = {0.,0.}; int r,c; for (r=0; r0) { return NOCONVER; @@ -940,7 +949,7 @@ integer* auxipiv = (integer*)malloc(mn*sizeof(integer)); memcpy(rp,ap,m*n*sizeof(doublecomplex)); integer res; - zgetrf_ (&m,&n,(doublecomplex*)rp,&m,auxipiv,&res); + zgetrf_ (&m,&n,rp,&m,auxipiv,&res); if(res>0) { res = 0; // fixme } @@ -990,7 +999,7 @@ } integer res; memcpy(xp,bp,mrhs*nrhs*sizeof(doublecomplex)); - zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&m,auxipiv,(doublecomplex*)xp,&mrhs,&res); + zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&m,auxipiv,xp,&mrhs,&res); CHECK(res,res); free(auxipiv); OK @@ -1004,6 +1013,7 @@ int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) { //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); + DEBUGMSG("dgemm_"); integer m = ta?ac:ar; integer n = tb?br:bc; integer k = ta?ar:ac; @@ -1022,6 +1032,7 @@ int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) { //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); + DEBUGMSG("zgemm_"); integer m = ta?ac:ar; integer n = tb?br:bc; integer k = ta?ar:ac; @@ -1031,14 +1042,67 @@ doublecomplex alpha = {1,0}; doublecomplex beta = {0,0}; zgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha, - (doublecomplex*)ap,&lda, - (doublecomplex*)bp,&ldb,&beta, - (doublecomplex*)rp,&ldc); + ap,&lda, + bp,&ldb,&beta, + rp,&ldc); OK } +void sgemm_(char *, char *, integer *, integer *, integer *, + float *, const float *, integer *, const float *, + integer *, float *, float *, integer *); + +int multiplyF(int ta, int tb, KFMAT(a),KFMAT(b),FMAT(r)) { + //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); + DEBUGMSG("sgemm_"); + integer m = ta?ac:ar; + integer n = tb?br:bc; + integer k = ta?ar:ac; + integer lda = ar; + integer ldb = br; + integer ldc = rr; + float alpha = 1; + float beta = 0; + sgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc); + OK +} + +void cgemm_(char *, char *, integer *, integer *, integer *, + complex *, const complex *, integer *, const complex *, + integer *, complex *, complex *, integer *); + +int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) { + //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); + DEBUGMSG("cgemm_"); + integer m = ta?ac:ar; + integer n = tb?br:bc; + integer k = ta?ar:ac; + integer lda = ar; + integer ldb = br; + integer ldc = rr; + complex alpha = {1,0}; + complex beta = {0,0}; + cgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha, + ap,&lda, + bp,&ldb,&beta, + rp,&ldc); + OK +} + //////////////////// transpose ///////////////////////// +int transF(KFMAT(x),FMAT(t)) { + REQUIRES(xr==tc && xc==tr,BAD_SIZE); + DEBUGMSG("transF"); + int i,j; + for (i=0; i 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 - -- | element by element multiplication - mul :: c e -> c e -> c e - -- | element by element division - divide :: c e -> c e -> c e - equal :: c e -> c e -> Bool - - -instance Linear Vector Double where - scale = vectorMapValR Scale - scaleRecip = vectorMapValR Recip - addConstant = vectorMapValR AddConstant - add = vectorZipR Add - sub = vectorZipR Sub - 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 - scaleRecip = vectorMapValC Recip - addConstant = vectorMapValC AddConstant - add = vectorZipC Add - sub = vectorZipC Sub - 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) - scaleRecip x = liftMatrix (scaleRecip x) - addConstant x = liftMatrix (addConstant x) - add = liftMatrix2 add - sub = liftMatrix2 sub - mul = liftMatrix2 mul - divide = liftMatrix2 divide - equal a b = cols a == cols b && flatten a `equal` flatten b - scalar x = (1><1) [x] diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ b/lib/Numeric/LinearAlgebra/Tests.hs @@ -21,11 +21,12 @@ --, runBigTests ) where +import Data.Packed.Random import Numeric.LinearAlgebra import Numeric.LinearAlgebra.LAPACK import Numeric.LinearAlgebra.Tests.Instances import Numeric.LinearAlgebra.Tests.Properties -import Test.HUnit hiding ((~:),test,Testable) +import Test.HUnit hiding ((~:),test,Testable,State) import System.Info import Data.List(foldl1') import Numeric.GSL @@ -33,9 +34,14 @@ import qualified Prelude import System.CPUTime import Text.Printf +import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr) +import Control.Arrow((***)) +import Debug.Trace #include "Tests/quickCheckCompat.h" +debug x = trace (show x) x + a ^ b = a Prelude.^ (b :: Int) utest str b = TestCase $ assertBool str b @@ -44,6 +50,8 @@ feye n = flipud (ident n) :: Matrix Double +----------------------------------------------------------- + detTest1 = det m == 26 && det mc == 38 :+ (-3) && det (feye 2) == -1 @@ -164,7 +172,7 @@ sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d - ok2 = pnorm PNorm2 (fromList (map fst sols) - fromList sol) < 1E-5 + ok2 = norm2 (fromList (map fst sols) - fromList sol) < 1E-5 ----------------------------------------------------- @@ -201,6 +209,155 @@ where fun n = foldl1' (<>) (map rot angles) where angles = toList $ linspace n (0,1) +--------------------------------------------------------------------- +-- vector <= 0.6.0.2 bug discovered by Patrick Perry +-- http://trac.haskell.org/vector/ticket/31 + +offsetTest = y == y' where + x = fromList [0..3 :: Double] + y = subVector 1 3 x + (f,o,n) = unsafeToForeignPtr y + y' = unsafeFromForeignPtr f o n + +--------------------------------------------------------------------- + +normsVTest = TestList [ + utest "normv2CD" $ norm2PropC v + , utest "normv2CF" $ norm2PropC (single v) + , utest "normv2D" $ norm2PropR x + , utest "normv2F" $ norm2PropR (single x) + + , utest "normv1CD" $ norm1 v == 8 + , utest "normv1CF" $ norm1 (single v) == 8 + , utest "normv1D" $ norm1 x == 6 + , utest "normv1F" $ norm1 (single x) == 6 + + , utest "normvInfCD" $ normInf v == 5 + , utest "normvInfCF" $ normInf (single v) == 5 + , utest "normvInfD" $ normInf x == 3 + , utest "normvInfF" $ normInf (single x) == 3 + + ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double) + x = fromList [1,2,-3] :: Vector Double + norm2PropR a = norm2 a =~= sqrt (dot a a) + norm2PropC a = norm2 a =~= realPart (sqrt (dot a (conj a))) + a =~= b = fromList [a] |~| fromList [b] + +normsMTest = TestList [ + utest "norm2mCD" $ pnorm PNorm2 v =~= 8.86164970498005 + , utest "norm2mCF" $ pnorm PNorm2 (single v) =~= 8.86164970498005 + , utest "norm2mD" $ pnorm PNorm2 x =~= 5.96667765076216 + , utest "norm2mF" $ pnorm PNorm2 (single x) =~= 5.96667765076216 + + , utest "norm1mCD" $ pnorm PNorm1 v == 9 + , utest "norm1mCF" $ pnorm PNorm1 (single v) == 9 + , utest "norm1mD" $ pnorm PNorm1 x == 7 + , utest "norm1mF" $ pnorm PNorm1 (single x) == 7 + + , utest "normmInfCD" $ pnorm Infinity v == 12 + , utest "normmInfCF" $ pnorm Infinity (single v) == 12 + , utest "normmInfD" $ pnorm Infinity x == 8 + , utest "normmInfF" $ pnorm Infinity (single x) == 8 + + , utest "normmFroCD" $ pnorm Frobenius v =~= 8.88819441731559 + , utest "normmFroCF" $ pnorm Frobenius (single v) =~~= 8.88819441731559 + , utest "normmFroD" $ pnorm Frobenius x =~= 6.24499799839840 + , utest "normmFroF" $ pnorm Frobenius (single x) =~~= 6.24499799839840 + + ] where v = (2><2) [1,-2*i,3:+4,7] :: Matrix (Complex Double) + x = (2><2) [1,2,-3,5] :: Matrix Double + a =~= b = fromList [a] :~10~: fromList [b] + a =~~= b = fromList [a] :~5~: fromList [b] + +--------------------------------------------------------------------- + +sumprodTest = TestList [ + utest "sumCD" $ sumElements z == 6 + , utest "sumCF" $ sumElements (single z) == 6 + , utest "sumD" $ sumElements v == 6 + , utest "sumF" $ sumElements (single v) == 6 + + , utest "prodCD" $ prodProp z + , utest "prodCF" $ prodProp (single z) + , utest "prodD" $ prodProp v + , utest "prodF" $ prodProp (single v) + ] where v = fromList [1,2,3] :: Vector Double + z = fromList [1,2-i,3+i] + prodProp x = prodElements x == product (toList x) + +--------------------------------------------------------------------- + +chainTest = utest "chain" $ foldl1' (<>) ms |~| optimiseMult ms where + ms = [ diag (fromList [1,2,3 :: Double]) + , konst 3 (3,5) + , (5><10) [1 .. ] + , konst 5 (10,2) + ] + +--------------------------------------------------------------------- + +conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m) + +--------------------------------------------------------------------- + +newtype State s a = State { runState :: s -> (a,s) } + +instance Monad (State s) where + return a = State $ \s -> (a,s) + m >>= f = State $ \s -> let (a,s') = runState m s + in runState (f a) s' + +state_get :: State s s +state_get = State $ \s -> (s,s) + +state_put :: s -> State s () +state_put s = State $ \_ -> ((),s) + +evalState :: State s a -> s -> a +evalState m s = let (a,s') = runState m s + in seq s' a + +newtype MaybeT m a = MaybeT { runMaybeT :: m (Maybe a) } + +instance Monad m => Monad (MaybeT m) where + return a = MaybeT $ return $ Just a + m >>= f = MaybeT $ do + res <- runMaybeT m + case res of + Nothing -> return Nothing + Just r -> runMaybeT (f r) + fail _ = MaybeT $ return Nothing + +lift_maybe m = MaybeT $ do + res <- m + return $ Just res + +-- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs +--successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool +successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (dim v - 1) v))) (v @> 0) + where step e = do + ep <- lift_maybe $ state_get + if t e ep + then lift_maybe $ state_put e + else (fail "successive_ test failed") + +-- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input +--successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b +successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0) + where step e = do + ep <- state_get + state_put e + return $ f ep e + + +succTest = utest "successive" $ + successive_ (>) (fromList [1 :: Double,2,3,4]) == True + && successive_ (>) (fromList [1 :: Double,3,2,4]) == False + && successive (+) (fromList [1..10 :: Double]) == 9 |> [3,5,7,9,11,13,15,17,19] + +--------------------------------------------------------------------- + + -- | All tests must pass with a maximum dimension of about 20 -- (some tests may fail with bigger sizes due to precision loss). runTests :: Int -- ^ maximum dimension @@ -208,14 +365,22 @@ runTests n = do setErrorHandlerOff let test p = qCheck n p - putStrLn "------ mult" - test (multProp1 . rConsist) - test (multProp1 . cConsist) - test (multProp2 . rConsist) - test (multProp2 . cConsist) + putStrLn "------ mult Double" + test (multProp1 10 . rConsist) + test (multProp1 10 . cConsist) + test (multProp2 10 . rConsist) + test (multProp2 10 . cConsist) + putStrLn "------ mult Float" + test (multProp1 6 . (single *** single) . rConsist) + test (multProp1 6 . (single *** single) . cConsist) + test (multProp2 6 . (single *** single) . rConsist) + test (multProp2 6 . (single *** single) . cConsist) putStrLn "------ sub-trans" test (subProp . rM) test (subProp . cM) + putStrLn "------ ctrans" + test (conjuTest . cM) + test (conjuTest . zM) putStrLn "------ lu" test (luProp . rM) test (luProp . cM) @@ -286,6 +451,9 @@ test (qrProp . cM) test (rqProp . rM) test (rqProp . cM) + test (rqProp1 . cM) + test (rqProp2 . cM) + test (rqProp3 . cM) putStrLn "------ hess" test (hessProp . rSq) test (hessProp . cSq) @@ -296,21 +464,31 @@ test (cholProp . rPosDef) test (cholProp . cPosDef) putStrLn "------ expm" - test (expmDiagProp . rSqWC) + test (expmDiagProp . complex. rSqWC) test (expmDiagProp . cSqWC) putStrLn "------ fft" test (\v -> ifft (fft v) |~| v) - putStrLn "------ vector operations" + putStrLn "------ vector operations - Double" test (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::RM)) test $ (\u -> sin u ^ 2 + cos u ^ 2 |~| (1::CM)) . liftMatrix makeUnitary test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM)) test (\u -> cos u * tan u |~| sin (u::RM)) test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary + putStrLn "------ vector operations - Float" + test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM)) + test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary + test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM)) + test (\u -> cos u * tan u |~~| sin (u::FM)) + test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary putStrLn "------ read . show" test (\m -> (m::RM) == read (show m)) test (\m -> (m::CM) == read (show m)) test (\m -> toRows (m::RM) == read (show (toRows m))) test (\m -> toRows (m::CM) == read (show (toRows m))) + test (\m -> (m::FM) == read (show m)) + test (\m -> (m::ZM) == read (show m)) + test (\m -> toRows (m::FM) == read (show (toRows m))) + test (\m -> toRows (m::ZM) == read (show (toRows m))) putStrLn "------ some unit tests" _ <- runTestTT $ TestList [ utest "1E5 rots" rotTest @@ -332,7 +510,7 @@ , utest "randomGaussian" randomTestGaussian , utest "randomUniform" randomTestUniform , utest "buildVector/Matrix" $ - comp (10 |> [0::Double ..]) == buildVector 10 fromIntegral + complex (10 |> [0::Double ..]) == buildVector 10 fromIntegral && 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 @@ -340,9 +518,20 @@ , odeTest , fittingTest , mbCholTest + , utest "offset" offsetTest + , normsVTest + , normsMTest + , sumprodTest + , chainTest + , succTest ] return () + +-- single precision approximate equality +infixl 4 |~~| +a |~~| b = a :~6~: b + makeUnitary v | realPart n > 1 = v / scalar n | otherwise = v where n = sqrt (conj v <.> v) @@ -356,6 +545,7 @@ -- | Performance measurements. runBenchmarks :: IO () runBenchmarks = do + --cholBench solveBench subBench multBench @@ -455,3 +645,18 @@ solveBenchN 500 solveBenchN 1000 -- solveBenchN 1500 + +-------------------------------- + +cholBenchN n = do + let x = uniformSample 777 (2*n) (replicate n (-1,1)) + a = trans x <> x + a `seq` putStrLn "" + time ("chol " ++ show n) (chol a) + +cholBench = do + cholBenchN 1200 + cholBenchN 600 + cholBenchN 300 +-- cholBenchN 150 +-- cholBenchN 50 diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs --- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP #-} +{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-} {-# OPTIONS_GHC -fno-warn-unused-imports #-} ----------------------------------------------------------------------------- {- | @@ -22,17 +22,16 @@ SqWC(..), rSqWC, cSqWC, PosDef(..), rPosDef, cPosDef, Consistent(..), rConsist, cConsist, - RM,CM, rM,cM + RM,CM, rM,cM, + FM,ZM, fM,zM ) where - - +import System.Random import Numeric.LinearAlgebra import Control.Monad(replicateM) #include "quickCheckCompat.h" - #if MIN_VERSION_QuickCheck(2,0,0) shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]] shrinkListElementwise [] = [] @@ -43,7 +42,8 @@ shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ] #endif - +#if MIN_VERSION_QuickCheck(2,1,1) +#else instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where arbitrary = do re <- arbitrary @@ -58,6 +58,8 @@ coarbitrary = undefined #endif +#endif + chooseDim = sized $ \m -> choose (1,max 1 m) instance (Field a, Arbitrary a) => Arbitrary (Vector a) where @@ -68,7 +70,7 @@ #if MIN_VERSION_QuickCheck(2,0,0) -- shrink any one of the components shrink = map fromList . shrinkListElementwise . toList - + #else coarbitrary = undefined #endif @@ -133,10 +135,14 @@ coarbitrary = undefined #endif +class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a +instance ArbitraryField Double +instance ArbitraryField (Complex Double) + -- a well-conditioned general matrix (the singular values are between 1 and 100) newtype (WC a) = WC (Matrix a) deriving Show -instance (Field a, Arbitrary a) => Arbitrary (WC a) where +instance (ArbitraryField a) => Arbitrary (WC a) where arbitrary = do m <- arbitrary let (u,_,v) = svd m @@ -144,7 +150,7 @@ c = cols m n = min r c sv' <- replicateM n (choose (1,100)) - let s = diagRect (fromList sv') r c + let s = diagRect 0 (fromList sv') r c return $ WC (u <> real s <> trans v) #if MIN_VERSION_QuickCheck(2,0,0) @@ -155,7 +161,7 @@ -- a well-conditioned square matrix (the singular values are between 1 and 100) newtype (SqWC a) = SqWC (Matrix a) deriving Show -instance (Field a, Arbitrary a) => Arbitrary (SqWC a) where +instance (ArbitraryField a) => Arbitrary (SqWC a) where arbitrary = do Sq m <- arbitrary let (u,_,v) = svd m @@ -172,7 +178,8 @@ -- a positive definite square matrix (the eigenvalues are between 0 and 100) newtype (PosDef a) = PosDef (Matrix a) deriving Show -instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (PosDef a) where +instance (ArbitraryField a, Num (Vector a)) + => Arbitrary (PosDef a) where arbitrary = do Her m <- arbitrary let (_,v) = eigSH m @@ -209,10 +216,16 @@ type RM = Matrix Double type CM = Matrix (Complex Double) +type FM = Matrix Float +type ZM = Matrix (Complex Float) + rM m = m :: RM cM m = m :: CM +fM m = m :: FM +zM m = m :: ZM + rHer (Her m) = m :: RM cHer (Her m) = m :: CM @@ -233,3 +246,4 @@ rConsist (Consistent (a,b)) = (a,b::RM) cConsist (Consistent (a,b)) = (a,b::CM) + diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs --- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs @@ -32,7 +32,7 @@ svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, eigProp, eigSHProp, eigProp2, eigSHProp2, - qrProp, rqProp, + qrProp, rqProp, rqProp1, rqProp2, rqProp3, hessProp, schurProp1, schurProp2, cholProp, @@ -42,24 +42,27 @@ linearSolveProp, linearSolveProp2 ) where -import Numeric.LinearAlgebra +import Numeric.LinearAlgebra --hiding (real,complex) import Numeric.LinearAlgebra.LAPACK import Debug.Trace #include "quickCheckCompat.h" +--real x = real'' x +--complex x = complex'' x + debug x = trace (show x) x -- relative error -dist :: (Normed t, Num t) => t -> t -> Double -dist a b = r +dist :: (Normed c t, Num (c t)) => c t -> c t -> Double +dist a b = realToFrac r where norm = pnorm Infinity na = norm a nb = norm b nab = norm (a-b) mx = max na nb mn = min na nb - r = if mn < eps + r = if mn < peps then mx else nab/mx @@ -68,7 +71,7 @@ --a |~| b = dist a b < 10^^(-10) data Aprox a = (:~) a Int -(~:) :: (Normed a, Num a) => Aprox a -> a -> Bool +-- (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool a :~n~: b = dist a b < 10^^(-n) ------------------------------------------------------ @@ -135,7 +138,7 @@ svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where (u,s,v) = svdfun m - d = diagRect s (rows m) (cols m) + d = diagRect 0 s (rows m) (cols m) svdProp1b svdfun m = unitary u && unitary v where (u,_,v) = svdfun m @@ -207,16 +210,22 @@ qrProp m = q <> r |~| m && unitary q && upperTriang r where (q,r) = qr m -rqProp m = r <> q |~| m && unitary q && utr +rqProp m = r <> q |~| m && unitary q && upperTriang' r where (r,q) = rq m - upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1 - where t = f-c - utr = upptr (rows r) (cols r) * r |~| r -upperTriang' m = rows m == 1 || down |~| z - where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m)) - z = constant 0 (dim down) +rqProp1 m = r <> q |~| m + where (r,q) = rq m +rqProp2 m = unitary q + where (r,q) = rq m + +rqProp3 m = upperTriang' r + where (r,q) = rq m + +upperTriang' r = upptr (rows r) (cols r) * r |~| r + where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1 + where t = f-c + hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h where (p,h) = hess m @@ -237,9 +246,9 @@ mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ] where doth u v = sum $ zipWith (*) (toList u) (toList v) -multProp1 (a,b) = a <> b |~| mulH a b +multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) -multProp2 (a,b) = ctrans (a <> b) |~| ctrans b <> ctrans a +multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a) linearSolveProp f m = f m m |~| ident (rows m) diff --git a/lib/Numeric/Matrix.hs b/lib/Numeric/Matrix.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Matrix.hs @@ -0,0 +1,72 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Matrix +-- Copyright : (c) Alberto Ruiz 2010 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Provides instances of standard classes 'Show', 'Read', 'Eq', +-- 'Num', 'Fractional', and 'Floating' for 'Matrix'. +-- +-- In arithmetic operations one-component +-- vectors and matrices automatically expand to match the dimensions of the other operand. + +----------------------------------------------------------------------------- + +module Numeric.Matrix ( + ) where + +------------------------------------------------------------------- + +import Numeric.Container + +------------------------------------------------------------------- + +instance Container Matrix a => Eq (Matrix a) where + (==) = equal + +instance (Container Matrix a, Num (Vector a)) => Num (Matrix a) where + (+) = liftMatrix2Auto (+) + (-) = liftMatrix2Auto (-) + negate = liftMatrix negate + (*) = liftMatrix2Auto (*) + signum = liftMatrix signum + abs = liftMatrix abs + fromInteger = (1><1) . return . fromInteger + +--------------------------------------------------- + +instance (Container Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where + fromRational n = (1><1) [fromRational n] + (/) = liftMatrix2Auto (/) + +--------------------------------------------------------- + +instance (Floating a, Container Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where + sin = liftMatrix sin + cos = liftMatrix cos + tan = liftMatrix tan + asin = liftMatrix asin + acos = liftMatrix acos + atan = liftMatrix atan + sinh = liftMatrix sinh + cosh = liftMatrix cosh + tanh = liftMatrix tanh + asinh = liftMatrix asinh + acosh = liftMatrix acosh + atanh = liftMatrix atanh + exp = liftMatrix exp + log = liftMatrix log + (**) = liftMatrix2Auto (**) + sqrt = liftMatrix sqrt + pi = (1><1) [pi] diff --git a/lib/Numeric/Vector.hs b/lib/Numeric/Vector.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Vector.hs @@ -0,0 +1,223 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Vector +-- Copyright : (c) Alberto Ruiz 2010 +-- License : GPL-style +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- Portability : portable +-- +-- Provides instances of standard classes 'Show', 'Read', 'Eq', +-- 'Num', 'Fractional', and 'Floating' for 'Vector'. +-- +----------------------------------------------------------------------------- + +module Numeric.Vector ( + ) where + +import Numeric.GSL.Vector +import Numeric.Container + +------------------------------------------------------------------- + +#ifndef VECTOR +import Foreign(Storable) +#endif + +------------------------------------------------------------------ + +#ifndef VECTOR + +instance (Show a, Storable a) => (Show (Vector a)) where + show v = (show (dim v))++" |> " ++ show (toList v) + +#endif + +#ifdef VECTOR + +instance (Element a, Read a) => Read (Vector a) where + readsPrec _ s = [(fromList . read $ listnums, rest)] + where (thing,trest) = breakAt ']' s + (dims,listnums) = breakAt ' ' (dropWhile (==' ') thing) + rest = drop 31 trest +#else + +instance (Element a, Read a) => Read (Vector a) where + readsPrec _ s = [((d |>) . read $ listnums, rest)] + where (thing,rest) = breakAt ']' s + (dims,listnums) = breakAt '>' thing + d = read . init . fst . breakAt '|' $ dims + +#endif + +breakAt c l = (a++[c],tail b) where + (a,b) = break (==c) l + + +------------------------------------------------------------------ + +adaptScalar f1 f2 f3 x y + | dim x == 1 = f1 (x@>0) y + | dim y == 1 = f3 x (y@>0) + | otherwise = f2 x y + +------------------------------------------------------------------ + +#ifndef VECTOR + +instance Container Vector a => Eq (Vector a) where + (==) = equal + +#endif + +instance Num (Vector Float) where + (+) = adaptScalar addConstant add (flip addConstant) + negate = scale (-1) + (*) = adaptScalar scale mul (flip scale) + signum = vectorMapF Sign + abs = vectorMapF Abs + fromInteger = fromList . return . fromInteger + +instance Num (Vector Double) where + (+) = adaptScalar addConstant add (flip addConstant) + negate = scale (-1) + (*) = adaptScalar scale mul (flip scale) + signum = vectorMapR Sign + abs = vectorMapR Abs + fromInteger = fromList . return . fromInteger + +instance Num (Vector (Complex Double)) where + (+) = adaptScalar addConstant add (flip addConstant) + negate = scale (-1) + (*) = adaptScalar scale mul (flip scale) + signum = vectorMapC Sign + abs = vectorMapC Abs + fromInteger = fromList . return . fromInteger + +instance Num (Vector (Complex Float)) where + (+) = adaptScalar addConstant add (flip addConstant) + negate = scale (-1) + (*) = adaptScalar scale mul (flip scale) + signum = vectorMapQ Sign + abs = vectorMapQ Abs + fromInteger = fromList . return . fromInteger + +--------------------------------------------------- + +instance (Container Vector a, Num (Vector a)) => Fractional (Vector a) where + fromRational n = fromList [fromRational n] + (/) = adaptScalar f divide g where + r `f` v = scaleRecip r v + v `g` r = scale (recip r) v + +------------------------------------------------------- + +instance Floating (Vector Float) where + sin = vectorMapF Sin + cos = vectorMapF Cos + tan = vectorMapF Tan + asin = vectorMapF ASin + acos = vectorMapF ACos + atan = vectorMapF ATan + sinh = vectorMapF Sinh + cosh = vectorMapF Cosh + tanh = vectorMapF Tanh + asinh = vectorMapF ASinh + acosh = vectorMapF ACosh + atanh = vectorMapF ATanh + exp = vectorMapF Exp + log = vectorMapF Log + sqrt = vectorMapF Sqrt + (**) = adaptScalar (vectorMapValF PowSV) (vectorZipF Pow) (flip (vectorMapValF PowVS)) + pi = fromList [pi] + +------------------------------------------------------------- + +instance Floating (Vector Double) where + sin = vectorMapR Sin + cos = vectorMapR Cos + tan = vectorMapR Tan + asin = vectorMapR ASin + acos = vectorMapR ACos + atan = vectorMapR ATan + sinh = vectorMapR Sinh + cosh = vectorMapR Cosh + tanh = vectorMapR Tanh + asinh = vectorMapR ASinh + acosh = vectorMapR ACosh + atanh = vectorMapR ATanh + exp = vectorMapR Exp + log = vectorMapR Log + sqrt = vectorMapR Sqrt + (**) = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS)) + pi = fromList [pi] + +------------------------------------------------------------- + +instance Floating (Vector (Complex Double)) where + sin = vectorMapC Sin + cos = vectorMapC Cos + tan = vectorMapC Tan + asin = vectorMapC ASin + acos = vectorMapC ACos + atan = vectorMapC ATan + sinh = vectorMapC Sinh + cosh = vectorMapC Cosh + tanh = vectorMapC Tanh + asinh = vectorMapC ASinh + acosh = vectorMapC ACosh + atanh = vectorMapC ATanh + exp = vectorMapC Exp + log = vectorMapC Log + sqrt = vectorMapC Sqrt + (**) = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS)) + pi = fromList [pi] + +----------------------------------------------------------- + +instance Floating (Vector (Complex Float)) where + sin = vectorMapQ Sin + cos = vectorMapQ Cos + tan = vectorMapQ Tan + asin = vectorMapQ ASin + acos = vectorMapQ ACos + atan = vectorMapQ ATan + sinh = vectorMapQ Sinh + cosh = vectorMapQ Cosh + tanh = vectorMapQ Tanh + asinh = vectorMapQ ASinh + acosh = vectorMapQ ACosh + atanh = vectorMapQ ATanh + exp = vectorMapQ Exp + log = vectorMapQ Log + sqrt = vectorMapQ Sqrt + (**) = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS)) + pi = fromList [pi] + +----------------------------------------------------------- + + +-- instance (Storable a, Num (Vector a)) => Monoid (Vector a) where +-- mempty = 0 { idim = 0 } +-- mappend a b = mconcat [a,b] +-- mconcat = j . filter ((>0).dim) +-- where j [] = mempty +-- j l = join l + +--------------------------------------------------------------- + +-- instance (NFData a, Storable a) => NFData (Vector a) where +-- rnf = rnf . (@>0) +-- +-- instance (NFData a, Element a) => NFData (Matrix a) where +-- rnf = rnf . flatten + + +-------------------------------------------------------------------------- +