diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-16 09:20:47 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-16 09:20:47 +0200 |
commit | 0ab93b8eb934167e16dbae193c4fd2b5359a184b (patch) | |
tree | e5d7236890c8310dae90952bfb5f6195dd97295c /packages/base | |
parent | a86c60a5fbfc73ff3080c88007625c2cd094e80f (diff) |
instances moved to base
Diffstat (limited to 'packages/base')
-rw-r--r-- | packages/base/hmatrix-base.cabal | 3 | ||||
-rw-r--r-- | packages/base/src/Numeric/Chain.hs | 143 | ||||
-rw-r--r-- | packages/base/src/Numeric/Matrix.hs | 100 | ||||
-rw-r--r-- | packages/base/src/Numeric/Vector.hs | 159 |
4 files changed, 405 insertions, 0 deletions
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 | |||
39 | Numeric.LinearAlgebra.LAPACK | 39 | Numeric.LinearAlgebra.LAPACK |
40 | Data.Packed.Numeric | 40 | Data.Packed.Numeric |
41 | Numeric.Vectorized | 41 | Numeric.Vectorized |
42 | Numeric.Vector | ||
43 | Numeric.Matrix | ||
42 | 44 | ||
43 | other-modules: Data.Packed.Internal, | 45 | other-modules: Data.Packed.Internal, |
44 | Data.Packed.Internal.Common, | 46 | Data.Packed.Internal.Common, |
45 | Data.Packed.Internal.Signatures, | 47 | Data.Packed.Internal.Signatures, |
46 | Data.Packed.Internal.Vector, | 48 | Data.Packed.Internal.Vector, |
47 | Data.Packed.Internal.Matrix | 49 | Data.Packed.Internal.Matrix |
50 | Numeric.Chain | ||
48 | 51 | ||
49 | C-sources: src/C/lapack-aux.c | 52 | C-sources: src/C/lapack-aux.c |
50 | src/C/vector-aux.c | 53 | 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 @@ | |||
1 | ----------------------------------------------------------------------------- | ||
2 | -- | | ||
3 | -- Module : Numeric.Chain | ||
4 | -- Copyright : (c) Vivian McPhail 2010 | ||
5 | -- License : GPL-style | ||
6 | -- | ||
7 | -- Maintainer : Vivian McPhail <haskell.vivian.mcphail <at> gmail.com> | ||
8 | -- Stability : provisional | ||
9 | -- Portability : portable | ||
10 | -- | ||
11 | -- optimisation of association order for chains of matrix multiplication | ||
12 | -- | ||
13 | ----------------------------------------------------------------------------- | ||
14 | |||
15 | module Numeric.Chain ( | ||
16 | optimiseMult, | ||
17 | ) where | ||
18 | |||
19 | import Data.Maybe | ||
20 | |||
21 | import Data.Packed.Matrix | ||
22 | import Data.Packed.Numeric | ||
23 | |||
24 | import qualified Data.Array.IArray as A | ||
25 | |||
26 | ----------------------------------------------------------------------------- | ||
27 | {- | | ||
28 | Provide optimal association order for a chain of matrix multiplications | ||
29 | and apply the multiplications. | ||
30 | |||
31 | The algorithm is the well-known O(n\^3) dynamic programming algorithm | ||
32 | that builds a pyramid of optimal associations. | ||
33 | |||
34 | > m1, m2, m3, m4 :: Matrix Double | ||
35 | > m1 = (10><15) [1..] | ||
36 | > m2 = (15><20) [1..] | ||
37 | > m3 = (20><5) [1..] | ||
38 | > m4 = (5><10) [1..] | ||
39 | |||
40 | > >>> optimiseMult [m1,m2,m3,m4] | ||
41 | |||
42 | will perform @((m1 `multiply` (m2 `multiply` m3)) `multiply` m4)@ | ||
43 | |||
44 | The naive left-to-right multiplication would take @4500@ scalar multiplications | ||
45 | whereas the optimised version performs @2750@ scalar multiplications. The complexity | ||
46 | in this case is 32 (= 4^3/2) * (2 comparisons, 3 scalar multiplications, 3 scalar additions, | ||
47 | 5 lookups, 2 updates) + a constant (= three table allocations) | ||
48 | -} | ||
49 | optimiseMult :: Product t => [Matrix t] -> Matrix t | ||
50 | optimiseMult = chain | ||
51 | |||
52 | ----------------------------------------------------------------------------- | ||
53 | |||
54 | type Matrices a = A.Array Int (Matrix a) | ||
55 | type Sizes = A.Array Int (Int,Int) | ||
56 | type Cost = A.Array Int (A.Array Int (Maybe Int)) | ||
57 | type Indexes = A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int)))) | ||
58 | |||
59 | update :: A.Array Int (A.Array Int a) -> (Int,Int) -> a -> A.Array Int (A.Array Int a) | ||
60 | update a (r,c) e = a A.// [(r,(a A.! r) A.// [(c,e)])] | ||
61 | |||
62 | newWorkSpaceCost :: Int -> A.Array Int (A.Array Int (Maybe Int)) | ||
63 | newWorkSpaceCost n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n] | ||
64 | where subArray i = A.listArray (1,i) (repeat Nothing) | ||
65 | |||
66 | newWorkSpaceIndexes :: Int -> A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int)))) | ||
67 | newWorkSpaceIndexes n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n] | ||
68 | where subArray i = A.listArray (1,i) (repeat Nothing) | ||
69 | |||
70 | matricesToSizes :: [Matrix a] -> Sizes | ||
71 | matricesToSizes ms = A.listArray (1,length ms) $ map (\m -> (rows m,cols m)) ms | ||
72 | |||
73 | chain :: Product a => [Matrix a] -> Matrix a | ||
74 | chain [] = error "chain: zero matrices to multiply" | ||
75 | chain [m] = m | ||
76 | chain [ml,mr] = ml `multiply` mr | ||
77 | chain ms = let ln = length ms | ||
78 | ma = A.listArray (1,ln) ms | ||
79 | mz = matricesToSizes ms | ||
80 | i = chain_cost mz | ||
81 | in chain_paren (ln,ln) i ma | ||
82 | |||
83 | chain_cost :: Sizes -> Indexes | ||
84 | chain_cost mz = let (_,u) = A.bounds mz | ||
85 | cost = newWorkSpaceCost u | ||
86 | ixes = newWorkSpaceIndexes u | ||
87 | (_,_,i) = foldl chain_cost' (mz,cost,ixes) (order u) | ||
88 | in i | ||
89 | |||
90 | chain_cost' :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes) | ||
91 | chain_cost' sci@(mz,cost,ixes) (r,c) | ||
92 | | c == 1 = let cost' = update cost (r,c) (Just 0) | ||
93 | ixes' = update ixes (r,c) (Just ((r,c),(r,c))) | ||
94 | in (mz,cost',ixes') | ||
95 | | otherwise = minimum_cost sci (r,c) | ||
96 | |||
97 | minimum_cost :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes) | ||
98 | minimum_cost sci fu = foldl (smaller_cost fu) sci (fulcrum_order fu) | ||
99 | |||
100 | smaller_cost :: (Int,Int) -> (Sizes,Cost,Indexes) -> ((Int,Int),(Int,Int)) -> (Sizes,Cost,Indexes) | ||
101 | smaller_cost (r,c) (mz,cost,ixes) ix@((lr,lc),(rr,rc)) = let op_cost = fromJust ((cost A.! lr) A.! lc) | ||
102 | + fromJust ((cost A.! rr) A.! rc) | ||
103 | + fst (mz A.! (lr-lc+1)) | ||
104 | * snd (mz A.! lc) | ||
105 | * snd (mz A.! rr) | ||
106 | cost' = (cost A.! r) A.! c | ||
107 | in case cost' of | ||
108 | Nothing -> let cost'' = update cost (r,c) (Just op_cost) | ||
109 | ixes'' = update ixes (r,c) (Just ix) | ||
110 | in (mz,cost'',ixes'') | ||
111 | Just ct -> if op_cost < ct then | ||
112 | let cost'' = update cost (r,c) (Just op_cost) | ||
113 | ixes'' = update ixes (r,c) (Just ix) | ||
114 | in (mz,cost'',ixes'') | ||
115 | else (mz,cost,ixes) | ||
116 | |||
117 | |||
118 | fulcrum_order (r,c) = let fs' = zip (repeat r) [1..(c-1)] | ||
119 | in map (partner (r,c)) fs' | ||
120 | |||
121 | partner (r,c) (a,b) = ((r-b, c-b), (a,b)) | ||
122 | |||
123 | order 0 = [] | ||
124 | order n = order (n-1) ++ zip (repeat n) [1..n] | ||
125 | |||
126 | chain_paren :: Product a => (Int,Int) -> Indexes -> Matrices a -> Matrix a | ||
127 | chain_paren (r,c) ixes ma = let ((lr,lc),(rr,rc)) = fromJust $ (ixes A.! r) A.! c | ||
128 | in if lr == rr && lc == rc then (ma A.! lr) | ||
129 | else (chain_paren (lr,lc) ixes ma) `multiply` (chain_paren (rr,rc) ixes ma) | ||
130 | |||
131 | -------------------------------------------------------------------------- | ||
132 | |||
133 | {- TESTS | ||
134 | |||
135 | -- optimal association is ((m1*(m2*m3))*m4) | ||
136 | m1, m2, m3, m4 :: Matrix Double | ||
137 | m1 = (10><15) [1..] | ||
138 | m2 = (15><20) [1..] | ||
139 | m3 = (20><5) [1..] | ||
140 | m4 = (5><10) [1..] | ||
141 | |||
142 | -} | ||
143 | |||
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 @@ | |||
1 | {-# LANGUAGE TypeFamilies #-} | ||
2 | {-# LANGUAGE FlexibleContexts #-} | ||
3 | {-# LANGUAGE FlexibleInstances #-} | ||
4 | {-# LANGUAGE UndecidableInstances #-} | ||
5 | {-# LANGUAGE MultiParamTypeClasses #-} | ||
6 | |||
7 | ----------------------------------------------------------------------------- | ||
8 | -- | | ||
9 | -- Module : Numeric.Matrix | ||
10 | -- Copyright : (c) Alberto Ruiz 2010 | ||
11 | -- License : GPL-style | ||
12 | -- | ||
13 | -- Maintainer : Alberto Ruiz <aruiz@um.es> | ||
14 | -- Stability : provisional | ||
15 | -- Portability : portable | ||
16 | -- | ||
17 | -- Provides instances of standard classes 'Show', 'Read', 'Eq', | ||
18 | -- 'Num', 'Fractional', and 'Floating' for 'Matrix'. | ||
19 | -- | ||
20 | -- In arithmetic operations one-component | ||
21 | -- vectors and matrices automatically expand to match the dimensions of the other operand. | ||
22 | |||
23 | ----------------------------------------------------------------------------- | ||
24 | |||
25 | module Numeric.Matrix ( | ||
26 | ) where | ||
27 | |||
28 | ------------------------------------------------------------------- | ||
29 | |||
30 | import Data.Packed | ||
31 | import Data.Packed.Numeric | ||
32 | import qualified Data.Monoid as M | ||
33 | import Data.List(partition) | ||
34 | import Numeric.Chain | ||
35 | |||
36 | ------------------------------------------------------------------- | ||
37 | |||
38 | instance Container Matrix a => Eq (Matrix a) where | ||
39 | (==) = equal | ||
40 | |||
41 | instance (Container Matrix a, Num (Vector a)) => Num (Matrix a) where | ||
42 | (+) = liftMatrix2Auto (+) | ||
43 | (-) = liftMatrix2Auto (-) | ||
44 | negate = liftMatrix negate | ||
45 | (*) = liftMatrix2Auto (*) | ||
46 | signum = liftMatrix signum | ||
47 | abs = liftMatrix abs | ||
48 | fromInteger = (1><1) . return . fromInteger | ||
49 | |||
50 | --------------------------------------------------- | ||
51 | |||
52 | instance (Container Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where | ||
53 | fromRational n = (1><1) [fromRational n] | ||
54 | (/) = liftMatrix2Auto (/) | ||
55 | |||
56 | --------------------------------------------------------- | ||
57 | |||
58 | instance (Floating a, Container Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where | ||
59 | sin = liftMatrix sin | ||
60 | cos = liftMatrix cos | ||
61 | tan = liftMatrix tan | ||
62 | asin = liftMatrix asin | ||
63 | acos = liftMatrix acos | ||
64 | atan = liftMatrix atan | ||
65 | sinh = liftMatrix sinh | ||
66 | cosh = liftMatrix cosh | ||
67 | tanh = liftMatrix tanh | ||
68 | asinh = liftMatrix asinh | ||
69 | acosh = liftMatrix acosh | ||
70 | atanh = liftMatrix atanh | ||
71 | exp = liftMatrix exp | ||
72 | log = liftMatrix log | ||
73 | (**) = liftMatrix2Auto (**) | ||
74 | sqrt = liftMatrix sqrt | ||
75 | pi = (1><1) [pi] | ||
76 | |||
77 | -------------------------------------------------------------------------------- | ||
78 | |||
79 | isScalar m = rows m == 1 && cols m == 1 | ||
80 | |||
81 | adaptScalarM f1 f2 f3 x y | ||
82 | | isScalar x = f1 (x @@>(0,0) ) y | ||
83 | | isScalar y = f3 x (y @@>(0,0) ) | ||
84 | | otherwise = f2 x y | ||
85 | |||
86 | instance (Container Vector t, Eq t, Num (Vector t), Product t) => M.Monoid (Matrix t) | ||
87 | where | ||
88 | mempty = 1 | ||
89 | mappend = adaptScalarM scale mXm (flip scale) | ||
90 | |||
91 | mconcat xs = work (partition isScalar xs) | ||
92 | where | ||
93 | work (ss,[]) = product ss | ||
94 | work (ss,ms) = scale' (product ss) (optimiseMult ms) | ||
95 | scale' x m | ||
96 | | isScalar x && x00 == 1 = m | ||
97 | | otherwise = scale x00 m | ||
98 | where | ||
99 | x00 = x @@> (0,0) | ||
100 | |||
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 @@ | |||
1 | {-# LANGUAGE TypeFamilies #-} | ||
2 | {-# LANGUAGE FlexibleContexts #-} | ||
3 | {-# LANGUAGE FlexibleInstances #-} | ||
4 | {-# LANGUAGE UndecidableInstances #-} | ||
5 | {-# LANGUAGE MultiParamTypeClasses #-} | ||
6 | ----------------------------------------------------------------------------- | ||
7 | -- | | ||
8 | -- Module : Numeric.Vector | ||
9 | -- Copyright : (c) Alberto Ruiz 2011 | ||
10 | -- License : GPL-style | ||
11 | -- | ||
12 | -- Maintainer : Alberto Ruiz <aruiz@um.es> | ||
13 | -- Stability : provisional | ||
14 | -- Portability : portable | ||
15 | -- | ||
16 | -- Provides instances of standard classes 'Show', 'Read', 'Eq', | ||
17 | -- 'Num', 'Fractional', and 'Floating' for 'Vector'. | ||
18 | -- | ||
19 | ----------------------------------------------------------------------------- | ||
20 | |||
21 | module Numeric.Vector () where | ||
22 | |||
23 | import Numeric.Vectorized | ||
24 | import Data.Packed.Vector | ||
25 | import Data.Packed.Numeric | ||
26 | |||
27 | ------------------------------------------------------------------- | ||
28 | |||
29 | adaptScalar f1 f2 f3 x y | ||
30 | | dim x == 1 = f1 (x@>0) y | ||
31 | | dim y == 1 = f3 x (y@>0) | ||
32 | | otherwise = f2 x y | ||
33 | |||
34 | ------------------------------------------------------------------ | ||
35 | |||
36 | instance Num (Vector Float) where | ||
37 | (+) = adaptScalar addConstant add (flip addConstant) | ||
38 | negate = scale (-1) | ||
39 | (*) = adaptScalar scale mul (flip scale) | ||
40 | signum = vectorMapF Sign | ||
41 | abs = vectorMapF Abs | ||
42 | fromInteger = fromList . return . fromInteger | ||
43 | |||
44 | instance Num (Vector Double) where | ||
45 | (+) = adaptScalar addConstant add (flip addConstant) | ||
46 | negate = scale (-1) | ||
47 | (*) = adaptScalar scale mul (flip scale) | ||
48 | signum = vectorMapR Sign | ||
49 | abs = vectorMapR Abs | ||
50 | fromInteger = fromList . return . fromInteger | ||
51 | |||
52 | instance Num (Vector (Complex Double)) where | ||
53 | (+) = adaptScalar addConstant add (flip addConstant) | ||
54 | negate = scale (-1) | ||
55 | (*) = adaptScalar scale mul (flip scale) | ||
56 | signum = vectorMapC Sign | ||
57 | abs = vectorMapC Abs | ||
58 | fromInteger = fromList . return . fromInteger | ||
59 | |||
60 | instance Num (Vector (Complex Float)) where | ||
61 | (+) = adaptScalar addConstant add (flip addConstant) | ||
62 | negate = scale (-1) | ||
63 | (*) = adaptScalar scale mul (flip scale) | ||
64 | signum = vectorMapQ Sign | ||
65 | abs = vectorMapQ Abs | ||
66 | fromInteger = fromList . return . fromInteger | ||
67 | |||
68 | --------------------------------------------------- | ||
69 | |||
70 | instance (Container Vector a, Num (Vector a)) => Fractional (Vector a) where | ||
71 | fromRational n = fromList [fromRational n] | ||
72 | (/) = adaptScalar f divide g where | ||
73 | r `f` v = scaleRecip r v | ||
74 | v `g` r = scale (recip r) v | ||
75 | |||
76 | ------------------------------------------------------- | ||
77 | |||
78 | instance Floating (Vector Float) where | ||
79 | sin = vectorMapF Sin | ||
80 | cos = vectorMapF Cos | ||
81 | tan = vectorMapF Tan | ||
82 | asin = vectorMapF ASin | ||
83 | acos = vectorMapF ACos | ||
84 | atan = vectorMapF ATan | ||
85 | sinh = vectorMapF Sinh | ||
86 | cosh = vectorMapF Cosh | ||
87 | tanh = vectorMapF Tanh | ||
88 | asinh = vectorMapF ASinh | ||
89 | acosh = vectorMapF ACosh | ||
90 | atanh = vectorMapF ATanh | ||
91 | exp = vectorMapF Exp | ||
92 | log = vectorMapF Log | ||
93 | sqrt = vectorMapF Sqrt | ||
94 | (**) = adaptScalar (vectorMapValF PowSV) (vectorZipF Pow) (flip (vectorMapValF PowVS)) | ||
95 | pi = fromList [pi] | ||
96 | |||
97 | ------------------------------------------------------------- | ||
98 | |||
99 | instance Floating (Vector Double) where | ||
100 | sin = vectorMapR Sin | ||
101 | cos = vectorMapR Cos | ||
102 | tan = vectorMapR Tan | ||
103 | asin = vectorMapR ASin | ||
104 | acos = vectorMapR ACos | ||
105 | atan = vectorMapR ATan | ||
106 | sinh = vectorMapR Sinh | ||
107 | cosh = vectorMapR Cosh | ||
108 | tanh = vectorMapR Tanh | ||
109 | asinh = vectorMapR ASinh | ||
110 | acosh = vectorMapR ACosh | ||
111 | atanh = vectorMapR ATanh | ||
112 | exp = vectorMapR Exp | ||
113 | log = vectorMapR Log | ||
114 | sqrt = vectorMapR Sqrt | ||
115 | (**) = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS)) | ||
116 | pi = fromList [pi] | ||
117 | |||
118 | ------------------------------------------------------------- | ||
119 | |||
120 | instance Floating (Vector (Complex Double)) where | ||
121 | sin = vectorMapC Sin | ||
122 | cos = vectorMapC Cos | ||
123 | tan = vectorMapC Tan | ||
124 | asin = vectorMapC ASin | ||
125 | acos = vectorMapC ACos | ||
126 | atan = vectorMapC ATan | ||
127 | sinh = vectorMapC Sinh | ||
128 | cosh = vectorMapC Cosh | ||
129 | tanh = vectorMapC Tanh | ||
130 | asinh = vectorMapC ASinh | ||
131 | acosh = vectorMapC ACosh | ||
132 | atanh = vectorMapC ATanh | ||
133 | exp = vectorMapC Exp | ||
134 | log = vectorMapC Log | ||
135 | sqrt = vectorMapC Sqrt | ||
136 | (**) = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS)) | ||
137 | pi = fromList [pi] | ||
138 | |||
139 | ----------------------------------------------------------- | ||
140 | |||
141 | instance Floating (Vector (Complex Float)) where | ||
142 | sin = vectorMapQ Sin | ||
143 | cos = vectorMapQ Cos | ||
144 | tan = vectorMapQ Tan | ||
145 | asin = vectorMapQ ASin | ||
146 | acos = vectorMapQ ACos | ||
147 | atan = vectorMapQ ATan | ||
148 | sinh = vectorMapQ Sinh | ||
149 | cosh = vectorMapQ Cosh | ||
150 | tanh = vectorMapQ Tanh | ||
151 | asinh = vectorMapQ ASinh | ||
152 | acosh = vectorMapQ ACosh | ||
153 | atanh = vectorMapQ ATanh | ||
154 | exp = vectorMapQ Exp | ||
155 | log = vectorMapQ Log | ||
156 | sqrt = vectorMapQ Sqrt | ||
157 | (**) = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS)) | ||
158 | pi = fromList [pi] | ||
159 | |||