From 0ab93b8eb934167e16dbae193c4fd2b5359a184b Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 16 May 2014 09:20:47 +0200 Subject: instances moved to base --- packages/base/hmatrix-base.cabal | 3 + packages/base/src/Numeric/Chain.hs | 143 +++++++++++++++++++++++++++ packages/base/src/Numeric/Matrix.hs | 100 +++++++++++++++++++ packages/base/src/Numeric/Vector.hs | 159 ++++++++++++++++++++++++++++++ packages/hmatrix/Config.hs | 154 ----------------------------- packages/hmatrix/Setup.lhs | 15 +-- packages/hmatrix/hmatrix.cabal | 3 - packages/hmatrix/src/Numeric/Chain.hs | 140 -------------------------- packages/hmatrix/src/Numeric/Container.hs | 6 +- packages/hmatrix/src/Numeric/Matrix.hs | 98 ------------------ packages/hmatrix/src/Numeric/Vector.hs | 158 ----------------------------- 11 files changed, 411 insertions(+), 568 deletions(-) create mode 100644 packages/base/src/Numeric/Chain.hs create mode 100644 packages/base/src/Numeric/Matrix.hs create mode 100644 packages/base/src/Numeric/Vector.hs delete mode 100644 packages/hmatrix/Config.hs delete mode 100644 packages/hmatrix/src/Numeric/Chain.hs delete mode 100644 packages/hmatrix/src/Numeric/Matrix.hs delete mode 100644 packages/hmatrix/src/Numeric/Vector.hs diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal index 8571e7a..8eeb97e 100644 --- a/packages/base/hmatrix-base.cabal +++ b/packages/base/hmatrix-base.cabal @@ -39,12 +39,15 @@ library Numeric.LinearAlgebra.LAPACK Data.Packed.Numeric Numeric.Vectorized + Numeric.Vector + Numeric.Matrix other-modules: Data.Packed.Internal, Data.Packed.Internal.Common, Data.Packed.Internal.Signatures, Data.Packed.Internal.Vector, Data.Packed.Internal.Matrix + Numeric.Chain C-sources: src/C/lapack-aux.c src/C/vector-aux.c diff --git a/packages/base/src/Numeric/Chain.hs b/packages/base/src/Numeric/Chain.hs new file mode 100644 index 0000000..fbdb01b --- /dev/null +++ b/packages/base/src/Numeric/Chain.hs @@ -0,0 +1,143 @@ +----------------------------------------------------------------------------- +-- | +-- 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 Data.Packed.Numeric + +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/packages/base/src/Numeric/Matrix.hs b/packages/base/src/Numeric/Matrix.hs new file mode 100644 index 0000000..3478aae --- /dev/null +++ b/packages/base/src/Numeric/Matrix.hs @@ -0,0 +1,100 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} + +----------------------------------------------------------------------------- +-- | +-- 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 Data.Packed +import Data.Packed.Numeric +import qualified Data.Monoid as M +import Data.List(partition) +import Numeric.Chain + +------------------------------------------------------------------- + +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] + +-------------------------------------------------------------------------------- + +isScalar m = rows m == 1 && cols m == 1 + +adaptScalarM f1 f2 f3 x y + | isScalar x = f1 (x @@>(0,0) ) y + | isScalar y = f3 x (y @@>(0,0) ) + | otherwise = f2 x y + +instance (Container Vector t, Eq t, Num (Vector t), Product t) => M.Monoid (Matrix t) + where + mempty = 1 + mappend = adaptScalarM scale mXm (flip scale) + + mconcat xs = work (partition isScalar xs) + where + work (ss,[]) = product ss + work (ss,ms) = scale' (product ss) (optimiseMult ms) + scale' x m + | isScalar x && x00 == 1 = m + | otherwise = scale x00 m + where + x00 = x @@> (0,0) + diff --git a/packages/base/src/Numeric/Vector.hs b/packages/base/src/Numeric/Vector.hs new file mode 100644 index 0000000..2769cd9 --- /dev/null +++ b/packages/base/src/Numeric/Vector.hs @@ -0,0 +1,159 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Vector +-- Copyright : (c) Alberto Ruiz 2011 +-- 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.Vectorized +import Data.Packed.Vector +import Data.Packed.Numeric + +------------------------------------------------------------------- + +adaptScalar f1 f2 f3 x y + | dim x == 1 = f1 (x@>0) y + | dim y == 1 = f3 x (y@>0) + | otherwise = f2 x y + +------------------------------------------------------------------ + +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] + diff --git a/packages/hmatrix/Config.hs b/packages/hmatrix/Config.hs deleted file mode 100644 index 5145d1a..0000000 --- a/packages/hmatrix/Config.hs +++ /dev/null @@ -1,154 +0,0 @@ -{- - GSL and LAPACK may require auxiliary libraries which depend on OS, - distribution, and implementation. This script tries to to find out - the correct link command for your system. - Suggestions and contributions are welcome. - - By default we try to link -lgsl -llapack. This works in ubuntu/debian, - both with and without ATLAS. - If this fails we try different sets of additional libraries which are - known to work in some systems. - - The desired libraries can also be explicitly given by the user using cabal - flags (e.g., -fmkl, -faccelerate) or --configure-option=link:lib1,lib2,lib3,... - --} - -module Config(config) where - -import System.Process -import System.Exit -import System.Environment -import System.Directory(createDirectoryIfMissing) -import System.FilePath(()) -import Data.List(isPrefixOf, intercalate) -import Distribution.Simple.LocalBuildInfo -import Distribution.Simple.Configure -import Distribution.PackageDescription - --- possible additional dependencies for the desired libs (by default gsl lapack) - -opts = [ "" -- Ubuntu/Debian - , "blas" - , "blas cblas" - , "cblas" - , "gslcblas" - , "blas gslcblas" - , "f77blas" - , "f77blas cblas atlas gcc_s" -- Arch Linux (older version of atlas-lapack) - , "blas gslcblas gfortran" -- Arch Linux with normal blas and lapack - ] - --- location of test program -testProgLoc bInfo = buildDir bInfo "dummy.c" -testOutLoc bInfo = buildDir bInfo "dummy" - --- write test program -writeTestProg bInfo contents = writeFile (testProgLoc bInfo) contents - --- compile, discarding error messages -compile cmd = do - let processRecord = (shell $ unwords cmd) { std_out = CreatePipe - , std_err = CreatePipe } - ( _, _, _, h) <- createProcess processRecord - waitForProcess h - --- command to compile the test program -compileCmd bInfo buildInfo = [ "gcc " - , (unwords $ ccOptions buildInfo) - , (unwords $ cppOptions buildInfo) - , (unwords $ map ("-I"++) $ includeDirs buildInfo) - , testProgLoc bInfo - , "-o" - , testOutLoc bInfo - , (unwords $ map ("-L"++) $ extraLibDirs buildInfo) - ] - --- compile a simple program with symbols from GSL and LAPACK with the given libs -testprog bInfo buildInfo libs fmks = do - writeTestProg bInfo "#include \nint main(){dgemm_(); zgesvd_(); gsl_sf_gamma(5);}" - compile $ compileCmd bInfo - buildInfo - ++ [ (prepend "-l" $ libs) - , (prepend "-framework " fmks) ] - -prepend x = unwords . map (x++) . words - -check bInfo buildInfo libs fmks = (ExitSuccess ==) `fmap` testprog bInfo buildInfo libs fmks - --- simple test for GSL -gsl bInfo buildInfo = do - writeTestProg bInfo "#include \nint main(){gsl_sf_gamma(5);}" - compile $ compileCmd bInfo buildInfo ++ ["-lgsl", "-lgslcblas"] - --- test for gsl >= 1.12 -gsl112 bInfo buildInfo = do - writeTestProg bInfo "#include \nint main(){gsl_sf_exprel_n_CF_e(1,1,0);}" - compile $ compileCmd bInfo buildInfo ++ ["-lgsl", "-lgslcblas"] - --- test for odeiv2 -gslodeiv2 bInfo buildInfo = do - writeTestProg bInfo "#include \nint main(){return 0;}" - compile $ compileCmd bInfo buildInfo ++ ["-lgsl", "-lgslcblas"] - -checkCommand c = (ExitSuccess ==) `fmap` c - --- test different configurations until the first one works -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 l i b f rest - --- read --configure-option=link:lib1,lib2,lib3,etc -linkop = "--configure-option=link:" -getUserLink = concatMap (g . drop (length linkop)) . filter (isPrefixOf linkop) - where g = map cs - cs ',' = ' ' - cs x = x - -config :: LocalBuildInfo -> IO HookedBuildInfo - -config bInfo = do - putStr "Checking foreign libraries..." - args <- getArgs - - let Just lib = library . localPkgDescr $ bInfo - buildInfo = libBuildInfo lib - base = unwords . extraLibs $ buildInfo - fwks = unwords . frameworks $ buildInfo - auxpref = getUserLink args - - -- We extract the desired libs from hmatrix.cabal (using a cabal flags) - -- and from a posible --configure-option=link:lib1,lib2,lib3 - -- by default the desired libs are gsl lapack. - - let pref = if null (words (base ++ " " ++ auxpref)) then "gsl lapack" else auxpref - fullOpts = map ((pref++" ")++) opts - - -- 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 bInfo buildInfo - if g - then putStrLn " *** Sorry, I can't link LAPACK." - else putStrLn " *** Sorry, I can't link GSL." - putStrLn " *** Please make sure that the appropriate -dev packages are installed." - putStrLn " *** You can also specify the required libraries using" - putStrLn " *** cabal install hmatrix --configure-option=link:lib1,lib2,lib3,etc." - return (Just emptyBuildInfo { buildable = False }, []) - Just ops -> do - putStrLn $ " OK " ++ ops - g1 <- checkCommand $ gsl112 bInfo buildInfo - let op1 = if g1 then "" else "-DGSL110" - g2 <- checkCommand $ gslodeiv2 bInfo buildInfo - let op2 = if g2 then "" else "-DGSLODE1" - opts = filter (not.null) [op1,op2] - let hbi = emptyBuildInfo { extraLibs = words ops, ccOptions = opts } - return (Just hbi, []) - diff --git a/packages/hmatrix/Setup.lhs b/packages/hmatrix/Setup.lhs index e3f9847..4b19c19 100644 --- a/packages/hmatrix/Setup.lhs +++ b/packages/hmatrix/Setup.lhs @@ -1,18 +1,5 @@ #! /usr/bin/env runhaskell > import Distribution.Simple -> import Distribution.Simple.Setup -> import Distribution.PackageDescription -> import Distribution.Simple.LocalBuildInfo - -> import System.Process(system) -> import Config(config) - -> main = defaultMainWithHooks simpleUserHooks { confHook = c } - -> c x y = do -> binfo <- confHook simpleUserHooks x y -> pbi <- config binfo -> let pkg_descr = localPkgDescr binfo -> return $ binfo { localPkgDescr = updatePackageDescription pbi pkg_descr } +> main = defaultMain diff --git a/packages/hmatrix/hmatrix.cabal b/packages/hmatrix/hmatrix.cabal index 73edac4..369c412 100644 --- a/packages/hmatrix/hmatrix.cabal +++ b/packages/hmatrix/hmatrix.cabal @@ -102,9 +102,6 @@ library Numeric.GSL.Internal, Numeric.GSL.Vector, Numeric.IO, - Numeric.Chain, - Numeric.Vector, - Numeric.Matrix, Numeric.LinearAlgebra.Util.Convolution diff --git a/packages/hmatrix/src/Numeric/Chain.hs b/packages/hmatrix/src/Numeric/Chain.hs deleted file mode 100644 index de6a86f..0000000 --- a/packages/hmatrix/src/Numeric/Chain.hs +++ /dev/null @@ -1,140 +0,0 @@ ------------------------------------------------------------------------------ --- | --- 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 Data.Packed.Numeric - -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/packages/hmatrix/src/Numeric/Container.hs b/packages/hmatrix/src/Numeric/Container.hs index 645a83f..e7f23d4 100644 --- a/packages/hmatrix/src/Numeric/Container.hs +++ b/packages/hmatrix/src/Numeric/Container.hs @@ -66,11 +66,11 @@ module Numeric.Container ( import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ) import Data.Packed.Numeric -import Numeric.Chain import Numeric.IO import Data.Complex import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD) import Numeric.Random +import Data.Monoid(Monoid(mconcat)) ------------------------------------------------------------------ @@ -268,4 +268,8 @@ infixl 7 ◇ dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t dot u v = udot (conj u) v +-------------------------------------------------------------------------------- + +optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t +optimiseMult = mconcat diff --git a/packages/hmatrix/src/Numeric/Matrix.hs b/packages/hmatrix/src/Numeric/Matrix.hs deleted file mode 100644 index e285ff2..0000000 --- a/packages/hmatrix/src/Numeric/Matrix.hs +++ /dev/null @@ -1,98 +0,0 @@ -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} - ------------------------------------------------------------------------------ --- | --- 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 -import qualified Data.Monoid as M -import Data.List(partition) - -------------------------------------------------------------------- - -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] - --------------------------------------------------------------------------------- - -isScalar m = rows m == 1 && cols m == 1 - -adaptScalarM f1 f2 f3 x y - | isScalar x = f1 (x @@>(0,0) ) y - | isScalar y = f3 x (y @@>(0,0) ) - | otherwise = f2 x y - -instance (Container Vector t, Eq t, Num (Vector t), Product t) => M.Monoid (Matrix t) - where - mempty = 1 - mappend = adaptScalarM scale mXm (flip scale) - - mconcat xs = work (partition isScalar xs) - where - work (ss,[]) = product ss - work (ss,ms) = scale' (product ss) (optimiseMult ms) - scale' x m - | isScalar x && x00 == 1 = m - | otherwise = scale x00 m - where - x00 = x @@> (0,0) - diff --git a/packages/hmatrix/src/Numeric/Vector.hs b/packages/hmatrix/src/Numeric/Vector.hs deleted file mode 100644 index 4c59d32..0000000 --- a/packages/hmatrix/src/Numeric/Vector.hs +++ /dev/null @@ -1,158 +0,0 @@ -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} ------------------------------------------------------------------------------ --- | --- Module : Numeric.Vector --- Copyright : (c) Alberto Ruiz 2011 --- 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.Vectorized -import Numeric.Container - -------------------------------------------------------------------- - -adaptScalar f1 f2 f3 x y - | dim x == 1 = f1 (x@>0) y - | dim y == 1 = f3 x (y@>0) - | otherwise = f2 x y - ------------------------------------------------------------------- - -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] - -- cgit v1.2.3