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/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 ------------------------------ 4 files changed, 5 insertions(+), 397 deletions(-) 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 (limited to 'packages/hmatrix/src') 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