diff options
-rw-r--r-- | hmatrix.cabal | 2 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Matrix.hs | 16 | ||||
-rw-r--r-- | lib/Data/Packed/Matrix.hs | 1 | ||||
-rw-r--r-- | lib/Data/Packed/Vector.hs | 1 | ||||
-rw-r--r-- | lib/Numeric/Container.hs | 152 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Algorithms.hs | 2 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/LAPACK.hs | 1 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Linear.hs | 67 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests/Instances.hs | 14 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests/Properties.hs | 6 | ||||
-rw-r--r-- | lib/Numeric/Matrix.hs | 52 | ||||
-rw-r--r-- | lib/Numeric/Vector.hs | 98 |
12 files changed, 276 insertions, 136 deletions
diff --git a/hmatrix.cabal b/hmatrix.cabal index 736c443..de7b1cb 100644 --- a/hmatrix.cabal +++ b/hmatrix.cabal | |||
@@ -123,7 +123,7 @@ library | |||
123 | cpp-options: -DVECTOR | 123 | cpp-options: -DVECTOR |
124 | 124 | ||
125 | if flag(tests) | 125 | if flag(tests) |
126 | Build-Depends: QuickCheck, HUnit | 126 | Build-Depends: QuickCheck, HUnit, random |
127 | exposed-modules: Numeric.LinearAlgebra.Tests | 127 | exposed-modules: Numeric.LinearAlgebra.Tests |
128 | other-modules: Numeric.LinearAlgebra.Tests.Instances, | 128 | other-modules: Numeric.LinearAlgebra.Tests.Instances, |
129 | Numeric.LinearAlgebra.Tests.Properties | 129 | Numeric.LinearAlgebra.Tests.Properties |
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs index 100b8ea..e66817e 100644 --- a/lib/Data/Packed/Internal/Matrix.hs +++ b/lib/Data/Packed/Internal/Matrix.hs | |||
@@ -29,7 +29,6 @@ module Data.Packed.Internal.Matrix( | |||
29 | liftMatrix, liftMatrix2, | 29 | liftMatrix, liftMatrix2, |
30 | (@@>), | 30 | (@@>), |
31 | saveMatrix, | 31 | saveMatrix, |
32 | fromComplexV, toComplexV, conjV, | ||
33 | singleton | 32 | singleton |
34 | ) where | 33 | ) where |
35 | 34 | ||
@@ -386,21 +385,6 @@ subMatrix' (r0,c0) (rt,ct) m = trans $ subMatrix' (c0,r0) (ct,rt) (trans m) | |||
386 | 385 | ||
387 | -------------------------------------------------------------------------- | 386 | -------------------------------------------------------------------------- |
388 | 387 | ||
389 | -- | obtains the complex conjugate of a complex vector | ||
390 | conjV :: (Storable a, RealFloat a) => Vector (Complex a) -> Vector (Complex a) | ||
391 | conjV = mapVector conjugate | ||
392 | |||
393 | -- | creates a complex vector from vectors with real and imaginary parts | ||
394 | toComplexV :: (RealFloat a, Element a) => (Vector a, Vector a) -> Vector (Complex a) | ||
395 | toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i] | ||
396 | |||
397 | -- | the inverse of 'toComplex' | ||
398 | fromComplexV :: (RealFloat a, Element a) => Vector (Complex a) -> (Vector a, Vector a) | ||
399 | fromComplexV z = (r,i) where | ||
400 | [r,i] = toColumns $ reshape 2 $ asReal z | ||
401 | |||
402 | -------------------------------------------------------------------------- | ||
403 | |||
404 | -- | Saves a matrix as 2D ASCII table. | 388 | -- | Saves a matrix as 2D ASCII table. |
405 | saveMatrix :: FilePath | 389 | saveMatrix :: FilePath |
406 | -> String -- ^ format (%f, %g, %e) | 390 | -> String -- ^ format (%f, %g, %e) |
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs index b860cb1..03e5889 100644 --- a/lib/Data/Packed/Matrix.hs +++ b/lib/Data/Packed/Matrix.hs | |||
@@ -55,6 +55,7 @@ import Control.Monad(replicateM) | |||
55 | --import Control.Arrow((***)) | 55 | --import Control.Arrow((***)) |
56 | --import GHC.Float(double2Float,float2Double) | 56 | --import GHC.Float(double2Float,float2Double) |
57 | 57 | ||
58 | |||
58 | ------------------------------------------------------------------- | 59 | ------------------------------------------------------------------- |
59 | 60 | ||
60 | instance (Binary a, Element a, Storable a) => Binary (Matrix a) where | 61 | instance (Binary a, Element a, Storable a) => Binary (Matrix a) where |
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs index ad690f9..253f2fa 100644 --- a/lib/Data/Packed/Vector.hs +++ b/lib/Data/Packed/Vector.hs | |||
@@ -72,7 +72,6 @@ instance (Binary a, Storable a) => Binary (Vector a) where | |||
72 | 72 | ||
73 | ------------------------------------------------------------------- | 73 | ------------------------------------------------------------------- |
74 | 74 | ||
75 | |||
76 | {- | 75 | {- |
77 | vectorFMax :: Vector Float -> Float | 76 | vectorFMax :: Vector Float -> Float |
78 | vectorFMax = toScalarF Max | 77 | vectorFMax = toScalarF Max |
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs index 010235f..a0e489b 100644 --- a/lib/Numeric/Container.hs +++ b/lib/Numeric/Container.hs | |||
@@ -19,11 +19,13 @@ | |||
19 | ----------------------------------------------------------------------------- | 19 | ----------------------------------------------------------------------------- |
20 | 20 | ||
21 | module Numeric.Container ( | 21 | module Numeric.Container ( |
22 | Container(..), RealElement, Precision, NumericContainer(..), comp, | 22 | Linear(..), |
23 | Convert(..), AutoReal(..), | 23 | Container(..), RealElement, Precision(..), NumericContainer(..), comp, |
24 | RealOf, ComplexOf, SingleOf, DoubleOf, | 24 | -- Complexable(..), Precisionable(..), |
25 | Convert(..), --AutoReal(..), | ||
26 | RealOf, ComplexOf, SingleOf, DoubleOf, | ||
25 | 27 | ||
26 | -- ElementOf, | 28 | -- ElementOf, |
27 | 29 | ||
28 | IndexOf, | 30 | IndexOf, |
29 | 31 | ||
@@ -33,11 +35,10 @@ module Numeric.Container ( | |||
33 | import Data.Packed.Vector | 35 | import Data.Packed.Vector |
34 | import Data.Packed.Matrix | 36 | import Data.Packed.Matrix |
35 | import Data.Packed.Internal.Vector | 37 | import Data.Packed.Internal.Vector |
36 | import Data.Packed.Internal.Matrix | 38 | --import Data.Packed.Internal.Matrix |
37 | --import qualified Data.Packed.ST as ST | 39 | --import qualified Data.Packed.ST as ST |
38 | 40 | ||
39 | import Control.Arrow((***)) | 41 | --import Control.Arrow((***)) |
40 | |||
41 | import Data.Complex | 42 | import Data.Complex |
42 | 43 | ||
43 | ------------------------------------------------------------------- | 44 | ------------------------------------------------------------------- |
@@ -71,7 +72,7 @@ class NumericContainer c where | |||
71 | fromComplex :: (RealElement e) => c (Complex e) -> (c e, c e) | 72 | fromComplex :: (RealElement e) => c (Complex e) -> (c e, c e) |
72 | complex' :: (RealElement e) => c e -> c (Complex e) | 73 | complex' :: (RealElement e) => c e -> c (Complex e) |
73 | conj :: (RealElement e) => c (Complex e) -> c (Complex e) | 74 | conj :: (RealElement e) => c (Complex e) -> c (Complex e) |
74 | cmap :: (Element a, Element b) => (a -> b) -> c a -> c b | 75 | -- cmap :: (Element a, Element b) => (a -> b) -> c a -> c b |
75 | single' :: Precision a b => c b -> c a | 76 | single' :: Precision a b => c b -> c a |
76 | double' :: Precision a b => c a -> c b | 77 | double' :: Precision a b => c a -> c b |
77 | 78 | ||
@@ -79,25 +80,6 @@ class NumericContainer c where | |||
79 | comp :: (NumericContainer c, RealElement e) => c e -> c (Complex e) | 80 | comp :: (NumericContainer c, RealElement e) => c e -> c (Complex e) |
80 | comp x = complex' x | 81 | comp x = complex' x |
81 | 82 | ||
82 | instance NumericContainer Vector where | ||
83 | toComplex = toComplexV | ||
84 | fromComplex = fromComplexV | ||
85 | complex' v = toComplex (v,constantD 0 (dim v)) | ||
86 | conj = conjV | ||
87 | cmap = mapVector | ||
88 | single' = double2FloatG | ||
89 | double' = float2DoubleG | ||
90 | |||
91 | instance NumericContainer Matrix where | ||
92 | toComplex = uncurry $ liftMatrix2 $ curry toComplex | ||
93 | fromComplex z = (reshape c *** reshape c) . fromComplex . flatten $ z | ||
94 | where c = cols z | ||
95 | complex' = liftMatrix complex' | ||
96 | conj = liftMatrix conj | ||
97 | cmap f = liftMatrix (cmap f) | ||
98 | single' = liftMatrix single' | ||
99 | double' = liftMatrix double' | ||
100 | |||
101 | ------------------------------------------------------------------- | 83 | ------------------------------------------------------------------- |
102 | 84 | ||
103 | type family RealOf x | 85 | type family RealOf x |
@@ -141,14 +123,78 @@ type instance IndexOf Vector = Int | |||
141 | type instance IndexOf Matrix = (Int,Int) | 123 | type instance IndexOf Matrix = (Int,Int) |
142 | 124 | ||
143 | ------------------------------------------------------------------- | 125 | ------------------------------------------------------------------- |
126 | {- | ||
127 | -- | Supported single-double precision type pairs | ||
128 | class (Element e) => V_Precision e where | ||
129 | v_double2FloatG :: Vector e -> Vector (SingleOf e) | ||
130 | v_float2DoubleG :: Vector (SingleOf e) -> Vector e | ||
131 | {- | ||
132 | instance V_Precision Float where | ||
133 | v_double2FloatG = double2FloatV | ||
134 | v_float2DoubleG = float2DoubleV | ||
135 | -} | ||
136 | instance V_Precision Double where | ||
137 | v_double2FloatG = double2FloatV | ||
138 | v_float2DoubleG = float2DoubleV | ||
139 | {- | ||
140 | instance V_Precision (Complex Float) where | ||
141 | v_double2FloatG = asComplex . double2FloatV . asReal | ||
142 | v_float2DoubleG = asComplex . float2DoubleV . asReal | ||
143 | -} | ||
144 | instance V_Precision (Complex Double) where | ||
145 | v_double2FloatG = asComplex . double2FloatV . asReal | ||
146 | v_float2DoubleG = asComplex . float2DoubleV . asReal | ||
147 | -} | ||
148 | ------------------------------------------------------------------- | ||
149 | {- | ||
150 | -- | converting to/from complex containers | ||
151 | class RealElement t => Complexable c t where | ||
152 | v_toComplex :: (c t, c t) -> c (Complex t) | ||
153 | v_fromComplex :: c (Complex t) -> (c t, c t) | ||
154 | v_complex' :: c t -> c (Complex t) | ||
155 | v_conj :: c (Complex t) -> c (Complex t) | ||
156 | |||
157 | -- | converting to/from single/double precision numbers | ||
158 | class (Element (SingleOf t), Element t, RealElement (RealOf t)) => Precisionable c t where | ||
159 | v_single' :: (V_Precision (DoubleOf t)) => c t -> c (SingleOf t) | ||
160 | v_double' :: (V_Precision (DoubleOf t)) => c t -> c (DoubleOf t) | ||
144 | 161 | ||
145 | -- | generic conversion functions | 162 | -- | generic conversion functions |
146 | class Convert t where | 163 | class (Element t, RealElement (RealOf t)) => V_Convert t where |
164 | -- | real/complex | ||
165 | v_real :: Complexable c (RealOf t) => c (RealOf t) -> c t -- from the instances, this looks like it turns a real object into a complex object WHEN the context is a complex object | ||
166 | v_complex :: Complexable c (RealOf t) => c t -> c (ComplexOf t) | ||
167 | -- | single/double | ||
168 | v_single :: Precisionable c t => c t -> c (SingleOf t) | ||
169 | v_double :: Precisionable c t => c t -> c (DoubleOf t) | ||
170 | -} | ||
171 | ------------------------------------------------------------------- | ||
172 | {- | ||
173 | instance Precisionable Vector Float where | ||
174 | v_single' = id | ||
175 | v_double' = float2DoubleG | ||
176 | |||
177 | instance Precisionable Vector Double where | ||
178 | v_single' = double2FloatG | ||
179 | v_double' = id | ||
180 | |||
181 | instance Precisionable Vector (Complex Float) where | ||
182 | v_single' = id | ||
183 | v_double' = float2DoubleG | ||
184 | |||
185 | instance Precisionable Vector (Complex Double) where | ||
186 | v_single' = double2FloatG | ||
187 | v_double' = id | ||
188 | -} | ||
189 | ------------------------------------------------------------------- | ||
190 | |||
191 | class (Element t, Element (RealOf t)) => Convert t where | ||
147 | real :: NumericContainer c => c (RealOf t) -> c t | 192 | real :: NumericContainer c => c (RealOf t) -> c t |
148 | complex :: NumericContainer c => c t -> c (ComplexOf t) | 193 | complex :: NumericContainer c => c t -> c (ComplexOf t) |
149 | single :: NumericContainer c => c t -> c (SingleOf t) | 194 | single :: NumericContainer c => c t -> c (SingleOf t) |
150 | double :: NumericContainer c => c t -> c (DoubleOf t) | 195 | double :: NumericContainer c => c t -> c (DoubleOf t) |
151 | 196 | ||
197 | |||
152 | instance Convert Double where | 198 | instance Convert Double where |
153 | real = id | 199 | real = id |
154 | complex = complex' | 200 | complex = complex' |
@@ -180,6 +226,7 @@ class Convert t => AutoReal t where | |||
180 | real'' :: NumericContainer c => c Double -> c t | 226 | real'' :: NumericContainer c => c Double -> c t |
181 | complex'' :: NumericContainer c => c t -> c (Complex Double) | 227 | complex'' :: NumericContainer c => c t -> c (Complex Double) |
182 | 228 | ||
229 | |||
183 | instance AutoReal Double where | 230 | instance AutoReal Double where |
184 | real'' = real | 231 | real'' = real |
185 | complex'' = complex | 232 | complex'' = complex |
@@ -198,13 +245,60 @@ instance AutoReal (Complex Float) where | |||
198 | 245 | ||
199 | ------------------------------------------------------------------- | 246 | ------------------------------------------------------------------- |
200 | 247 | ||
201 | -- | Basic element-by-element functions. | 248 | -- | Basic element-by-element functions for numeric containers |
202 | class (Element e) => Container c e where | 249 | class (Element e) => Container c e where |
250 | {- | ||
251 | -- | create a structure with a single element | ||
252 | scalar :: e -> c e | ||
253 | -- | multiply every element by a scalar | ||
254 | scale :: e -> c e -> c e | ||
255 | -- | scale the element by element reciprocal of the object: | ||
256 | -- | ||
257 | -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@ | ||
258 | scaleRecip :: e -> c e -> c e | ||
259 | -- | add a constant to each element | ||
260 | addConstant :: e -> c e -> c e | ||
261 | add :: c e -> c e -> c e | ||
262 | sub :: c e -> c e -> c e | ||
263 | -- | element by element multiplication | ||
264 | mul :: c e -> c e -> c e | ||
265 | -- | element by element division | ||
266 | divide :: c e -> c e -> c e | ||
267 | equal :: c e -> c e -> Bool | ||
268 | -} | ||
269 | -- | cannot implement instance Functor because of Element class constraint | ||
270 | cmap :: (Element a, Element b) => (a -> b) -> c a -> c b | ||
271 | -- | ||
272 | -- | indexing function | ||
273 | atIndex :: c e -> IndexOf c -> e | ||
274 | -- | index of min/max element | ||
203 | minIndex :: c e -> IndexOf c | 275 | minIndex :: c e -> IndexOf c |
204 | maxIndex :: c e -> IndexOf c | 276 | maxIndex :: c e -> IndexOf c |
277 | -- | value of min/max element | ||
205 | minElement :: c e -> e | 278 | minElement :: c e -> e |
206 | maxElement :: c e -> e | 279 | maxElement :: c e -> e |
280 | -- the C functions sumX/prodX are twice as fast as using foldVector | ||
281 | -- | the sum/product of elements (faster than using @fold@ | ||
282 | sumElements :: c e -> e | ||
283 | prodElements :: c e -> e | ||
207 | 284 | ||
285 | -- | Basic element-by-element functions. | ||
286 | class (Element e, Container c e) => Linear c e where | ||
287 | -- | create a structure with a single element | ||
288 | scalar :: e -> c e | ||
289 | scale :: e -> c e -> c e | ||
290 | -- | scale the element by element reciprocal of the object: | ||
291 | -- | ||
292 | -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@ | ||
293 | scaleRecip :: e -> c e -> c e | ||
294 | addConstant :: e -> c e -> c e | ||
295 | add :: c e -> c e -> c e | ||
296 | sub :: c e -> c e -> c e | ||
297 | -- | element by element multiplication | ||
298 | mul :: c e -> c e -> c e | ||
299 | -- | element by element division | ||
300 | divide :: c e -> c e -> c e | ||
301 | equal :: c e -> c e -> Bool | ||
208 | 302 | ||
209 | 303 | ||
210 | 304 | ||
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index ac46847..8306961 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs | |||
@@ -85,7 +85,7 @@ import Numeric.Vector | |||
85 | import Numeric.Matrix() | 85 | import Numeric.Matrix() |
86 | 86 | ||
87 | -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. | 87 | -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. |
88 | class (Product t, Linear Vector t, Linear Matrix t) => Field t where | 88 | class (Product t, Linear Vector t, Container Vector t, Container Matrix t) => Field t where |
89 | svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) | 89 | svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) |
90 | thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) | 90 | thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) |
91 | sv' :: Matrix t -> Vector Double | 91 | sv' :: Matrix t -> Vector Double |
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs index 8888712..5d0154d 100644 --- a/lib/Numeric/LinearAlgebra/LAPACK.hs +++ b/lib/Numeric/LinearAlgebra/LAPACK.hs | |||
@@ -44,6 +44,7 @@ module Numeric.LinearAlgebra.LAPACK ( | |||
44 | import Data.Packed.Internal | 44 | import Data.Packed.Internal |
45 | import Data.Packed.Matrix | 45 | import Data.Packed.Matrix |
46 | import Data.Complex | 46 | import Data.Complex |
47 | import Numeric.Vector() | ||
47 | import Numeric.Container | 48 | import Numeric.Container |
48 | import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) | 49 | import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) |
49 | import Foreign | 50 | import Foreign |
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs index 952661d..775060e 100644 --- a/lib/Numeric/LinearAlgebra/Linear.hs +++ b/lib/Numeric/LinearAlgebra/Linear.hs | |||
@@ -18,7 +18,7 @@ Basic optimized operations on vectors and matrices. | |||
18 | 18 | ||
19 | module Numeric.LinearAlgebra.Linear ( | 19 | module Numeric.LinearAlgebra.Linear ( |
20 | -- * Linear Algebra Typeclasses | 20 | -- * Linear Algebra Typeclasses |
21 | Vectors(..), Linear(..), | 21 | Vectors(..), |
22 | -- * Products | 22 | -- * Products |
23 | Product(..), | 23 | Product(..), |
24 | mXm,mXv,vXm, | 24 | mXm,mXv,vXm, |
@@ -34,22 +34,52 @@ import Data.Packed.Matrix | |||
34 | import Data.Packed.Vector | 34 | import Data.Packed.Vector |
35 | import Data.Complex | 35 | import Data.Complex |
36 | import Numeric.Container | 36 | import Numeric.Container |
37 | --import Numeric.Vector | 37 | import Numeric.Vector() |
38 | --import Numeric.Matrix | 38 | import Numeric.Matrix() |
39 | --import Numeric.GSL.Vector | 39 | import Numeric.GSL.Vector |
40 | import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) | 40 | import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) |
41 | 41 | ||
42 | -- | basic Vector functions | 42 | -- | Linear algebraic properties of objects |
43 | class Num e => Vectors a e where | 43 | class Num e => Vectors a e where |
44 | -- the C functions sumX are twice as fast as using foldVector | 44 | -- | dot (inner) product |
45 | vectorSum :: a e -> e | ||
46 | vectorProd :: a e -> e | ||
47 | absSum :: a e -> e | ||
48 | dot :: a e -> a e -> e | 45 | dot :: a e -> a e -> e |
46 | -- | sum of absolute value of elements (differs in complex case from @norm1@ | ||
47 | absSum :: a e -> e | ||
48 | -- | sum of absolute value of elements | ||
49 | norm1 :: a e -> e | 49 | norm1 :: a e -> e |
50 | -- | euclidean norm | ||
50 | norm2 :: a e -> e | 51 | norm2 :: a e -> e |
52 | -- | element of maximum magnitude | ||
51 | normInf :: a e -> e | 53 | normInf :: a e -> e |
52 | 54 | ||
55 | instance Vectors Vector Float where | ||
56 | norm2 = toScalarF Norm2 | ||
57 | absSum = toScalarF AbsSum | ||
58 | dot = dotF | ||
59 | norm1 = toScalarF AbsSum | ||
60 | normInf = maxElement . vectorMapF Abs | ||
61 | |||
62 | instance Vectors Vector Double where | ||
63 | norm2 = toScalarR Norm2 | ||
64 | absSum = toScalarR AbsSum | ||
65 | dot = dotR | ||
66 | norm1 = toScalarR AbsSum | ||
67 | normInf = maxElement . vectorMapR Abs | ||
68 | |||
69 | instance Vectors Vector (Complex Float) where | ||
70 | norm2 = (:+ 0) . toScalarQ Norm2 | ||
71 | absSum = (:+ 0) . toScalarQ AbsSum | ||
72 | dot = dotQ | ||
73 | norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapQ Abs | ||
74 | normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapQ Abs | ||
75 | |||
76 | instance Vectors Vector (Complex Double) where | ||
77 | norm2 = (:+ 0) . toScalarC Norm2 | ||
78 | absSum = (:+ 0) . toScalarC AbsSum | ||
79 | dot = dotC | ||
80 | norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapC Abs | ||
81 | normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapC Abs | ||
82 | |||
53 | ---------------------------------------------------- | 83 | ---------------------------------------------------- |
54 | 84 | ||
55 | class Element t => Product t where | 85 | class Element t => Product t where |
@@ -128,22 +158,3 @@ kronecker a b = fromBlocks | |||
128 | 158 | ||
129 | 159 | ||
130 | ------------------------------------------------------------------- | 160 | ------------------------------------------------------------------- |
131 | |||
132 | |||
133 | -- | Basic element-by-element functions. | ||
134 | class (Element e, Container c e) => Linear c e where | ||
135 | -- | create a structure with a single element | ||
136 | scalar :: e -> c e | ||
137 | scale :: e -> c e -> c e | ||
138 | -- | scale the element by element reciprocal of the object: | ||
139 | -- | ||
140 | -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@ | ||
141 | scaleRecip :: e -> c e -> c e | ||
142 | addConstant :: e -> c e -> c e | ||
143 | add :: c e -> c e -> c e | ||
144 | sub :: c e -> c e -> c e | ||
145 | -- | element by element multiplication | ||
146 | mul :: c e -> c e -> c e | ||
147 | -- | element by element division | ||
148 | divide :: c e -> c e -> c e | ||
149 | equal :: c e -> c e -> Bool | ||
diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs index 21a6f88..6046ccb 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs | |||
@@ -26,6 +26,7 @@ module Numeric.LinearAlgebra.Tests.Instances( | |||
26 | FM,ZM, fM,zM | 26 | FM,ZM, fM,zM |
27 | ) where | 27 | ) where |
28 | 28 | ||
29 | import System.Random | ||
29 | 30 | ||
30 | import Numeric.LinearAlgebra | 31 | import Numeric.LinearAlgebra |
31 | import Control.Monad(replicateM) | 32 | import Control.Monad(replicateM) |
@@ -137,7 +138,7 @@ instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where | |||
137 | 138 | ||
138 | -- a well-conditioned general matrix (the singular values are between 1 and 100) | 139 | -- a well-conditioned general matrix (the singular values are between 1 and 100) |
139 | newtype (WC a) = WC (Matrix a) deriving Show | 140 | newtype (WC a) = WC (Matrix a) deriving Show |
140 | instance (AutoReal a, Field a, Arbitrary a) => Arbitrary (WC a) where | 141 | instance (Convert a, Field a, Arbitrary a, Random (RealOf a)) => Arbitrary (WC a) where |
141 | arbitrary = do | 142 | arbitrary = do |
142 | m <- arbitrary | 143 | m <- arbitrary |
143 | let (u,_,v) = svd m | 144 | let (u,_,v) = svd m |
@@ -146,7 +147,7 @@ instance (AutoReal a, Field a, Arbitrary a) => Arbitrary (WC a) where | |||
146 | n = min r c | 147 | n = min r c |
147 | sv' <- replicateM n (choose (1,100)) | 148 | sv' <- replicateM n (choose (1,100)) |
148 | let s = diagRect (fromList sv') r c | 149 | let s = diagRect (fromList sv') r c |
149 | return $ WC (u <> real'' s <> trans v) | 150 | return $ WC (u <> real s <> trans v) |
150 | 151 | ||
151 | #if MIN_VERSION_QuickCheck(2,0,0) | 152 | #if MIN_VERSION_QuickCheck(2,0,0) |
152 | #else | 153 | #else |
@@ -156,14 +157,14 @@ instance (AutoReal a, Field a, Arbitrary a) => Arbitrary (WC a) where | |||
156 | 157 | ||
157 | -- a well-conditioned square matrix (the singular values are between 1 and 100) | 158 | -- a well-conditioned square matrix (the singular values are between 1 and 100) |
158 | newtype (SqWC a) = SqWC (Matrix a) deriving Show | 159 | newtype (SqWC a) = SqWC (Matrix a) deriving Show |
159 | instance (AutoReal a, Field a, Arbitrary a) => Arbitrary (SqWC a) where | 160 | instance (Convert a, Field a, Arbitrary a, Random (RealOf a)) => Arbitrary (SqWC a) where |
160 | arbitrary = do | 161 | arbitrary = do |
161 | Sq m <- arbitrary | 162 | Sq m <- arbitrary |
162 | let (u,_,v) = svd m | 163 | let (u,_,v) = svd m |
163 | n = rows m | 164 | n = rows m |
164 | sv' <- replicateM n (choose (1,100)) | 165 | sv' <- replicateM n (choose (1,100)) |
165 | let s = diag (fromList sv') | 166 | let s = diag (fromList sv') |
166 | return $ SqWC (u <> real'' s <> trans v) | 167 | return $ SqWC (u <> real s <> trans v) |
167 | 168 | ||
168 | #if MIN_VERSION_QuickCheck(2,0,0) | 169 | #if MIN_VERSION_QuickCheck(2,0,0) |
169 | #else | 170 | #else |
@@ -173,14 +174,15 @@ instance (AutoReal a, Field a, Arbitrary a) => Arbitrary (SqWC a) where | |||
173 | 174 | ||
174 | -- a positive definite square matrix (the eigenvalues are between 0 and 100) | 175 | -- a positive definite square matrix (the eigenvalues are between 0 and 100) |
175 | newtype (PosDef a) = PosDef (Matrix a) deriving Show | 176 | newtype (PosDef a) = PosDef (Matrix a) deriving Show |
176 | instance (AutoReal a, Field a, Arbitrary a, Num (Vector a)) => Arbitrary (PosDef a) where | 177 | instance (Convert a, Field a, Arbitrary a, Num (Vector a), Random (RealOf a)) |
178 | => Arbitrary (PosDef a) where | ||
177 | arbitrary = do | 179 | arbitrary = do |
178 | Her m <- arbitrary | 180 | Her m <- arbitrary |
179 | let (_,v) = eigSH m | 181 | let (_,v) = eigSH m |
180 | n = rows m | 182 | n = rows m |
181 | l <- replicateM n (choose (0,100)) | 183 | l <- replicateM n (choose (0,100)) |
182 | let s = diag (fromList l) | 184 | let s = diag (fromList l) |
183 | p = v <> real'' s <> ctrans v | 185 | p = v <> real s <> ctrans v |
184 | return $ PosDef (0.5 * p + 0.5 * ctrans p) | 186 | return $ PosDef (0.5 * p + 0.5 * ctrans p) |
185 | 187 | ||
186 | #if MIN_VERSION_QuickCheck(2,0,0) | 188 | #if MIN_VERSION_QuickCheck(2,0,0) |
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs index d312e52..b96f53e 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs | |||
@@ -42,14 +42,14 @@ module Numeric.LinearAlgebra.Tests.Properties ( | |||
42 | linearSolveProp, linearSolveProp2 | 42 | linearSolveProp, linearSolveProp2 |
43 | ) where | 43 | ) where |
44 | 44 | ||
45 | import Numeric.LinearAlgebra hiding (real,complex) | 45 | import Numeric.LinearAlgebra --hiding (real,complex) |
46 | import Numeric.LinearAlgebra.LAPACK | 46 | import Numeric.LinearAlgebra.LAPACK |
47 | import Debug.Trace | 47 | import Debug.Trace |
48 | #include "quickCheckCompat.h" | 48 | #include "quickCheckCompat.h" |
49 | 49 | ||
50 | 50 | ||
51 | real x = real'' x | 51 | --real x = real'' x |
52 | complex x = complex'' x | 52 | --complex x = complex'' x |
53 | 53 | ||
54 | debug x = trace (show x) x | 54 | debug x = trace (show x) x |
55 | 55 | ||
diff --git a/lib/Numeric/Matrix.hs b/lib/Numeric/Matrix.hs index f240384..73515c1 100644 --- a/lib/Numeric/Matrix.hs +++ b/lib/Numeric/Matrix.hs | |||
@@ -27,11 +27,13 @@ module Numeric.Matrix ( | |||
27 | import Data.Packed.Vector | 27 | import Data.Packed.Vector |
28 | import Data.Packed.Matrix | 28 | import Data.Packed.Matrix |
29 | import Numeric.Container | 29 | import Numeric.Container |
30 | import Numeric.LinearAlgebra.Linear | 30 | --import Numeric.LinearAlgebra.Linear |
31 | --import Numeric.Vector | 31 | import Numeric.Vector() |
32 | 32 | ||
33 | import Control.Monad(ap) | 33 | import Control.Monad(ap) |
34 | 34 | ||
35 | import Control.Arrow((***)) | ||
36 | |||
35 | ------------------------------------------------------------------- | 37 | ------------------------------------------------------------------- |
36 | 38 | ||
37 | instance Linear Matrix a => Eq (Matrix a) where | 39 | instance Linear Matrix a => Eq (Matrix a) where |
@@ -75,7 +77,34 @@ instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floati | |||
75 | 77 | ||
76 | --------------------------------------------------------------- | 78 | --------------------------------------------------------------- |
77 | 79 | ||
78 | instance (Linear Vector a, NumericContainer Matrix) => (Linear Matrix a) where | 80 | instance NumericContainer Matrix where |
81 | toComplex = uncurry $ liftMatrix2 $ curry toComplex | ||
82 | fromComplex z = (reshape c *** reshape c) . fromComplex . flatten $ z | ||
83 | where c = cols z | ||
84 | complex' = liftMatrix complex' | ||
85 | conj = liftMatrix conj | ||
86 | -- cmap f = liftMatrix (cmap f) | ||
87 | single' = liftMatrix single' | ||
88 | double' = liftMatrix double' | ||
89 | |||
90 | --------------------------------------------------------------- | ||
91 | {- | ||
92 | instance (RealElement e, Complexable Vector e) => Complexable Matrix e where | ||
93 | v_toComplex = uncurry $ liftMatrix2 $ curry toComplex | ||
94 | v_fromComplex z = (reshape c *** reshape c) . fromComplex . flatten $ z | ||
95 | where c = cols z | ||
96 | v_conj = liftMatrix conj | ||
97 | v_complex' = liftMatrix complex' | ||
98 | |||
99 | --------------------------------------------------------------- | ||
100 | |||
101 | instance (Precisionable Vector e) => Precisionable Matrix e where | ||
102 | v_single' = liftMatrix single' | ||
103 | v_double' = liftMatrix double' | ||
104 | -} | ||
105 | --------------------------------------------------------------- | ||
106 | |||
107 | instance (Linear Vector a, Container Matrix a) => Linear Matrix a where | ||
79 | scale x = liftMatrix (scale x) | 108 | scale x = liftMatrix (scale x) |
80 | scaleRecip x = liftMatrix (scaleRecip x) | 109 | scaleRecip x = liftMatrix (scaleRecip x) |
81 | addConstant x = liftMatrix (addConstant x) | 110 | addConstant x = liftMatrix (addConstant x) |
@@ -85,16 +114,19 @@ instance (Linear Vector a, NumericContainer Matrix) => (Linear Matrix a) where | |||
85 | divide = liftMatrix2 divide | 114 | divide = liftMatrix2 divide |
86 | equal a b = cols a == cols b && flatten a `equal` flatten b | 115 | equal a b = cols a == cols b && flatten a `equal` flatten b |
87 | scalar x = (1><1) [x] | 116 | scalar x = (1><1) [x] |
88 | 117 | -- | |
89 | 118 | instance (Container Vector a) => Container Matrix a where | |
90 | instance (Linear Vector a, NumericContainer Matrix) => (Container Matrix a) where | 119 | cmap f = liftMatrix (mapVector f) |
120 | atIndex = (@@>) | ||
91 | minIndex m = let (r,c) = (rows m,cols m) | 121 | minIndex m = let (r,c) = (rows m,cols m) |
92 | i = 1 + (minIndex $ flatten m) | 122 | i = (minIndex $ flatten m) |
93 | in (i `div` r,i `mod` r) | 123 | in (i `div` c,(i `mod` c) + 1) |
94 | maxIndex m = let (r,c) = (rows m,cols m) | 124 | maxIndex m = let (r,c) = (rows m,cols m) |
95 | i = 1 + (maxIndex $ flatten m) | 125 | i = (maxIndex $ flatten m) |
96 | in (i `div` r,i `mod` r) | 126 | in (i `div` c,(i `mod` c) + 1) |
97 | minElement = ap (@@>) minIndex | 127 | minElement = ap (@@>) minIndex |
98 | maxElement = ap (@@>) maxIndex | 128 | maxElement = ap (@@>) maxIndex |
129 | sumElements = sumElements . flatten | ||
130 | prodElements = prodElements . flatten | ||
99 | 131 | ||
100 | ---------------------------------------------------- | 132 | ---------------------------------------------------- |
diff --git a/lib/Numeric/Vector.hs b/lib/Numeric/Vector.hs index d92a5e4..5cc51ac 100644 --- a/lib/Numeric/Vector.hs +++ b/lib/Numeric/Vector.hs | |||
@@ -30,10 +30,12 @@ import Control.Monad(ap) | |||
30 | 30 | ||
31 | import Data.Packed.Vector | 31 | import Data.Packed.Vector |
32 | import Data.Packed.Internal.Matrix(Element(..)) | 32 | import Data.Packed.Internal.Matrix(Element(..)) |
33 | import Data.Packed.Internal.Vector(asComplex,asReal) | ||
34 | import Data.Packed.Matrix(toColumns,fromColumns,flatten,reshape) | ||
33 | import Numeric.GSL.Vector | 35 | import Numeric.GSL.Vector |
34 | 36 | ||
35 | import Numeric.Container | 37 | import Numeric.Container |
36 | import Numeric.LinearAlgebra.Linear | 38 | --import Numeric.LinearAlgebra.Linear |
37 | 39 | ||
38 | ------------------------------------------------------------------- | 40 | ------------------------------------------------------------------- |
39 | 41 | ||
@@ -254,6 +256,40 @@ instance Floating (Vector (Complex Float)) where | |||
254 | 256 | ||
255 | --------------------------------------------------------------- | 257 | --------------------------------------------------------------- |
256 | 258 | ||
259 | -- | obtains the complex conjugate of a complex vector | ||
260 | conjV :: (RealElement a) => Vector (Complex a) -> Vector (Complex a) | ||
261 | conjV = mapVector conjugate | ||
262 | |||
263 | -- | creates a complex vector from vectors with real and imaginary parts | ||
264 | toComplexV :: (RealElement a) => (Vector a, Vector a) -> Vector (Complex a) | ||
265 | toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i] | ||
266 | |||
267 | -- | the inverse of 'toComplex' | ||
268 | fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a) | ||
269 | fromComplexV z = (r,i) where | ||
270 | [r,i] = toColumns $ reshape 2 $ asReal z | ||
271 | |||
272 | -------------------------------------------------------------------------- | ||
273 | |||
274 | instance NumericContainer Vector where | ||
275 | toComplex = toComplexV | ||
276 | fromComplex = fromComplexV | ||
277 | complex' v = toComplex (v,constant 0 (dim v)) | ||
278 | conj = conjV | ||
279 | -- cmap = mapVector | ||
280 | single' = double2FloatG | ||
281 | double' = float2DoubleG | ||
282 | |||
283 | -------------------------------------------------------------------------- | ||
284 | {- | ||
285 | instance RealElement e => Complexable Vector e where | ||
286 | v_toComplex = toComplexV | ||
287 | v_fromComplex = fromComplexV | ||
288 | v_conj = conjV | ||
289 | v_complex' v = toComplex (v,constantD 0 (dim v)) | ||
290 | -} | ||
291 | ------------------------------------------------------------------- | ||
292 | |||
257 | instance Linear Vector Float where | 293 | instance Linear Vector Float where |
258 | scale = vectorMapValF Scale | 294 | scale = vectorMapValF Scale |
259 | scaleRecip = vectorMapValF Recip | 295 | scaleRecip = vectorMapValF Recip |
@@ -264,12 +300,16 @@ instance Linear Vector Float where | |||
264 | divide = vectorZipF Div | 300 | divide = vectorZipF Div |
265 | equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 | 301 | equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 |
266 | scalar x = fromList [x] | 302 | scalar x = fromList [x] |
267 | 303 | -- | |
268 | instance Container Vector Float where | 304 | instance Container Vector Float where |
305 | cmap = mapVector | ||
306 | atIndex = (@>) | ||
269 | minIndex = round . toScalarF MinIdx | 307 | minIndex = round . toScalarF MinIdx |
270 | maxIndex = round . toScalarF MaxIdx | 308 | maxIndex = round . toScalarF MaxIdx |
271 | minElement = toScalarF Min | 309 | minElement = toScalarF Min |
272 | maxElement = toScalarF Max | 310 | maxElement = toScalarF Max |
311 | sumElements = sumF | ||
312 | prodElements = prodF | ||
273 | 313 | ||
274 | instance Linear Vector Double where | 314 | instance Linear Vector Double where |
275 | scale = vectorMapValR Scale | 315 | scale = vectorMapValR Scale |
@@ -281,12 +321,16 @@ instance Linear Vector Double where | |||
281 | divide = vectorZipR Div | 321 | divide = vectorZipR Div |
282 | equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 | 322 | equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 |
283 | scalar x = fromList [x] | 323 | scalar x = fromList [x] |
284 | 324 | -- | |
285 | instance Container Vector Double where | 325 | instance Container Vector Double where |
326 | cmap = mapVector | ||
327 | atIndex = (@>) | ||
286 | minIndex = round . toScalarR MinIdx | 328 | minIndex = round . toScalarR MinIdx |
287 | maxIndex = round . toScalarR MaxIdx | 329 | maxIndex = round . toScalarR MaxIdx |
288 | minElement = toScalarR Min | 330 | minElement = toScalarR Min |
289 | maxElement = toScalarR Max | 331 | maxElement = toScalarR Max |
332 | sumElements = sumR | ||
333 | prodElements = prodR | ||
290 | 334 | ||
291 | instance Linear Vector (Complex Double) where | 335 | instance Linear Vector (Complex Double) where |
292 | scale = vectorMapValC Scale | 336 | scale = vectorMapValC Scale |
@@ -298,12 +342,16 @@ instance Linear Vector (Complex Double) where | |||
298 | divide = vectorZipC Div | 342 | divide = vectorZipC Div |
299 | equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 | 343 | equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 |
300 | scalar x = fromList [x] | 344 | scalar x = fromList [x] |
301 | 345 | -- | |
302 | instance Container Vector (Complex Double) where | 346 | instance Container Vector (Complex Double) where |
347 | cmap = mapVector | ||
348 | atIndex = (@>) | ||
303 | minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) | 349 | minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) |
304 | maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) | 350 | maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) |
305 | minElement = ap (@>) minIndex | 351 | minElement = ap (@>) minIndex |
306 | maxElement = ap (@>) maxIndex | 352 | maxElement = ap (@>) maxIndex |
353 | sumElements = sumC | ||
354 | prodElements = prodC | ||
307 | 355 | ||
308 | instance Linear Vector (Complex Float) where | 356 | instance Linear Vector (Complex Float) where |
309 | scale = vectorMapValQ Scale | 357 | scale = vectorMapValQ Scale |
@@ -315,47 +363,15 @@ instance Linear Vector (Complex Float) where | |||
315 | divide = vectorZipQ Div | 363 | divide = vectorZipQ Div |
316 | equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 | 364 | equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 |
317 | scalar x = fromList [x] | 365 | scalar x = fromList [x] |
318 | 366 | -- | |
319 | instance Container Vector (Complex Float) where | 367 | instance Container Vector (Complex Float) where |
368 | cmap = mapVector | ||
369 | atIndex = (@>) | ||
320 | minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) | 370 | minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) |
321 | maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) | 371 | maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) |
322 | minElement = ap (@>) minIndex | 372 | minElement = ap (@>) minIndex |
323 | maxElement = ap (@>) maxIndex | 373 | maxElement = ap (@>) maxIndex |
374 | sumElements = sumQ | ||
375 | prodElements = prodQ | ||
324 | 376 | ||
325 | --------------------------------------------------------------- | 377 | --------------------------------------------------------------- |
326 | |||
327 | instance Vectors Vector Float where | ||
328 | vectorSum = sumF | ||
329 | vectorProd = prodF | ||
330 | norm2 = toScalarF Norm2 | ||
331 | absSum = toScalarF AbsSum | ||
332 | dot = dotF | ||
333 | norm1 = toScalarF AbsSum | ||
334 | normInf = maxElement . vectorMapF Abs | ||
335 | |||
336 | instance Vectors Vector Double where | ||
337 | vectorSum = sumR | ||
338 | vectorProd = prodR | ||
339 | norm2 = toScalarR Norm2 | ||
340 | absSum = toScalarR AbsSum | ||
341 | dot = dotR | ||
342 | norm1 = toScalarR AbsSum | ||
343 | normInf = maxElement . vectorMapR Abs | ||
344 | |||
345 | instance Vectors Vector (Complex Float) where | ||
346 | vectorSum = sumQ | ||
347 | vectorProd = prodQ | ||
348 | norm2 = (:+ 0) . toScalarQ Norm2 | ||
349 | absSum = (:+ 0) . toScalarQ AbsSum | ||
350 | dot = dotQ | ||
351 | norm1 = (:+ 0) . vectorSum . fst . fromComplex . vectorMapQ Abs | ||
352 | normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapQ Abs | ||
353 | |||
354 | instance Vectors Vector (Complex Double) where | ||
355 | vectorSum = sumC | ||
356 | vectorProd = prodC | ||
357 | norm2 = (:+ 0) . toScalarC Norm2 | ||
358 | absSum = (:+ 0) . toScalarC AbsSum | ||
359 | dot = dotC | ||
360 | norm1 = (:+ 0) . vectorSum . fst . fromComplex . vectorMapC Abs | ||
361 | normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapC Abs \ No newline at end of file | ||