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