summaryrefslogtreecommitdiff
path: root/lib/Numeric/ContainerBoot.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/ContainerBoot.hs')
-rw-r--r--lib/Numeric/ContainerBoot.hs611
1 files changed, 0 insertions, 611 deletions
diff --git a/lib/Numeric/ContainerBoot.hs b/lib/Numeric/ContainerBoot.hs
deleted file mode 100644
index ea4262c..0000000
--- a/lib/Numeric/ContainerBoot.hs
+++ /dev/null
@@ -1,611 +0,0 @@
1{-# LANGUAGE CPP #-}
2{-# LANGUAGE TypeFamilies #-}
3{-# LANGUAGE FlexibleContexts #-}
4{-# LANGUAGE FlexibleInstances #-}
5{-# LANGUAGE MultiParamTypeClasses #-}
6{-# LANGUAGE UndecidableInstances #-}
7
8-----------------------------------------------------------------------------
9-- |
10-- Module : Numeric.ContainerBoot
11-- Copyright : (c) Alberto Ruiz 2010
12-- License : GPL-style
13--
14-- Maintainer : Alberto Ruiz <aruiz@um.es>
15-- Stability : provisional
16-- Portability : portable
17--
18-- Module to avoid cyclyc dependencies.
19--
20-----------------------------------------------------------------------------
21
22module Numeric.ContainerBoot (
23 -- * Basic functions
24 ident, diag, ctrans,
25 -- * Generic operations
26 Container(..),
27 -- * Matrix product and related functions
28 Product(..), udot,
29 mXm,mXv,vXm,
30 outer, kronecker,
31 -- * Element conversion
32 Convert(..),
33 Complexable(),
34 RealElement(),
35
36 RealOf, ComplexOf, SingleOf, DoubleOf,
37
38 IndexOf,
39 module Data.Complex
40) where
41
42import Data.Packed
43import Data.Packed.ST as ST
44import Numeric.Conversion
45import Data.Packed.Internal
46import Numeric.GSL.Vector
47import Data.Complex
48import Control.Applicative((<*>))
49
50import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
51
52-------------------------------------------------------------------
53
54type family IndexOf (c :: * -> *)
55
56type instance IndexOf Vector = Int
57type instance IndexOf Matrix = (Int,Int)
58
59type family ArgOf (c :: * -> *) a
60
61type instance ArgOf Vector a = a -> a
62type instance ArgOf Matrix a = a -> a -> a
63
64-------------------------------------------------------------------
65
66-- | Basic element-by-element functions for numeric containers
67class (Complexable c, Fractional e, Element e) => Container c e where
68 -- | create a structure with a single element
69 --
70 -- >>> let v = fromList [1..3::Double]
71 -- >>> v / scalar (norm2 v)
72 -- fromList [0.2672612419124244,0.5345224838248488,0.8017837257372732]
73 --
74 scalar :: e -> c e
75 -- | complex conjugate
76 conj :: c e -> c e
77 scale :: e -> c e -> c e
78 -- | scale the element by element reciprocal of the object:
79 --
80 -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
81 scaleRecip :: e -> c e -> c e
82 addConstant :: e -> c e -> c e
83 add :: c e -> c e -> c e
84 sub :: c e -> c e -> c e
85 -- | element by element multiplication
86 mul :: c e -> c e -> c e
87 -- | element by element division
88 divide :: c e -> c e -> c e
89 equal :: c e -> c e -> Bool
90 --
91 -- element by element inverse tangent
92 arctan2 :: c e -> c e -> c e
93 --
94 -- | cannot implement instance Functor because of Element class constraint
95 cmap :: (Element b) => (e -> b) -> c e -> c b
96 -- | constant structure of given size
97 konst' :: e -> IndexOf c -> c e
98 -- | create a structure using a function
99 --
100 -- Hilbert matrix of order N:
101 --
102 -- @hilb n = build' (n,n) (\\i j -> 1/(i+j+1))@
103 build' :: IndexOf c -> (ArgOf c e) -> c e
104 -- | indexing function
105 atIndex :: c e -> IndexOf c -> e
106 -- | index of min element
107 minIndex :: c e -> IndexOf c
108 -- | index of max element
109 maxIndex :: c e -> IndexOf c
110 -- | value of min element
111 minElement :: c e -> e
112 -- | value of max element
113 maxElement :: c e -> e
114 -- the C functions sumX/prodX are twice as fast as using foldVector
115 -- | the sum of elements (faster than using @fold@)
116 sumElements :: c e -> e
117 -- | the product of elements (faster than using @fold@)
118 prodElements :: c e -> e
119
120 -- | A more efficient implementation of @cmap (\\x -> if x>0 then 1 else 0)@
121 --
122 -- >>> step $ linspace 5 (-1,1::Double)
123 -- 5 |> [0.0,0.0,0.0,1.0,1.0]
124 --
125
126 step :: RealElement e => c e -> c e
127
128 -- | Element by element version of @case compare a b of {LT -> l; EQ -> e; GT -> g}@.
129 --
130 -- Arguments with any dimension = 1 are automatically expanded:
131 --
132 -- >>> cond ((1><4)[1..]) ((3><1)[1..]) 0 100 ((3><4)[1..]) :: Matrix Double
133 -- (3><4)
134 -- [ 100.0, 2.0, 3.0, 4.0
135 -- , 0.0, 100.0, 7.0, 8.0
136 -- , 0.0, 0.0, 100.0, 12.0 ]
137 --
138
139 cond :: RealElement e
140 => c e -- ^ a
141 -> c e -- ^ b
142 -> c e -- ^ l
143 -> c e -- ^ e
144 -> c e -- ^ g
145 -> c e -- ^ result
146
147 -- | Find index of elements which satisfy a predicate
148 --
149 -- >>> find (>0) (ident 3 :: Matrix Double)
150 -- [(0,0),(1,1),(2,2)]
151 --
152
153 find :: (e -> Bool) -> c e -> [IndexOf c]
154
155 -- | Create a structure from an association list
156 --
157 -- >>> assoc 5 0 [(3,7),(1,4)] :: Vector Double
158 -- fromList [0.0,4.0,0.0,7.0,0.0]
159 --
160 -- >>> assoc (2,3) 0 [((0,2),7),((1,0),2*i-3)] :: Matrix (Complex Double)
161 -- (2><3)
162 -- [ 0.0 :+ 0.0, 0.0 :+ 0.0, 7.0 :+ 0.0
163 -- , (-3.0) :+ 2.0, 0.0 :+ 0.0, 0.0 :+ 0.0 ]
164 --
165 assoc :: IndexOf c -- ^ size
166 -> e -- ^ default value
167 -> [(IndexOf c, e)] -- ^ association list
168 -> c e -- ^ result
169
170 -- | Modify a structure using an update function
171 --
172 -- >>> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix Double
173 -- (5><5)
174 -- [ 1.0, 0.0, 0.0, 3.0, 0.0
175 -- , 0.0, 6.0, 0.0, 0.0, 0.0
176 -- , 0.0, 0.0, 1.0, 0.0, 0.0
177 -- , 0.0, 0.0, 0.0, 1.0, 0.0
178 -- , 0.0, 0.0, 0.0, 0.0, 1.0 ]
179 --
180 -- computation of histogram:
181 --
182 -- >>> accum (konst 0 7) (+) (map (flip (,) 1) [4,5,4,1,5,2,5]) :: Vector Double
183 -- fromList [0.0,1.0,1.0,0.0,2.0,3.0,0.0]
184 --
185
186 accum :: c e -- ^ initial structure
187 -> (e -> e -> e) -- ^ update function
188 -> [(IndexOf c, e)] -- ^ association list
189 -> c e -- ^ result
190
191--------------------------------------------------------------------------
192
193instance Container Vector Float where
194 scale = vectorMapValF Scale
195 scaleRecip = vectorMapValF Recip
196 addConstant = vectorMapValF AddConstant
197 add = vectorZipF Add
198 sub = vectorZipF Sub
199 mul = vectorZipF Mul
200 divide = vectorZipF Div
201 equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
202 arctan2 = vectorZipF ATan2
203 scalar x = fromList [x]
204 konst' = constantD
205 build' = buildV
206 conj = id
207 cmap = mapVector
208 atIndex = (@>)
209 minIndex = emptyErrorV "minIndex" (round . toScalarF MinIdx)
210 maxIndex = emptyErrorV "maxIndex" (round . toScalarF MaxIdx)
211 minElement = emptyErrorV "minElement" (toScalarF Min)
212 maxElement = emptyErrorV "maxElement" (toScalarF Max)
213 sumElements = sumF
214 prodElements = prodF
215 step = stepF
216 find = findV
217 assoc = assocV
218 accum = accumV
219 cond = condV condF
220
221instance Container Vector Double where
222 scale = vectorMapValR Scale
223 scaleRecip = vectorMapValR Recip
224 addConstant = vectorMapValR AddConstant
225 add = vectorZipR Add
226 sub = vectorZipR Sub
227 mul = vectorZipR Mul
228 divide = vectorZipR Div
229 equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
230 arctan2 = vectorZipR ATan2
231 scalar x = fromList [x]
232 konst' = constantD
233 build' = buildV
234 conj = id
235 cmap = mapVector
236 atIndex = (@>)
237 minIndex = emptyErrorV "minIndex" (round . toScalarR MinIdx)
238 maxIndex = emptyErrorV "maxIndex" (round . toScalarR MaxIdx)
239 minElement = emptyErrorV "minElement" (toScalarR Min)
240 maxElement = emptyErrorV "maxElement" (toScalarR Max)
241 sumElements = sumR
242 prodElements = prodR
243 step = stepD
244 find = findV
245 assoc = assocV
246 accum = accumV
247 cond = condV condD
248
249instance Container Vector (Complex Double) where
250 scale = vectorMapValC Scale
251 scaleRecip = vectorMapValC Recip
252 addConstant = vectorMapValC AddConstant
253 add = vectorZipC Add
254 sub = vectorZipC Sub
255 mul = vectorZipC Mul
256 divide = vectorZipC Div
257 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
258 arctan2 = vectorZipC ATan2
259 scalar x = fromList [x]
260 konst' = constantD
261 build' = buildV
262 conj = conjugateC
263 cmap = mapVector
264 atIndex = (@>)
265 minIndex = emptyErrorV "minIndex" (minIndex . fst . fromComplex . (mul <*> conj))
266 maxIndex = emptyErrorV "maxIndex" (maxIndex . fst . fromComplex . (mul <*> conj))
267 minElement = emptyErrorV "minElement" (atIndex <*> minIndex)
268 maxElement = emptyErrorV "maxElement" (atIndex <*> maxIndex)
269 sumElements = sumC
270 prodElements = prodC
271 step = undefined -- cannot match
272 find = findV
273 assoc = assocV
274 accum = accumV
275 cond = undefined -- cannot match
276
277instance Container Vector (Complex Float) where
278 scale = vectorMapValQ Scale
279 scaleRecip = vectorMapValQ Recip
280 addConstant = vectorMapValQ AddConstant
281 add = vectorZipQ Add
282 sub = vectorZipQ Sub
283 mul = vectorZipQ Mul
284 divide = vectorZipQ Div
285 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
286 arctan2 = vectorZipQ ATan2
287 scalar x = fromList [x]
288 konst' = constantD
289 build' = buildV
290 conj = conjugateQ
291 cmap = mapVector
292 atIndex = (@>)
293 minIndex = emptyErrorV "minIndex" (minIndex . fst . fromComplex . (mul <*> conj))
294 maxIndex = emptyErrorV "maxIndex" (maxIndex . fst . fromComplex . (mul <*> conj))
295 minElement = emptyErrorV "minElement" (atIndex <*> minIndex)
296 maxElement = emptyErrorV "maxElement" (atIndex <*> maxIndex)
297 sumElements = sumQ
298 prodElements = prodQ
299 step = undefined -- cannot match
300 find = findV
301 assoc = assocV
302 accum = accumV
303 cond = undefined -- cannot match
304
305---------------------------------------------------------------
306
307instance (Container Vector a) => Container Matrix a where
308 scale x = liftMatrix (scale x)
309 scaleRecip x = liftMatrix (scaleRecip x)
310 addConstant x = liftMatrix (addConstant x)
311 add = liftMatrix2 add
312 sub = liftMatrix2 sub
313 mul = liftMatrix2 mul
314 divide = liftMatrix2 divide
315 equal a b = cols a == cols b && flatten a `equal` flatten b
316 arctan2 = liftMatrix2 arctan2
317 scalar x = (1><1) [x]
318 konst' v (r,c) = matrixFromVector RowMajor r c (konst' v (r*c))
319 build' = buildM
320 conj = liftMatrix conj
321 cmap f = liftMatrix (mapVector f)
322 atIndex = (@@>)
323 minIndex = emptyErrorM "minIndex of Matrix" $
324 \m -> divMod (minIndex $ flatten m) (cols m)
325 maxIndex = emptyErrorM "maxIndex of Matrix" $
326 \m -> divMod (maxIndex $ flatten m) (cols m)
327 minElement = emptyErrorM "minElement of Matrix" (atIndex <*> minIndex)
328 maxElement = emptyErrorM "maxElement of Matrix" (atIndex <*> maxIndex)
329 sumElements = sumElements . flatten
330 prodElements = prodElements . flatten
331 step = liftMatrix step
332 find = findM
333 assoc = assocM
334 accum = accumM
335 cond = condM
336
337
338emptyErrorV msg f v =
339 if dim v > 0
340 then f v
341 else error $ msg ++ " of Vector with dim = 0"
342
343emptyErrorM msg f m =
344 if rows m > 0 && cols m > 0
345 then f m
346 else error $ msg++" "++shSize m
347
348----------------------------------------------------
349
350-- | Matrix product and related functions
351class (Num e, Element e) => Product e where
352 -- | matrix product
353 multiply :: Matrix e -> Matrix e -> Matrix e
354 -- | sum of absolute value of elements (differs in complex case from @norm1@)
355 absSum :: Vector e -> RealOf e
356 -- | sum of absolute value of elements
357 norm1 :: Vector e -> RealOf e
358 -- | euclidean norm
359 norm2 :: Vector e -> RealOf e
360 -- | element of maximum magnitude
361 normInf :: Vector e -> RealOf e
362
363instance Product Float where
364 norm2 = emptyVal (toScalarF Norm2)
365 absSum = emptyVal (toScalarF AbsSum)
366 norm1 = emptyVal (toScalarF AbsSum)
367 normInf = emptyVal (maxElement . vectorMapF Abs)
368 multiply = emptyMul multiplyF
369
370instance Product Double where
371 norm2 = emptyVal (toScalarR Norm2)
372 absSum = emptyVal (toScalarR AbsSum)
373 norm1 = emptyVal (toScalarR AbsSum)
374 normInf = emptyVal (maxElement . vectorMapR Abs)
375 multiply = emptyMul multiplyR
376
377instance Product (Complex Float) where
378 norm2 = emptyVal (toScalarQ Norm2)
379 absSum = emptyVal (toScalarQ AbsSum)
380 norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapQ Abs)
381 normInf = emptyVal (maxElement . fst . fromComplex . vectorMapQ Abs)
382 multiply = emptyMul multiplyQ
383
384instance Product (Complex Double) where
385 norm2 = emptyVal (toScalarC Norm2)
386 absSum = emptyVal (toScalarC AbsSum)
387 norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapC Abs)
388 normInf = emptyVal (maxElement . fst . fromComplex . vectorMapC Abs)
389 multiply = emptyMul multiplyC
390
391emptyMul m a b
392 | x1 == 0 && x2 == 0 || r == 0 || c == 0 = konst' 0 (r,c)
393 | otherwise = m a b
394 where
395 r = rows a
396 x1 = cols a
397 x2 = rows b
398 c = cols b
399
400emptyVal f v =
401 if dim v > 0
402 then f v
403 else 0
404
405-- FIXME remove unused C wrappers
406-- | (unconjugated) dot product
407udot :: Product e => Vector e -> Vector e -> e
408udot u v
409 | dim u == dim v = val (asRow u `multiply` asColumn v)
410 | otherwise = error $ "different dimensions "++show (dim u)++" and "++show (dim v)++" in dot product"
411 where
412 val m | dim u > 0 = m@@>(0,0)
413 | otherwise = 0
414
415----------------------------------------------------------
416
417-- synonym for matrix product
418mXm :: Product t => Matrix t -> Matrix t -> Matrix t
419mXm = multiply
420
421-- matrix - vector product
422mXv :: Product t => Matrix t -> Vector t -> Vector t
423mXv m v = flatten $ m `mXm` (asColumn v)
424
425-- vector - matrix product
426vXm :: Product t => Vector t -> Matrix t -> Vector t
427vXm v m = flatten $ (asRow v) `mXm` m
428
429{- | Outer product of two vectors.
430
431>>> fromList [1,2,3] `outer` fromList [5,2,3]
432(3><3)
433 [ 5.0, 2.0, 3.0
434 , 10.0, 4.0, 6.0
435 , 15.0, 6.0, 9.0 ]
436
437-}
438outer :: (Product t) => Vector t -> Vector t -> Matrix t
439outer u v = asColumn u `multiply` asRow v
440
441{- | Kronecker product of two matrices.
442
443@m1=(2><3)
444 [ 1.0, 2.0, 0.0
445 , 0.0, -1.0, 3.0 ]
446m2=(4><3)
447 [ 1.0, 2.0, 3.0
448 , 4.0, 5.0, 6.0
449 , 7.0, 8.0, 9.0
450 , 10.0, 11.0, 12.0 ]@
451
452>>> kronecker m1 m2
453(8><9)
454 [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0
455 , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0
456 , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0
457 , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0
458 , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0
459 , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0
460 , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0
461 , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]
462
463-}
464kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
465kronecker a b = fromBlocks
466 . splitEvery (cols a)
467 . map (reshape (cols b))
468 . toRows
469 $ flatten a `outer` flatten b
470
471-------------------------------------------------------------------
472
473
474class Convert t where
475 real :: Container c t => c (RealOf t) -> c t
476 complex :: Container c t => c t -> c (ComplexOf t)
477 single :: Container c t => c t -> c (SingleOf t)
478 double :: Container c t => c t -> c (DoubleOf t)
479 toComplex :: (Container c t, RealElement t) => (c t, c t) -> c (Complex t)
480 fromComplex :: (Container c t, RealElement t) => c (Complex t) -> (c t, c t)
481
482
483instance Convert Double where
484 real = id
485 complex = comp'
486 single = single'
487 double = id
488 toComplex = toComplex'
489 fromComplex = fromComplex'
490
491instance Convert Float where
492 real = id
493 complex = comp'
494 single = id
495 double = double'
496 toComplex = toComplex'
497 fromComplex = fromComplex'
498
499instance Convert (Complex Double) where
500 real = comp'
501 complex = id
502 single = single'
503 double = id
504 toComplex = toComplex'
505 fromComplex = fromComplex'
506
507instance Convert (Complex Float) where
508 real = comp'
509 complex = id
510 single = id
511 double = double'
512 toComplex = toComplex'
513 fromComplex = fromComplex'
514
515-------------------------------------------------------------------
516
517type family RealOf x
518
519type instance RealOf Double = Double
520type instance RealOf (Complex Double) = Double
521
522type instance RealOf Float = Float
523type instance RealOf (Complex Float) = Float
524
525type family ComplexOf x
526
527type instance ComplexOf Double = Complex Double
528type instance ComplexOf (Complex Double) = Complex Double
529
530type instance ComplexOf Float = Complex Float
531type instance ComplexOf (Complex Float) = Complex Float
532
533type family SingleOf x
534
535type instance SingleOf Double = Float
536type instance SingleOf Float = Float
537
538type instance SingleOf (Complex a) = Complex (SingleOf a)
539
540type family DoubleOf x
541
542type instance DoubleOf Double = Double
543type instance DoubleOf Float = Double
544
545type instance DoubleOf (Complex a) = Complex (DoubleOf a)
546
547type family ElementOf c
548
549type instance ElementOf (Vector a) = a
550type instance ElementOf (Matrix a) = a
551
552------------------------------------------------------------
553
554buildM (rc,cc) f = fromLists [ [f r c | c <- cs] | r <- rs ]
555 where rs = map fromIntegral [0 .. (rc-1)]
556 cs = map fromIntegral [0 .. (cc-1)]
557
558buildV n f = fromList [f k | k <- ks]
559 where ks = map fromIntegral [0 .. (n-1)]
560
561--------------------------------------------------------
562-- | conjugate transpose
563ctrans :: (Container Vector e, Element e) => Matrix e -> Matrix e
564ctrans = liftMatrix conj . trans
565
566-- | Creates a square matrix with a given diagonal.
567diag :: (Num a, Element a) => Vector a -> Matrix a
568diag v = diagRect 0 v n n where n = dim v
569
570-- | creates the identity matrix of given dimension
571ident :: (Num a, Element a) => Int -> Matrix a
572ident n = diag (constantD 1 n)
573
574--------------------------------------------------------
575
576findV p x = foldVectorWithIndex g [] x where
577 g k z l = if p z then k:l else l
578
579findM p x = map ((`divMod` cols x)) $ findV p (flatten x)
580
581assocV n z xs = ST.runSTVector $ do
582 v <- ST.newVector z n
583 mapM_ (\(k,x) -> ST.writeVector v k x) xs
584 return v
585
586assocM (r,c) z xs = ST.runSTMatrix $ do
587 m <- ST.newMatrix z r c
588 mapM_ (\((i,j),x) -> ST.writeMatrix m i j x) xs
589 return m
590
591accumV v0 f xs = ST.runSTVector $ do
592 v <- ST.thawVector v0
593 mapM_ (\(k,x) -> ST.modifyVector v k (f x)) xs
594 return v
595
596accumM m0 f xs = ST.runSTMatrix $ do
597 m <- ST.thawMatrix m0
598 mapM_ (\((i,j),x) -> ST.modifyMatrix m i j (f x)) xs
599 return m
600
601----------------------------------------------------------------------
602
603condM a b l e t = matrixFromVector RowMajor (rows a'') (cols a'') $ cond a' b' l' e' t'
604 where
605 args@(a'':_) = conformMs [a,b,l,e,t]
606 [a', b', l', e', t'] = map flatten args
607
608condV f a b l e t = f a' b' l' e' t'
609 where
610 [a', b', l', e', t'] = conformVs [a,b,l,e,t]
611
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-15 09:25:46 +0000