summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/Data/Packed/Internal/Matrix.hs2
-rw-r--r--lib/Numeric/Container.hs103
-rw-r--r--lib/Numeric/Conversion.hs137
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs2
-rw-r--r--lib/Numeric/LinearAlgebra/LAPACK.hs8
-rw-r--r--lib/Numeric/LinearAlgebra/Tests.hs2
-rw-r--r--lib/Numeric/LinearAlgebra/Tests/Instances.hs12
7 files changed, 125 insertions, 141 deletions
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index d39481d..7a17ef0 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -245,7 +245,7 @@ compat m1 m2 = rows m1 == rows m2 && cols m1 == cols m2
245 245
246------------------------------------------------------------------ 246------------------------------------------------------------------
247 247
248-- | Auxiliary class. 248-- | Supported element types for basic matrix operations.
249class (Storable a, Floating a) => Element a where 249class (Storable a, Floating a) => Element a where
250 subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position 250 subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position
251 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix 251 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs
index aaa068f..e9a8f22 100644
--- a/lib/Numeric/Container.hs
+++ b/lib/Numeric/Container.hs
@@ -25,9 +25,10 @@ module Numeric.Container (
25 mXm,mXv,vXm, 25 mXm,mXv,vXm,
26 outer, kronecker, 26 outer, kronecker,
27 27
28 RealElement, --Precision, 28 Convert(..),
29 ComplexContainer(toComplex,fromComplex,comp,conj), 29 Complexable(),
30 Convert(..), --AutoReal(..), 30 RealElement(),
31
31 RealOf, ComplexOf, SingleOf, DoubleOf, 32 RealOf, ComplexOf, SingleOf, DoubleOf,
32 33
33 IndexOf, 34 IndexOf,
@@ -54,10 +55,11 @@ type instance IndexOf Matrix = (Int,Int)
54------------------------------------------------------------------- 55-------------------------------------------------------------------
55 56
56-- | Basic element-by-element functions for numeric containers 57-- | Basic element-by-element functions for numeric containers
57class (Element e) => Container c e where 58class (Complexable c, Element e) => Container c e where
58
59 -- | create a structure with a single element 59 -- | create a structure with a single element
60 scalar :: e -> c e 60 scalar :: e -> c e
61 -- | complex conjugate
62 conj :: c e -> c e
61 scale :: e -> c e -> c e 63 scale :: e -> c e -> c e
62 -- | scale the element by element reciprocal of the object: 64 -- | scale the element by element reciprocal of the object:
63 -- 65 --
@@ -75,7 +77,7 @@ class (Element e) => Container c e where
75 -- | cannot implement instance Functor because of Element class constraint 77 -- | cannot implement instance Functor because of Element class constraint
76 cmap :: (Element a, Element b) => (a -> b) -> c a -> c b 78 cmap :: (Element a, Element b) => (a -> b) -> c a -> c b
77 -- | constant structure of given size 79 -- | constant structure of given size
78 konst :: e -> IndexOf c -> c e 80 konst :: e -> IndexOf c -> c e
79 -- 81 --
80 -- | indexing function 82 -- | indexing function
81 atIndex :: c e -> IndexOf c -> e 83 atIndex :: c e -> IndexOf c -> e
@@ -110,6 +112,7 @@ instance Container Vector Float where
110 equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 112 equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
111 scalar x = fromList [x] 113 scalar x = fromList [x]
112 konst = constantD 114 konst = constantD
115 conj = conjugateD
113 cmap = mapVector 116 cmap = mapVector
114 atIndex = (@>) 117 atIndex = (@>)
115 minIndex = round . toScalarF MinIdx 118 minIndex = round . toScalarF MinIdx
@@ -130,6 +133,7 @@ instance Container Vector Double where
130 equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 133 equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
131 scalar x = fromList [x] 134 scalar x = fromList [x]
132 konst = constantD 135 konst = constantD
136 conj = conjugateD
133 cmap = mapVector 137 cmap = mapVector
134 atIndex = (@>) 138 atIndex = (@>)
135 minIndex = round . toScalarR MinIdx 139 minIndex = round . toScalarR MinIdx
@@ -150,6 +154,7 @@ instance Container Vector (Complex Double) where
150 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 154 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
151 scalar x = fromList [x] 155 scalar x = fromList [x]
152 konst = constantD 156 konst = constantD
157 conj = conjugateD
153 cmap = mapVector 158 cmap = mapVector
154 atIndex = (@>) 159 atIndex = (@>)
155 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) 160 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
@@ -170,6 +175,7 @@ instance Container Vector (Complex Float) where
170 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 175 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
171 scalar x = fromList [x] 176 scalar x = fromList [x]
172 konst = constantD 177 konst = constantD
178 conj = conjugateD
173 cmap = mapVector 179 cmap = mapVector
174 atIndex = (@>) 180 atIndex = (@>)
175 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) 181 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
@@ -192,6 +198,7 @@ instance (Container Vector a) => Container Matrix a where
192 equal a b = cols a == cols b && flatten a `equal` flatten b 198 equal a b = cols a == cols b && flatten a `equal` flatten b
193 scalar x = (1><1) [x] 199 scalar x = (1><1) [x]
194 konst v (r,c) = reshape c (konst v (r*c)) 200 konst v (r,c) = reshape c (konst v (r*c))
201 conj = liftMatrix conjugateD
195 cmap f = liftMatrix (mapVector f) 202 cmap f = liftMatrix (mapVector f)
196 atIndex = (@@>) 203 atIndex = (@@>)
197 minIndex m = let (r,c) = (rows m,cols m) 204 minIndex m = let (r,c) = (rows m,cols m)
@@ -208,7 +215,7 @@ instance (Container Vector a) => Container Matrix a where
208---------------------------------------------------- 215----------------------------------------------------
209 216
210 217
211-- | Linear algebraic properties of objects 218-- | Matrix product and related functions
212class Element e => Product e where 219class Element e => Product e where
213 -- | matrix product 220 -- | matrix product
214 multiply :: Matrix e -> Matrix e -> Matrix e 221 multiply :: Matrix e -> Matrix e -> Matrix e
@@ -309,4 +316,84 @@ kronecker a b = fromBlocks
309 . toRows 316 . toRows
310 $ flatten a `outer` flatten b 317 $ flatten a `outer` flatten b
311 318
312---------------------------------------------------------- 319-------------------------------------------------------------------
320
321
322class Convert t where
323 real :: Container c t => c (RealOf t) -> c t
324 complex :: Container c t => c t -> c (ComplexOf t)
325 single :: Container c t => c t -> c (SingleOf t)
326 double :: Container c t => c t -> c (DoubleOf t)
327 toComplex :: (Container c t, RealElement t) => (c t, c t) -> c (Complex t)
328 fromComplex :: (Container c t, RealElement t) => c (Complex t) -> (c t, c t)
329
330
331instance Convert Double where
332 real = id
333 complex = comp'
334 single = single'
335 double = id
336 toComplex = toComplex'
337 fromComplex = fromComplex'
338
339instance Convert Float where
340 real = id
341 complex = comp'
342 single = id
343 double = double'
344 toComplex = toComplex'
345 fromComplex = fromComplex'
346
347instance Convert (Complex Double) where
348 real = comp'
349 complex = id
350 single = single'
351 double = id
352 toComplex = toComplex'
353 fromComplex = fromComplex'
354
355instance Convert (Complex Float) where
356 real = comp'
357 complex = id
358 single = id
359 double = double'
360 toComplex = toComplex'
361 fromComplex = fromComplex'
362
363-------------------------------------------------------------------
364
365type family RealOf x
366
367type instance RealOf Double = Double
368type instance RealOf (Complex Double) = Double
369
370type instance RealOf Float = Float
371type instance RealOf (Complex Float) = Float
372
373type family ComplexOf x
374
375type instance ComplexOf Double = Complex Double
376type instance ComplexOf (Complex Double) = Complex Double
377
378type instance ComplexOf Float = Complex Float
379type instance ComplexOf (Complex Float) = Complex Float
380
381type family SingleOf x
382
383type instance SingleOf Double = Float
384type instance SingleOf Float = Float
385
386type instance SingleOf (Complex a) = Complex (SingleOf a)
387
388type family DoubleOf x
389
390type instance DoubleOf Double = Double
391type instance DoubleOf Float = Double
392
393type instance DoubleOf (Complex a) = Complex (DoubleOf a)
394
395type family ElementOf c
396
397type instance ElementOf (Vector a) = a
398type instance ElementOf (Matrix a) = a
399
diff --git a/lib/Numeric/Conversion.hs b/lib/Numeric/Conversion.hs
index 809ac51..fbf608a 100644
--- a/lib/Numeric/Conversion.hs
+++ b/lib/Numeric/Conversion.hs
@@ -20,24 +20,15 @@
20----------------------------------------------------------------------------- 20-----------------------------------------------------------------------------
21 21
22module Numeric.Conversion ( 22module Numeric.Conversion (
23 RealElement, --Precision, 23 Complexable(..), RealElement,
24 ComplexContainer(toComplex,fromComplex,conj,comp),
25 Convert(..), --AutoReal(..),
26 RealOf, ComplexOf, SingleOf, DoubleOf,
27 module Data.Complex 24 module Data.Complex
28) where 25) where
29 26
30
31import Data.Packed.Internal.Vector 27import Data.Packed.Internal.Vector
32import Data.Packed.Internal.Matrix 28import Data.Packed.Internal.Matrix
33--import Numeric.GSL.Vector
34
35import Data.Complex 29import Data.Complex
36--import Control.Monad(ap)
37import Control.Arrow((***)) 30import Control.Arrow((***))
38 31
39--import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
40
41------------------------------------------------------------------- 32-------------------------------------------------------------------
42 33
43-- | Supported single-double precision type pairs 34-- | Supported single-double precision type pairs
@@ -60,24 +51,22 @@ class (Element t, Element (Complex t), RealFloat t
60 => RealElement t 51 => RealElement t
61 52
62instance RealElement Double 53instance RealElement Double
63
64instance RealElement Float 54instance RealElement Float
65 55
66-- | Conversion utilities 56
67class ComplexContainer c where 57-- | Structures that may contain complex numbers
68 toComplex :: (RealElement e) => (c e, c e) -> c (Complex e) 58class Complexable c where
69 fromComplex :: (RealElement e) => c (Complex e) -> (c e, c e) 59 toComplex' :: (RealElement e) => (c e, c e) -> c (Complex e)
70 comp :: (RealElement e) => c e -> c (Complex e) 60 fromComplex' :: (RealElement e) => c (Complex e) -> (c e, c e)
71 conj :: (RealElement e) => c (Complex e) -> c (Complex e) 61 comp' :: (RealElement e) => c e -> c (Complex e)
72 single' :: Precision a b => c b -> c a 62 single' :: Precision a b => c b -> c a
73 double' :: Precision a b => c a -> c b 63 double' :: Precision a b => c a -> c b
74 64
75 65
76instance ComplexContainer Vector where 66instance Complexable Vector where
77 toComplex = toComplexV 67 toComplex' = toComplexV
78 fromComplex = fromComplexV 68 fromComplex' = fromComplexV
79 comp v = toComplex (v,constantD 0 (dim v)) 69 comp' v = toComplex' (v,constantD 0 (dim v))
80 conj = conjugateD
81 single' = double2FloatG 70 single' = double2FloatG
82 double' = float2DoubleG 71 double' = float2DoubleG
83 72
@@ -92,107 +81,11 @@ fromComplexV z = (r,i) where
92 [r,i] = toColumns $ reshape 2 $ asReal z 81 [r,i] = toColumns $ reshape 2 $ asReal z
93 82
94 83
95instance ComplexContainer Matrix where 84instance Complexable Matrix where
96 toComplex = uncurry $ liftMatrix2 $ curry toComplex 85 toComplex' = uncurry $ liftMatrix2 $ curry toComplex'
97 fromComplex z = (reshape c *** reshape c) . fromComplex . flatten $ z 86 fromComplex' z = (reshape c *** reshape c) . fromComplex' . flatten $ z
98 where c = cols z 87 where c = cols z
99 comp = liftMatrix comp 88 comp' = liftMatrix comp'
100 conj = liftMatrix conj
101 single' = liftMatrix single' 89 single' = liftMatrix single'
102 double' = liftMatrix double' 90 double' = liftMatrix double'
103 91
104-------------------------------------------------------------------
105
106type family RealOf x
107
108type instance RealOf Double = Double
109type instance RealOf (Complex Double) = Double
110
111type instance RealOf Float = Float
112type instance RealOf (Complex Float) = Float
113
114type family ComplexOf x
115
116type instance ComplexOf Double = Complex Double
117type instance ComplexOf (Complex Double) = Complex Double
118
119type instance ComplexOf Float = Complex Float
120type instance ComplexOf (Complex Float) = Complex Float
121
122type family SingleOf x
123
124type instance SingleOf Double = Float
125type instance SingleOf Float = Float
126
127type instance SingleOf (Complex a) = Complex (SingleOf a)
128
129type family DoubleOf x
130
131type instance DoubleOf Double = Double
132type instance DoubleOf Float = Double
133
134type instance DoubleOf (Complex a) = Complex (DoubleOf a)
135
136type family ElementOf c
137
138type instance ElementOf (Vector a) = a
139type instance ElementOf (Matrix a) = a
140
141
142-------------------------------------------------------------------
143
144class (Element t, Element (RealOf t)) => Convert t where
145 real :: ComplexContainer c => c (RealOf t) -> c t
146 complex :: ComplexContainer c => c t -> c (ComplexOf t)
147 single :: ComplexContainer c => c t -> c (SingleOf t)
148 double :: ComplexContainer c => c t -> c (DoubleOf t)
149
150
151instance Convert Double where
152 real = id
153 complex = comp
154 single = single'
155 double = id
156
157instance Convert Float where
158 real = id
159 complex = comp
160 single = id
161 double = double'
162
163instance Convert (Complex Double) where
164 real = comp
165 complex = id
166 single = single'
167 double = id
168
169instance Convert (Complex Float) where
170 real = comp
171 complex = id
172 single = id
173 double = double'
174
175-------------------------------------------------------------------
176
177-- | to be replaced by Convert
178class Convert t => AutoReal t where
179 real'' :: ComplexContainer c => c Double -> c t
180 complex'' :: ComplexContainer c => c t -> c (Complex Double)
181
182
183instance AutoReal Double where
184 real'' = real
185 complex'' = complex
186
187instance AutoReal (Complex Double) where
188 real'' = real
189 complex'' = complex
190
191instance AutoReal Float where
192 real'' = real . single
193 complex'' = double . complex
194
195instance AutoReal (Complex Float) where
196 real'' = real . single
197 complex'' = double . complex
198
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index f2f5473..3cd200e 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -82,7 +82,7 @@ import Data.Array
82import Numeric.Container 82import Numeric.Container
83 83
84-- | Auxiliary typeclass used to define generic computations for both real and complex matrices. 84-- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
85class (Product t, Container Vector t, Container Matrix t) => Field t where 85class (Product t, Convert t, Container Vector t, Container Matrix t) => Field t where
86 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) 86 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
87 thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) 87 thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
88 sv' :: Matrix t -> Vector Double 88 sv' :: Matrix t -> Vector Double
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index cb48571..fbc5460 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -259,14 +259,14 @@ eigRaux m = unsafePerformIO $ do
259 where r = rows m 259 where r = rows m
260 g ra ca pa = dgeev ra ca pa 0 0 nullPtr 260 g ra ca pa = dgeev ra ca pa 0 0 nullPtr
261 261
262fixeig1 s = toComplex (subVector 0 r (asReal s), subVector r r (asReal s)) 262fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s))
263 where r = dim s 263 where r = dim s
264 264
265fixeig [] _ = [] 265fixeig [] _ = []
266fixeig [_] [v] = [comp v] 266fixeig [_] [v] = [comp' v]
267fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) 267fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
268 | r1 == r2 && i1 == (-i2) = toComplex (v1,v2) : toComplex (v1,scale (-1) v2) : fixeig r vs 268 | r1 == r2 && i1 == (-i2) = toComplex' (v1,v2) : toComplex' (v1,scale (-1) v2) : fixeig r vs
269 | otherwise = comp v1 : fixeig ((r2:+i2):r) (v2:vs) 269 | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs)
270 where scale = vectorMapValR Scale 270 where scale = vectorMapValR Scale
271fixeig _ _ = error "fixeig with impossible inputs" 271fixeig _ _ = error "fixeig with impossible inputs"
272 272
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index aa7b01c..0df29a8 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -453,7 +453,7 @@ runTests n = do
453 , utest "randomGaussian" randomTestGaussian 453 , utest "randomGaussian" randomTestGaussian
454 , utest "randomUniform" randomTestUniform 454 , utest "randomUniform" randomTestUniform
455 , utest "buildVector/Matrix" $ 455 , utest "buildVector/Matrix" $
456 comp (10 |> [0::Double ..]) == buildVector 10 fromIntegral 456 complex (10 |> [0::Double ..]) == buildVector 10 fromIntegral
457 && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0) 457 && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0)
458 , utest "rank" $ rank ((2><3)[1,0,0,1,6*eps,0]) == 1 458 , utest "rank" $ rank ((2><3)[1,0,0,1,6*eps,0]) == 1
459 && rank ((2><3)[1,0,0,1,7*eps,0]) == 2 459 && rank ((2><3)[1,0,0,1,7*eps,0]) == 2
diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
index 6046ccb..804c481 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -1,4 +1,4 @@
1{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP #-} 1{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-}
2{-# OPTIONS_GHC -fno-warn-unused-imports #-} 2{-# OPTIONS_GHC -fno-warn-unused-imports #-}
3----------------------------------------------------------------------------- 3-----------------------------------------------------------------------------
4{- | 4{- |
@@ -135,10 +135,14 @@ instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where
135 coarbitrary = undefined 135 coarbitrary = undefined
136#endif 136#endif
137 137
138class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a
139instance ArbitraryField Double
140instance ArbitraryField (Complex Double)
141
138 142
139-- a well-conditioned general matrix (the singular values are between 1 and 100) 143-- a well-conditioned general matrix (the singular values are between 1 and 100)
140newtype (WC a) = WC (Matrix a) deriving Show 144newtype (WC a) = WC (Matrix a) deriving Show
141instance (Convert a, Field a, Arbitrary a, Random (RealOf a)) => Arbitrary (WC a) where 145instance (ArbitraryField a) => Arbitrary (WC a) where
142 arbitrary = do 146 arbitrary = do
143 m <- arbitrary 147 m <- arbitrary
144 let (u,_,v) = svd m 148 let (u,_,v) = svd m
@@ -157,7 +161,7 @@ instance (Convert a, Field a, Arbitrary a, Random (RealOf a)) => Arbitrary (WC a
157 161
158-- a well-conditioned square matrix (the singular values are between 1 and 100) 162-- a well-conditioned square matrix (the singular values are between 1 and 100)
159newtype (SqWC a) = SqWC (Matrix a) deriving Show 163newtype (SqWC a) = SqWC (Matrix a) deriving Show
160instance (Convert a, Field a, Arbitrary a, Random (RealOf a)) => Arbitrary (SqWC a) where 164instance (ArbitraryField a) => Arbitrary (SqWC a) where
161 arbitrary = do 165 arbitrary = do
162 Sq m <- arbitrary 166 Sq m <- arbitrary
163 let (u,_,v) = svd m 167 let (u,_,v) = svd m
@@ -174,7 +178,7 @@ instance (Convert a, Field a, Arbitrary a, Random (RealOf a)) => Arbitrary (SqWC
174 178
175-- a positive definite square matrix (the eigenvalues are between 0 and 100) 179-- a positive definite square matrix (the eigenvalues are between 0 and 100)
176newtype (PosDef a) = PosDef (Matrix a) deriving Show 180newtype (PosDef a) = PosDef (Matrix a) deriving Show
177instance (Convert a, Field a, Arbitrary a, Num (Vector a), Random (RealOf a)) 181instance (ArbitraryField a, Num (Vector a))
178 => Arbitrary (PosDef a) where 182 => Arbitrary (PosDef a) where
179 arbitrary = do 183 arbitrary = do
180 Her m <- arbitrary 184 Her m <- arbitrary