diff options
Diffstat (limited to 'lib/Numeric/Container.hs')
-rw-r--r-- | lib/Numeric/Container.hs | 418 |
1 files changed, 254 insertions, 164 deletions
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs index d6b1342..fc4fb0d 100644 --- a/lib/Numeric/Container.hs +++ b/lib/Numeric/Container.hs | |||
@@ -3,6 +3,7 @@ | |||
3 | {-# LANGUAGE FlexibleInstances #-} | 3 | {-# LANGUAGE FlexibleInstances #-} |
4 | {-# LANGUAGE MultiParamTypeClasses #-} | 4 | {-# LANGUAGE MultiParamTypeClasses #-} |
5 | {-# LANGUAGE FunctionalDependencies #-} | 5 | {-# LANGUAGE FunctionalDependencies #-} |
6 | {-# LANGUAGE UndecidableInstances #-} | ||
6 | 7 | ||
7 | ----------------------------------------------------------------------------- | 8 | ----------------------------------------------------------------------------- |
8 | -- | | 9 | -- | |
@@ -19,103 +20,35 @@ | |||
19 | ----------------------------------------------------------------------------- | 20 | ----------------------------------------------------------------------------- |
20 | 21 | ||
21 | module Numeric.Container ( | 22 | module Numeric.Container ( |
22 | Linear(..), | 23 | --Linear(..), |
23 | Container(..), RealElement, Precision(..), NumericContainer(..), comp, | 24 | Container(..), |
24 | Convert(..), --AutoReal(..), | 25 | Vectors(..), |
25 | RealOf, ComplexOf, SingleOf, DoubleOf, | 26 | Product(..), |
26 | 27 | mXm,mXv,vXm, | |
27 | -- ElementOf, | 28 | outer, kronecker, |
29 | |||
30 | RealElement, --Precision, | ||
31 | ComplexContainer(toComplex,fromComplex,comp,conj), | ||
32 | Convert(..), --AutoReal(..), | ||
33 | RealOf, ComplexOf, SingleOf, DoubleOf, | ||
28 | 34 | ||
29 | IndexOf, | 35 | IndexOf, |
30 | |||
31 | module Data.Complex | 36 | module Data.Complex |
32 | ) where | 37 | ) where |
33 | 38 | ||
34 | import Data.Packed.Vector | 39 | import Data.Packed |
35 | import Data.Packed.Matrix | 40 | import Numeric.Conversion |
36 | import Data.Packed.Internal.Vector | 41 | import Data.Packed.Internal |
37 | --import Data.Packed.Internal.Matrix | 42 | import Numeric.GSL.Vector |
38 | --import qualified Data.Packed.ST as ST | ||
39 | 43 | ||
40 | --import Control.Arrow((***)) | ||
41 | import Data.Complex | 44 | import Data.Complex |
45 | import Control.Monad(ap) | ||
46 | --import Control.Arrow((***)) | ||
42 | 47 | ||
43 | ------------------------------------------------------------------- | 48 | import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) |
44 | |||
45 | -- | Supported single-double precision type pairs | ||
46 | class (Element s, Element d) => Precision s d | s -> d, d -> s where | ||
47 | double2FloatG :: Vector d -> Vector s | ||
48 | float2DoubleG :: Vector s -> Vector d | ||
49 | |||
50 | instance Precision Float Double where | ||
51 | double2FloatG = double2FloatV | ||
52 | float2DoubleG = float2DoubleV | ||
53 | |||
54 | instance Precision (Complex Float) (Complex Double) where | ||
55 | double2FloatG = asComplex . double2FloatV . asReal | ||
56 | float2DoubleG = asComplex . float2DoubleV . asReal | ||
57 | |||
58 | -- | Supported real types | ||
59 | class (Element t, Element (Complex t), RealFloat t | ||
60 | -- , RealOf t ~ t, RealOf (Complex t) ~ t | ||
61 | ) | ||
62 | => RealElement t | ||
63 | |||
64 | instance RealElement Double | ||
65 | |||
66 | instance RealElement Float | ||
67 | |||
68 | -- | Conversion utilities | ||
69 | class NumericContainer c where | ||
70 | toComplex :: (RealElement e) => (c e, c e) -> c (Complex e) | ||
71 | fromComplex :: (RealElement e) => c (Complex e) -> (c e, c e) | ||
72 | complex' :: (RealElement e) => c e -> c (Complex e) | ||
73 | conj :: (RealElement e) => c (Complex e) -> c (Complex e) | ||
74 | -- cmap :: (Element a, Element b) => (a -> b) -> c a -> c b | ||
75 | single' :: Precision a b => c b -> c a | ||
76 | double' :: Precision a b => c a -> c b | ||
77 | |||
78 | -- | a synonym for "complex'" | ||
79 | comp :: (NumericContainer c, RealElement e) => c e -> c (Complex e) | ||
80 | comp x = complex' x | ||
81 | 49 | ||
82 | ------------------------------------------------------------------- | 50 | ------------------------------------------------------------------- |
83 | 51 | ||
84 | type family RealOf x | ||
85 | |||
86 | type instance RealOf Double = Double | ||
87 | type instance RealOf (Complex Double) = Double | ||
88 | |||
89 | type instance RealOf Float = Float | ||
90 | type instance RealOf (Complex Float) = Float | ||
91 | |||
92 | type family ComplexOf x | ||
93 | |||
94 | type instance ComplexOf Double = Complex Double | ||
95 | type instance ComplexOf (Complex Double) = Complex Double | ||
96 | |||
97 | type instance ComplexOf Float = Complex Float | ||
98 | type instance ComplexOf (Complex Float) = Complex Float | ||
99 | |||
100 | type family SingleOf x | ||
101 | |||
102 | type instance SingleOf Double = Float | ||
103 | type instance SingleOf Float = Float | ||
104 | |||
105 | type instance SingleOf (Complex a) = Complex (SingleOf a) | ||
106 | |||
107 | type family DoubleOf x | ||
108 | |||
109 | type instance DoubleOf Double = Double | ||
110 | type instance DoubleOf Float = Double | ||
111 | |||
112 | type instance DoubleOf (Complex a) = Complex (DoubleOf a) | ||
113 | |||
114 | type family ElementOf c | ||
115 | |||
116 | type instance ElementOf (Vector a) = a | ||
117 | type instance ElementOf (Matrix a) = a | ||
118 | |||
119 | type family IndexOf c | 52 | type family IndexOf c |
120 | 53 | ||
121 | type instance IndexOf Vector = Int | 54 | type instance IndexOf Vector = Int |
@@ -123,75 +56,16 @@ type instance IndexOf Matrix = (Int,Int) | |||
123 | 56 | ||
124 | ------------------------------------------------------------------- | 57 | ------------------------------------------------------------------- |
125 | 58 | ||
126 | class (Element t, Element (RealOf t)) => Convert t where | ||
127 | real :: NumericContainer c => c (RealOf t) -> c t | ||
128 | complex :: NumericContainer c => c t -> c (ComplexOf t) | ||
129 | single :: NumericContainer c => c t -> c (SingleOf t) | ||
130 | double :: NumericContainer c => c t -> c (DoubleOf t) | ||
131 | |||
132 | |||
133 | instance Convert Double where | ||
134 | real = id | ||
135 | complex = complex' | ||
136 | single = single' | ||
137 | double = id | ||
138 | |||
139 | instance Convert Float where | ||
140 | real = id | ||
141 | complex = complex' | ||
142 | single = id | ||
143 | double = double' | ||
144 | |||
145 | instance Convert (Complex Double) where | ||
146 | real = complex' | ||
147 | complex = id | ||
148 | single = single' | ||
149 | double = id | ||
150 | |||
151 | instance Convert (Complex Float) where | ||
152 | real = complex' | ||
153 | complex = id | ||
154 | single = id | ||
155 | double = double' | ||
156 | |||
157 | ------------------------------------------------------------------- | ||
158 | |||
159 | -- | to be replaced by Convert | ||
160 | class Convert t => AutoReal t where | ||
161 | real'' :: NumericContainer c => c Double -> c t | ||
162 | complex'' :: NumericContainer c => c t -> c (Complex Double) | ||
163 | |||
164 | |||
165 | instance AutoReal Double where | ||
166 | real'' = real | ||
167 | complex'' = complex | ||
168 | |||
169 | instance AutoReal (Complex Double) where | ||
170 | real'' = real | ||
171 | complex'' = complex | ||
172 | |||
173 | instance AutoReal Float where | ||
174 | real'' = real . single | ||
175 | complex'' = double . complex | ||
176 | |||
177 | instance AutoReal (Complex Float) where | ||
178 | real'' = real . single | ||
179 | complex'' = double . complex | ||
180 | |||
181 | ------------------------------------------------------------------- | ||
182 | |||
183 | -- | Basic element-by-element functions for numeric containers | 59 | -- | Basic element-by-element functions for numeric containers |
184 | class (Element e) => Container c e where | 60 | class (Element e) => Container c e where |
185 | {- | 61 | |
186 | -- | create a structure with a single element | 62 | -- | create a structure with a single element |
187 | scalar :: e -> c e | 63 | scalar :: e -> c e |
188 | -- | multiply every element by a scalar | ||
189 | scale :: e -> c e -> c e | 64 | scale :: e -> c e -> c e |
190 | -- | scale the element by element reciprocal of the object: | 65 | -- | scale the element by element reciprocal of the object: |
191 | -- | 66 | -- |
192 | -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@ | 67 | -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@ |
193 | scaleRecip :: e -> c e -> c e | 68 | scaleRecip :: e -> c e -> c e |
194 | -- | add a constant to each element | ||
195 | addConstant :: e -> c e -> c e | 69 | addConstant :: e -> c e -> c e |
196 | add :: c e -> c e -> c e | 70 | add :: c e -> c e -> c e |
197 | sub :: c e -> c e -> c e | 71 | sub :: c e -> c e -> c e |
@@ -200,7 +74,7 @@ class (Element e) => Container c e where | |||
200 | -- | element by element division | 74 | -- | element by element division |
201 | divide :: c e -> c e -> c e | 75 | divide :: c e -> c e -> c e |
202 | equal :: c e -> c e -> Bool | 76 | equal :: c e -> c e -> Bool |
203 | -} | 77 | |
204 | -- | cannot implement instance Functor because of Element class constraint | 78 | -- | cannot implement instance Functor because of Element class constraint |
205 | cmap :: (Element a, Element b) => (a -> b) -> c a -> c b | 79 | cmap :: (Element a, Element b) => (a -> b) -> c a -> c b |
206 | -- | 80 | -- |
@@ -217,23 +91,239 @@ class (Element e) => Container c e where | |||
217 | sumElements :: c e -> e | 91 | sumElements :: c e -> e |
218 | prodElements :: c e -> e | 92 | prodElements :: c e -> e |
219 | 93 | ||
220 | -- | Basic element-by-element functions. | 94 | -- -- | Basic element-by-element functions. |
221 | class (Element e, Container c e) => Linear c e where | 95 | -- class (Element e, Container c e) => Linear c e where |
222 | -- | create a structure with a single element | ||
223 | scalar :: e -> c e | ||
224 | scale :: e -> c e -> c e | ||
225 | -- | scale the element by element reciprocal of the object: | ||
226 | -- | ||
227 | -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@ | ||
228 | scaleRecip :: e -> c e -> c e | ||
229 | addConstant :: e -> c e -> c e | ||
230 | add :: c e -> c e -> c e | ||
231 | sub :: c e -> c e -> c e | ||
232 | -- | element by element multiplication | ||
233 | mul :: c e -> c e -> c e | ||
234 | -- | element by element division | ||
235 | divide :: c e -> c e -> c e | ||
236 | equal :: c e -> c e -> Bool | ||
237 | 96 | ||
238 | 97 | ||
98 | -------------------------------------------------------------------------- | ||
239 | 99 | ||
100 | instance Container Vector Float where | ||
101 | scale = vectorMapValF Scale | ||
102 | scaleRecip = vectorMapValF Recip | ||
103 | addConstant = vectorMapValF AddConstant | ||
104 | add = vectorZipF Add | ||
105 | sub = vectorZipF Sub | ||
106 | mul = vectorZipF Mul | ||
107 | divide = vectorZipF Div | ||
108 | equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 | ||
109 | scalar x = fromList [x] | ||
110 | -- | ||
111 | --instance Container Vector Float where | ||
112 | cmap = mapVector | ||
113 | atIndex = (@>) | ||
114 | minIndex = round . toScalarF MinIdx | ||
115 | maxIndex = round . toScalarF MaxIdx | ||
116 | minElement = toScalarF Min | ||
117 | maxElement = toScalarF Max | ||
118 | sumElements = sumF | ||
119 | prodElements = prodF | ||
120 | |||
121 | instance Container Vector Double where | ||
122 | scale = vectorMapValR Scale | ||
123 | scaleRecip = vectorMapValR Recip | ||
124 | addConstant = vectorMapValR AddConstant | ||
125 | add = vectorZipR Add | ||
126 | sub = vectorZipR Sub | ||
127 | mul = vectorZipR Mul | ||
128 | divide = vectorZipR Div | ||
129 | equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 | ||
130 | scalar x = fromList [x] | ||
131 | -- | ||
132 | --instance Container Vector Double where | ||
133 | cmap = mapVector | ||
134 | atIndex = (@>) | ||
135 | minIndex = round . toScalarR MinIdx | ||
136 | maxIndex = round . toScalarR MaxIdx | ||
137 | minElement = toScalarR Min | ||
138 | maxElement = toScalarR Max | ||
139 | sumElements = sumR | ||
140 | prodElements = prodR | ||
141 | |||
142 | instance Container Vector (Complex Double) where | ||
143 | scale = vectorMapValC Scale | ||
144 | scaleRecip = vectorMapValC Recip | ||
145 | addConstant = vectorMapValC AddConstant | ||
146 | add = vectorZipC Add | ||
147 | sub = vectorZipC Sub | ||
148 | mul = vectorZipC Mul | ||
149 | divide = vectorZipC Div | ||
150 | equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 | ||
151 | scalar x = fromList [x] | ||
152 | -- | ||
153 | --instance Container Vector (Complex Double) where | ||
154 | cmap = mapVector | ||
155 | atIndex = (@>) | ||
156 | minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) | ||
157 | maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) | ||
158 | minElement = ap (@>) minIndex | ||
159 | maxElement = ap (@>) maxIndex | ||
160 | sumElements = sumC | ||
161 | prodElements = prodC | ||
162 | |||
163 | instance Container Vector (Complex Float) where | ||
164 | scale = vectorMapValQ Scale | ||
165 | scaleRecip = vectorMapValQ Recip | ||
166 | addConstant = vectorMapValQ AddConstant | ||
167 | add = vectorZipQ Add | ||
168 | sub = vectorZipQ Sub | ||
169 | mul = vectorZipQ Mul | ||
170 | divide = vectorZipQ Div | ||
171 | equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 | ||
172 | scalar x = fromList [x] | ||
173 | -- | ||
174 | --instance Container Vector (Complex Float) where | ||
175 | cmap = mapVector | ||
176 | atIndex = (@>) | ||
177 | minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) | ||
178 | maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) | ||
179 | minElement = ap (@>) minIndex | ||
180 | maxElement = ap (@>) maxIndex | ||
181 | sumElements = sumQ | ||
182 | prodElements = prodQ | ||
183 | |||
184 | --------------------------------------------------------------- | ||
185 | |||
186 | instance (Container Vector a) => Container Matrix a where | ||
187 | scale x = liftMatrix (scale x) | ||
188 | scaleRecip x = liftMatrix (scaleRecip x) | ||
189 | addConstant x = liftMatrix (addConstant x) | ||
190 | add = liftMatrix2 add | ||
191 | sub = liftMatrix2 sub | ||
192 | mul = liftMatrix2 mul | ||
193 | divide = liftMatrix2 divide | ||
194 | equal a b = cols a == cols b && flatten a `equal` flatten b | ||
195 | scalar x = (1><1) [x] | ||
196 | -- | ||
197 | --instance (Container Vector a) => Container Matrix a where | ||
198 | cmap f = liftMatrix (mapVector f) | ||
199 | atIndex = (@@>) | ||
200 | minIndex m = let (r,c) = (rows m,cols m) | ||
201 | i = (minIndex $ flatten m) | ||
202 | in (i `div` c,(i `mod` c) + 1) | ||
203 | maxIndex m = let (r,c) = (rows m,cols m) | ||
204 | i = (maxIndex $ flatten m) | ||
205 | in (i `div` c,(i `mod` c) + 1) | ||
206 | minElement = ap (@@>) minIndex | ||
207 | maxElement = ap (@@>) maxIndex | ||
208 | sumElements = sumElements . flatten | ||
209 | prodElements = prodElements . flatten | ||
210 | |||
211 | ---------------------------------------------------- | ||
212 | |||
213 | |||
214 | -- | Linear algebraic properties of objects | ||
215 | class Num e => Vectors a e where | ||
216 | -- | dot (inner) product | ||
217 | dot :: a e -> a e -> e | ||
218 | -- | sum of absolute value of elements (differs in complex case from @norm1@ | ||
219 | absSum :: a e -> e | ||
220 | -- | sum of absolute value of elements | ||
221 | norm1 :: a e -> e | ||
222 | -- | euclidean norm | ||
223 | norm2 :: a e -> e | ||
224 | -- | element of maximum magnitude | ||
225 | normInf :: a e -> e | ||
226 | |||
227 | instance Vectors Vector Float where | ||
228 | norm2 = toScalarF Norm2 | ||
229 | absSum = toScalarF AbsSum | ||
230 | dot = dotF | ||
231 | norm1 = toScalarF AbsSum | ||
232 | normInf = maxElement . vectorMapF Abs | ||
233 | |||
234 | instance Vectors Vector Double where | ||
235 | norm2 = toScalarR Norm2 | ||
236 | absSum = toScalarR AbsSum | ||
237 | dot = dotR | ||
238 | norm1 = toScalarR AbsSum | ||
239 | normInf = maxElement . vectorMapR Abs | ||
240 | |||
241 | instance Vectors Vector (Complex Float) where | ||
242 | norm2 = (:+ 0) . toScalarQ Norm2 | ||
243 | absSum = (:+ 0) . toScalarQ AbsSum | ||
244 | dot = dotQ | ||
245 | norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapQ Abs | ||
246 | normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapQ Abs | ||
247 | |||
248 | instance Vectors Vector (Complex Double) where | ||
249 | norm2 = (:+ 0) . toScalarC Norm2 | ||
250 | absSum = (:+ 0) . toScalarC AbsSum | ||
251 | dot = dotC | ||
252 | norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapC Abs | ||
253 | normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapC Abs | ||
254 | |||
255 | ---------------------------------------------------- | ||
256 | |||
257 | class Element t => Product t where | ||
258 | multiply :: Matrix t -> Matrix t -> Matrix t | ||
259 | ctrans :: Matrix t -> Matrix t | ||
260 | |||
261 | instance Product Double where | ||
262 | multiply = multiplyR | ||
263 | ctrans = trans | ||
264 | |||
265 | instance Product (Complex Double) where | ||
266 | multiply = multiplyC | ||
267 | ctrans = conj . trans | ||
268 | |||
269 | instance Product Float where | ||
270 | multiply = multiplyF | ||
271 | ctrans = trans | ||
272 | |||
273 | instance Product (Complex Float) where | ||
274 | multiply = multiplyQ | ||
275 | ctrans = conj . trans | ||
276 | |||
277 | ---------------------------------------------------------- | ||
278 | |||
279 | -- synonym for matrix product | ||
280 | mXm :: Product t => Matrix t -> Matrix t -> Matrix t | ||
281 | mXm = multiply | ||
282 | |||
283 | -- matrix - vector product | ||
284 | mXv :: Product t => Matrix t -> Vector t -> Vector t | ||
285 | mXv m v = flatten $ m `mXm` (asColumn v) | ||
286 | |||
287 | -- vector - matrix product | ||
288 | vXm :: Product t => Vector t -> Matrix t -> Vector t | ||
289 | vXm v m = flatten $ (asRow v) `mXm` m | ||
290 | |||
291 | {- | Outer product of two vectors. | ||
292 | |||
293 | @\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3] | ||
294 | (3><3) | ||
295 | [ 5.0, 2.0, 3.0 | ||
296 | , 10.0, 4.0, 6.0 | ||
297 | , 15.0, 6.0, 9.0 ]@ | ||
298 | -} | ||
299 | outer :: (Product t) => Vector t -> Vector t -> Matrix t | ||
300 | outer u v = asColumn u `multiply` asRow v | ||
301 | |||
302 | {- | Kronecker product of two matrices. | ||
303 | |||
304 | @m1=(2><3) | ||
305 | [ 1.0, 2.0, 0.0 | ||
306 | , 0.0, -1.0, 3.0 ] | ||
307 | m2=(4><3) | ||
308 | [ 1.0, 2.0, 3.0 | ||
309 | , 4.0, 5.0, 6.0 | ||
310 | , 7.0, 8.0, 9.0 | ||
311 | , 10.0, 11.0, 12.0 ]@ | ||
312 | |||
313 | @\> kronecker m1 m2 | ||
314 | (8><9) | ||
315 | [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0 | ||
316 | , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0 | ||
317 | , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0 | ||
318 | , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0 | ||
319 | , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0 | ||
320 | , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0 | ||
321 | , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 | ||
322 | , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@ | ||
323 | -} | ||
324 | kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t | ||
325 | kronecker a b = fromBlocks | ||
326 | . splitEvery (cols a) | ||
327 | . map (reshape (cols b)) | ||
328 | . toRows | ||
329 | $ flatten a `outer` flatten b | ||