summaryrefslogtreecommitdiff
path: root/lib/Numeric
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric')
-rw-r--r--lib/Numeric/Container.hs418
-rw-r--r--lib/Numeric/LinearAlgebra.hs10
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs9
-rw-r--r--lib/Numeric/LinearAlgebra/Instances.hs276
-rw-r--r--lib/Numeric/LinearAlgebra/Interface.hs118
-rw-r--r--lib/Numeric/LinearAlgebra/LAPACK.hs3
-rw-r--r--lib/Numeric/LinearAlgebra/Linear.hs160
-rw-r--r--lib/Numeric/Matrix.hs157
-rw-r--r--lib/Numeric/Vector.hs141
9 files changed, 392 insertions, 900 deletions
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs
index d6b1342..fc4fb0d 100644
--- a/lib/Numeric/Container.hs
+++ b/lib/Numeric/Container.hs
@@ -3,6 +3,7 @@
3{-# LANGUAGE FlexibleInstances #-} 3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE MultiParamTypeClasses #-} 4{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE FunctionalDependencies #-} 5{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE UndecidableInstances #-}
6 7
7----------------------------------------------------------------------------- 8-----------------------------------------------------------------------------
8-- | 9-- |
@@ -19,103 +20,35 @@
19----------------------------------------------------------------------------- 20-----------------------------------------------------------------------------
20 21
21module Numeric.Container ( 22module Numeric.Container (
22 Linear(..), 23 --Linear(..),
23 Container(..), RealElement, Precision(..), NumericContainer(..), comp, 24 Container(..),
24 Convert(..), --AutoReal(..), 25 Vectors(..),
25 RealOf, ComplexOf, SingleOf, DoubleOf, 26 Product(..),
26 27 mXm,mXv,vXm,
27-- ElementOf, 28 outer, kronecker,
29
30 RealElement, --Precision,
31 ComplexContainer(toComplex,fromComplex,comp,conj),
32 Convert(..), --AutoReal(..),
33 RealOf, ComplexOf, SingleOf, DoubleOf,
28 34
29 IndexOf, 35 IndexOf,
30
31 module Data.Complex 36 module Data.Complex
32) where 37) where
33 38
34import Data.Packed.Vector 39import Data.Packed
35import Data.Packed.Matrix 40import Numeric.Conversion
36import Data.Packed.Internal.Vector 41import Data.Packed.Internal
37--import Data.Packed.Internal.Matrix 42import Numeric.GSL.Vector
38--import qualified Data.Packed.ST as ST
39 43
40--import Control.Arrow((***))
41import Data.Complex 44import Data.Complex
45import Control.Monad(ap)
46--import Control.Arrow((***))
42 47
43------------------------------------------------------------------- 48import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
44
45-- | Supported single-double precision type pairs
46class (Element s, Element d) => Precision s d | s -> d, d -> s where
47 double2FloatG :: Vector d -> Vector s
48 float2DoubleG :: Vector s -> Vector d
49
50instance Precision Float Double where
51 double2FloatG = double2FloatV
52 float2DoubleG = float2DoubleV
53
54instance Precision (Complex Float) (Complex Double) where
55 double2FloatG = asComplex . double2FloatV . asReal
56 float2DoubleG = asComplex . float2DoubleV . asReal
57
58-- | Supported real types
59class (Element t, Element (Complex t), RealFloat t
60-- , RealOf t ~ t, RealOf (Complex t) ~ t
61 )
62 => RealElement t
63
64instance RealElement Double
65
66instance RealElement Float
67
68-- | Conversion utilities
69class NumericContainer c where
70 toComplex :: (RealElement e) => (c e, c e) -> c (Complex e)
71 fromComplex :: (RealElement e) => c (Complex e) -> (c e, c e)
72 complex' :: (RealElement e) => c e -> c (Complex e)
73 conj :: (RealElement e) => c (Complex e) -> c (Complex e)
74-- cmap :: (Element a, Element b) => (a -> b) -> c a -> c b
75 single' :: Precision a b => c b -> c a
76 double' :: Precision a b => c a -> c b
77
78-- | a synonym for "complex'"
79comp :: (NumericContainer c, RealElement e) => c e -> c (Complex e)
80comp x = complex' x
81 49
82------------------------------------------------------------------- 50-------------------------------------------------------------------
83 51
84type family RealOf x
85
86type instance RealOf Double = Double
87type instance RealOf (Complex Double) = Double
88
89type instance RealOf Float = Float
90type instance RealOf (Complex Float) = Float
91
92type family ComplexOf x
93
94type instance ComplexOf Double = Complex Double
95type instance ComplexOf (Complex Double) = Complex Double
96
97type instance ComplexOf Float = Complex Float
98type instance ComplexOf (Complex Float) = Complex Float
99
100type family SingleOf x
101
102type instance SingleOf Double = Float
103type instance SingleOf Float = Float
104
105type instance SingleOf (Complex a) = Complex (SingleOf a)
106
107type family DoubleOf x
108
109type instance DoubleOf Double = Double
110type instance DoubleOf Float = Double
111
112type instance DoubleOf (Complex a) = Complex (DoubleOf a)
113
114type family ElementOf c
115
116type instance ElementOf (Vector a) = a
117type instance ElementOf (Matrix a) = a
118
119type family IndexOf c 52type family IndexOf c
120 53
121type instance IndexOf Vector = Int 54type instance IndexOf Vector = Int
@@ -123,75 +56,16 @@ type instance IndexOf Matrix = (Int,Int)
123 56
124------------------------------------------------------------------- 57-------------------------------------------------------------------
125 58
126class (Element t, Element (RealOf t)) => Convert t where
127 real :: NumericContainer c => c (RealOf t) -> c t
128 complex :: NumericContainer c => c t -> c (ComplexOf t)
129 single :: NumericContainer c => c t -> c (SingleOf t)
130 double :: NumericContainer c => c t -> c (DoubleOf t)
131
132
133instance Convert Double where
134 real = id
135 complex = complex'
136 single = single'
137 double = id
138
139instance Convert Float where
140 real = id
141 complex = complex'
142 single = id
143 double = double'
144
145instance Convert (Complex Double) where
146 real = complex'
147 complex = id
148 single = single'
149 double = id
150
151instance Convert (Complex Float) where
152 real = complex'
153 complex = id
154 single = id
155 double = double'
156
157-------------------------------------------------------------------
158
159-- | to be replaced by Convert
160class Convert t => AutoReal t where
161 real'' :: NumericContainer c => c Double -> c t
162 complex'' :: NumericContainer c => c t -> c (Complex Double)
163
164
165instance AutoReal Double where
166 real'' = real
167 complex'' = complex
168
169instance AutoReal (Complex Double) where
170 real'' = real
171 complex'' = complex
172
173instance AutoReal Float where
174 real'' = real . single
175 complex'' = double . complex
176
177instance AutoReal (Complex Float) where
178 real'' = real . single
179 complex'' = double . complex
180
181-------------------------------------------------------------------
182
183-- | Basic element-by-element functions for numeric containers 59-- | Basic element-by-element functions for numeric containers
184class (Element e) => Container c e where 60class (Element e) => Container c e where
185{- 61
186 -- | create a structure with a single element 62 -- | create a structure with a single element
187 scalar :: e -> c e 63 scalar :: e -> c e
188 -- | multiply every element by a scalar
189 scale :: e -> c e -> c e 64 scale :: e -> c e -> c e
190 -- | scale the element by element reciprocal of the object: 65 -- | scale the element by element reciprocal of the object:
191 -- 66 --
192 -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@ 67 -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
193 scaleRecip :: e -> c e -> c e 68 scaleRecip :: e -> c e -> c e
194 -- | add a constant to each element
195 addConstant :: e -> c e -> c e 69 addConstant :: e -> c e -> c e
196 add :: c e -> c e -> c e 70 add :: c e -> c e -> c e
197 sub :: c e -> c e -> c e 71 sub :: c e -> c e -> c e
@@ -200,7 +74,7 @@ class (Element e) => Container c e where
200 -- | element by element division 74 -- | element by element division
201 divide :: c e -> c e -> c e 75 divide :: c e -> c e -> c e
202 equal :: c e -> c e -> Bool 76 equal :: c e -> c e -> Bool
203-} 77
204 -- | cannot implement instance Functor because of Element class constraint 78 -- | cannot implement instance Functor because of Element class constraint
205 cmap :: (Element a, Element b) => (a -> b) -> c a -> c b 79 cmap :: (Element a, Element b) => (a -> b) -> c a -> c b
206 -- 80 --
@@ -217,23 +91,239 @@ class (Element e) => Container c e where
217 sumElements :: c e -> e 91 sumElements :: c e -> e
218 prodElements :: c e -> e 92 prodElements :: c e -> e
219 93
220-- | Basic element-by-element functions. 94-- -- | Basic element-by-element functions.
221class (Element e, Container c e) => Linear c e where 95-- class (Element e, Container c e) => Linear c e where
222 -- | create a structure with a single element
223 scalar :: e -> c e
224 scale :: e -> c e -> c e
225 -- | scale the element by element reciprocal of the object:
226 --
227 -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
228 scaleRecip :: e -> c e -> c e
229 addConstant :: e -> c e -> c e
230 add :: c e -> c e -> c e
231 sub :: c e -> c e -> c e
232 -- | element by element multiplication
233 mul :: c e -> c e -> c e
234 -- | element by element division
235 divide :: c e -> c e -> c e
236 equal :: c e -> c e -> Bool
237 96
238 97
98--------------------------------------------------------------------------
239 99
100instance Container Vector Float where
101 scale = vectorMapValF Scale
102 scaleRecip = vectorMapValF Recip
103 addConstant = vectorMapValF AddConstant
104 add = vectorZipF Add
105 sub = vectorZipF Sub
106 mul = vectorZipF Mul
107 divide = vectorZipF Div
108 equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
109 scalar x = fromList [x]
110 --
111--instance Container Vector Float where
112 cmap = mapVector
113 atIndex = (@>)
114 minIndex = round . toScalarF MinIdx
115 maxIndex = round . toScalarF MaxIdx
116 minElement = toScalarF Min
117 maxElement = toScalarF Max
118 sumElements = sumF
119 prodElements = prodF
120
121instance Container Vector Double where
122 scale = vectorMapValR Scale
123 scaleRecip = vectorMapValR Recip
124 addConstant = vectorMapValR AddConstant
125 add = vectorZipR Add
126 sub = vectorZipR Sub
127 mul = vectorZipR Mul
128 divide = vectorZipR Div
129 equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
130 scalar x = fromList [x]
131 --
132--instance Container Vector Double where
133 cmap = mapVector
134 atIndex = (@>)
135 minIndex = round . toScalarR MinIdx
136 maxIndex = round . toScalarR MaxIdx
137 minElement = toScalarR Min
138 maxElement = toScalarR Max
139 sumElements = sumR
140 prodElements = prodR
141
142instance Container Vector (Complex Double) where
143 scale = vectorMapValC Scale
144 scaleRecip = vectorMapValC Recip
145 addConstant = vectorMapValC AddConstant
146 add = vectorZipC Add
147 sub = vectorZipC Sub
148 mul = vectorZipC Mul
149 divide = vectorZipC Div
150 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
151 scalar x = fromList [x]
152 --
153--instance Container Vector (Complex Double) where
154 cmap = mapVector
155 atIndex = (@>)
156 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
157 maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
158 minElement = ap (@>) minIndex
159 maxElement = ap (@>) maxIndex
160 sumElements = sumC
161 prodElements = prodC
162
163instance Container Vector (Complex Float) where
164 scale = vectorMapValQ Scale
165 scaleRecip = vectorMapValQ Recip
166 addConstant = vectorMapValQ AddConstant
167 add = vectorZipQ Add
168 sub = vectorZipQ Sub
169 mul = vectorZipQ Mul
170 divide = vectorZipQ Div
171 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
172 scalar x = fromList [x]
173 --
174--instance Container Vector (Complex Float) where
175 cmap = mapVector
176 atIndex = (@>)
177 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
178 maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
179 minElement = ap (@>) minIndex
180 maxElement = ap (@>) maxIndex
181 sumElements = sumQ
182 prodElements = prodQ
183
184---------------------------------------------------------------
185
186instance (Container Vector a) => Container Matrix a where
187 scale x = liftMatrix (scale x)
188 scaleRecip x = liftMatrix (scaleRecip x)
189 addConstant x = liftMatrix (addConstant x)
190 add = liftMatrix2 add
191 sub = liftMatrix2 sub
192 mul = liftMatrix2 mul
193 divide = liftMatrix2 divide
194 equal a b = cols a == cols b && flatten a `equal` flatten b
195 scalar x = (1><1) [x]
196 --
197--instance (Container Vector a) => Container Matrix a where
198 cmap f = liftMatrix (mapVector f)
199 atIndex = (@@>)
200 minIndex m = let (r,c) = (rows m,cols m)
201 i = (minIndex $ flatten m)
202 in (i `div` c,(i `mod` c) + 1)
203 maxIndex m = let (r,c) = (rows m,cols m)
204 i = (maxIndex $ flatten m)
205 in (i `div` c,(i `mod` c) + 1)
206 minElement = ap (@@>) minIndex
207 maxElement = ap (@@>) maxIndex
208 sumElements = sumElements . flatten
209 prodElements = prodElements . flatten
210
211----------------------------------------------------
212
213
214-- | Linear algebraic properties of objects
215class Num e => Vectors a e where
216 -- | dot (inner) product
217 dot :: a e -> a e -> e
218 -- | sum of absolute value of elements (differs in complex case from @norm1@
219 absSum :: a e -> e
220 -- | sum of absolute value of elements
221 norm1 :: a e -> e
222 -- | euclidean norm
223 norm2 :: a e -> e
224 -- | element of maximum magnitude
225 normInf :: a e -> e
226
227instance Vectors Vector Float where
228 norm2 = toScalarF Norm2
229 absSum = toScalarF AbsSum
230 dot = dotF
231 norm1 = toScalarF AbsSum
232 normInf = maxElement . vectorMapF Abs
233
234instance Vectors Vector Double where
235 norm2 = toScalarR Norm2
236 absSum = toScalarR AbsSum
237 dot = dotR
238 norm1 = toScalarR AbsSum
239 normInf = maxElement . vectorMapR Abs
240
241instance Vectors Vector (Complex Float) where
242 norm2 = (:+ 0) . toScalarQ Norm2
243 absSum = (:+ 0) . toScalarQ AbsSum
244 dot = dotQ
245 norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapQ Abs
246 normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapQ Abs
247
248instance Vectors Vector (Complex Double) where
249 norm2 = (:+ 0) . toScalarC Norm2
250 absSum = (:+ 0) . toScalarC AbsSum
251 dot = dotC
252 norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapC Abs
253 normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapC Abs
254
255----------------------------------------------------
256
257class Element t => Product t where
258 multiply :: Matrix t -> Matrix t -> Matrix t
259 ctrans :: Matrix t -> Matrix t
260
261instance Product Double where
262 multiply = multiplyR
263 ctrans = trans
264
265instance Product (Complex Double) where
266 multiply = multiplyC
267 ctrans = conj . trans
268
269instance Product Float where
270 multiply = multiplyF
271 ctrans = trans
272
273instance Product (Complex Float) where
274 multiply = multiplyQ
275 ctrans = conj . trans
276
277----------------------------------------------------------
278
279-- synonym for matrix product
280mXm :: Product t => Matrix t -> Matrix t -> Matrix t
281mXm = multiply
282
283-- matrix - vector product
284mXv :: Product t => Matrix t -> Vector t -> Vector t
285mXv m v = flatten $ m `mXm` (asColumn v)
286
287-- vector - matrix product
288vXm :: Product t => Vector t -> Matrix t -> Vector t
289vXm v m = flatten $ (asRow v) `mXm` m
290
291{- | Outer product of two vectors.
292
293@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
294(3><3)
295 [ 5.0, 2.0, 3.0
296 , 10.0, 4.0, 6.0
297 , 15.0, 6.0, 9.0 ]@
298-}
299outer :: (Product t) => Vector t -> Vector t -> Matrix t
300outer u v = asColumn u `multiply` asRow v
301
302{- | Kronecker product of two matrices.
303
304@m1=(2><3)
305 [ 1.0, 2.0, 0.0
306 , 0.0, -1.0, 3.0 ]
307m2=(4><3)
308 [ 1.0, 2.0, 3.0
309 , 4.0, 5.0, 6.0
310 , 7.0, 8.0, 9.0
311 , 10.0, 11.0, 12.0 ]@
312
313@\> kronecker m1 m2
314(8><9)
315 [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0
316 , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0
317 , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0
318 , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0
319 , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0
320 , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0
321 , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0
322 , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
323-}
324kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
325kronecker a b = fromBlocks
326 . splitEvery (cols a)
327 . map (reshape (cols b))
328 . toRows
329 $ flatten a `outer` flatten b
diff --git a/lib/Numeric/LinearAlgebra.hs b/lib/Numeric/LinearAlgebra.hs
index da4a4b5..edfd87e 100644
--- a/lib/Numeric/LinearAlgebra.hs
+++ b/lib/Numeric/LinearAlgebra.hs
@@ -14,14 +14,12 @@ This module reexports all normally required functions for Linear Algebra applica
14----------------------------------------------------------------------------- 14-----------------------------------------------------------------------------
15module Numeric.LinearAlgebra ( 15module Numeric.LinearAlgebra (
16 module Numeric.LinearAlgebra.Algorithms, 16 module Numeric.LinearAlgebra.Algorithms,
17 module Numeric.LinearAlgebra.Interface, 17-- module Numeric.LinearAlgebra.Interface,
18 module Numeric.LinearAlgebra.Linear,
19 module Numeric.Matrix, 18 module Numeric.Matrix,
20 module Numeric.Vector 19 module Data.Packed.Random
21) where 20) where
22 21
23import Numeric.LinearAlgebra.Algorithms 22import Numeric.LinearAlgebra.Algorithms
24import Numeric.LinearAlgebra.Interface 23--import Numeric.LinearAlgebra.Interface
25import Numeric.LinearAlgebra.Linear
26import Numeric.Matrix 24import Numeric.Matrix
27import Numeric.Vector 25import Data.Packed.Random
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 8306961..fa6c475 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -77,15 +77,16 @@ module Numeric.LinearAlgebra.Algorithms (
77import Data.Packed.Internal hiding ((//)) 77import Data.Packed.Internal hiding ((//))
78import Data.Packed.Matrix 78import Data.Packed.Matrix
79import Data.Complex 79import Data.Complex
80import Numeric.LinearAlgebra.Linear 80--import Numeric.LinearAlgebra.Linear
81import Numeric.LinearAlgebra.LAPACK as LAPACK 81import Numeric.LinearAlgebra.LAPACK as LAPACK
82import Data.List(foldl1') 82import Data.List(foldl1')
83import Data.Array 83import Data.Array
84import Numeric.Vector 84import Numeric.Container
85import Numeric.Matrix() 85
86constant x = constantD x
86 87
87-- | Auxiliary typeclass used to define generic computations for both real and complex matrices. 88-- | Auxiliary typeclass used to define generic computations for both real and complex matrices.
88class (Product t, Linear Vector t, Container Vector t, Container Matrix t) => Field t where 89class (Product t, Container Vector t, Container Matrix t) => Field t where
89 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) 90 svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
90 thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) 91 thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t)
91 sv' :: Matrix t -> Vector Double 92 sv' :: Matrix t -> Vector Double
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
deleted file mode 100644
index 04a9d88..0000000
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ /dev/null
@@ -1,276 +0,0 @@
1{-# LANGUAGE UndecidableInstances, FlexibleInstances #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Numeric.LinearAlgebra.Instances
5Copyright : (c) Alberto Ruiz 2006
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : portable
11
12This module exports Show, Read, Eq, Num, Fractional, and Floating instances for Vector and Matrix.
13
14In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand.
15
16-}
17-----------------------------------------------------------------------------
18
19module Numeric.LinearAlgebra.Instances(
20) where
21
22import Numeric.LinearAlgebra.Linear
23import Numeric.GSL.Vector
24import Data.Packed.Matrix
25import Data.Complex
26import Data.List(transpose,intersperse)
27import Data.Packed.Internal.Vector
28
29#ifndef VECTOR
30import Foreign(Storable)
31#endif
32
33------------------------------------------------------------------
34
35instance (Show a, Element a) => (Show (Matrix a)) where
36 show m = (sizes++) . dsp . map (map show) . toLists $ m
37 where sizes = "("++show (rows m)++"><"++show (cols m)++")\n"
38
39dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
40 where
41 mt = transpose as
42 longs = map (maximum . map length) mt
43 mtp = zipWith (\a b -> map (pad a) b) longs mt
44 pad n str = replicate (n - length str) ' ' ++ str
45 unwords' = concat . intersperse ", "
46
47#ifndef VECTOR
48
49instance (Show a, Storable a) => (Show (Vector a)) where
50 show v = (show (dim v))++" |> " ++ show (toList v)
51
52#endif
53
54------------------------------------------------------------------
55
56instance (Element a, Read a) => Read (Matrix a) where
57 readsPrec _ s = [((rs><cs) . read $ listnums, rest)]
58 where (thing,rest) = breakAt ']' s
59 (dims,listnums) = breakAt ')' thing
60 cs = read . init . fst. breakAt ')' . snd . breakAt '<' $ dims
61 rs = read . snd . breakAt '(' .init . fst . breakAt '>' $ dims
62
63#ifdef VECTOR
64
65instance (Element a, Read a) => Read (Vector a) where
66 readsPrec _ s = [(fromList . read $ listnums, rest)]
67 where (thing,trest) = breakAt ']' s
68 (dims,listnums) = breakAt ' ' (dropWhile (==' ') thing)
69 rest = drop 31 trest
70#else
71
72instance (Element a, Read a) => Read (Vector a) where
73 readsPrec _ s = [((d |>) . read $ listnums, rest)]
74 where (thing,rest) = breakAt ']' s
75 (dims,listnums) = breakAt '>' thing
76 d = read . init . fst . breakAt '|' $ dims
77
78#endif
79
80breakAt c l = (a++[c],tail b) where
81 (a,b) = break (==c) l
82
83------------------------------------------------------------------
84
85adaptScalar f1 f2 f3 x y
86 | dim x == 1 = f1 (x@>0) y
87 | dim y == 1 = f3 x (y@>0)
88 | otherwise = f2 x y
89
90#ifndef VECTOR
91
92instance Linear Vector a => Eq (Vector a) where
93 (==) = equal
94
95#endif
96
97instance Num (Vector Float) where
98 (+) = adaptScalar addConstant add (flip addConstant)
99 negate = scale (-1)
100 (*) = adaptScalar scale mul (flip scale)
101 signum = vectorMapF Sign
102 abs = vectorMapF Abs
103 fromInteger = fromList . return . fromInteger
104
105instance Num (Vector Double) where
106 (+) = adaptScalar addConstant add (flip addConstant)
107 negate = scale (-1)
108 (*) = adaptScalar scale mul (flip scale)
109 signum = vectorMapR Sign
110 abs = vectorMapR Abs
111 fromInteger = fromList . return . fromInteger
112
113instance Num (Vector (Complex Double)) where
114 (+) = adaptScalar addConstant add (flip addConstant)
115 negate = scale (-1)
116 (*) = adaptScalar scale mul (flip scale)
117 signum = vectorMapC Sign
118 abs = vectorMapC Abs
119 fromInteger = fromList . return . fromInteger
120
121instance Num (Vector (Complex Float)) where
122 (+) = adaptScalar addConstant add (flip addConstant)
123 negate = scale (-1)
124 (*) = adaptScalar scale mul (flip scale)
125 signum = vectorMapQ Sign
126 abs = vectorMapQ Abs
127 fromInteger = fromList . return . fromInteger
128
129instance Linear Matrix a => Eq (Matrix a) where
130 (==) = equal
131
132instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where
133 (+) = liftMatrix2Auto (+)
134 (-) = liftMatrix2Auto (-)
135 negate = liftMatrix negate
136 (*) = liftMatrix2Auto (*)
137 signum = liftMatrix signum
138 abs = liftMatrix abs
139 fromInteger = (1><1) . return . fromInteger
140
141---------------------------------------------------
142
143instance (Linear Vector a, Num (Vector a)) => Fractional (Vector a) where
144 fromRational n = fromList [fromRational n]
145 (/) = adaptScalar f divide g where
146 r `f` v = scaleRecip r v
147 v `g` r = scale (recip r) v
148
149-------------------------------------------------------
150
151instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where
152 fromRational n = (1><1) [fromRational n]
153 (/) = liftMatrix2Auto (/)
154
155---------------------------------------------------------
156
157instance Floating (Vector Float) where
158 sin = vectorMapF Sin
159 cos = vectorMapF Cos
160 tan = vectorMapF Tan
161 asin = vectorMapF ASin
162 acos = vectorMapF ACos
163 atan = vectorMapF ATan
164 sinh = vectorMapF Sinh
165 cosh = vectorMapF Cosh
166 tanh = vectorMapF Tanh
167 asinh = vectorMapF ASinh
168 acosh = vectorMapF ACosh
169 atanh = vectorMapF ATanh
170 exp = vectorMapF Exp
171 log = vectorMapF Log
172 sqrt = vectorMapF Sqrt
173 (**) = adaptScalar (vectorMapValF PowSV) (vectorZipF Pow) (flip (vectorMapValF PowVS))
174 pi = fromList [pi]
175
176-------------------------------------------------------------
177
178instance Floating (Vector Double) where
179 sin = vectorMapR Sin
180 cos = vectorMapR Cos
181 tan = vectorMapR Tan
182 asin = vectorMapR ASin
183 acos = vectorMapR ACos
184 atan = vectorMapR ATan
185 sinh = vectorMapR Sinh
186 cosh = vectorMapR Cosh
187 tanh = vectorMapR Tanh
188 asinh = vectorMapR ASinh
189 acosh = vectorMapR ACosh
190 atanh = vectorMapR ATanh
191 exp = vectorMapR Exp
192 log = vectorMapR Log
193 sqrt = vectorMapR Sqrt
194 (**) = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS))
195 pi = fromList [pi]
196
197-------------------------------------------------------------
198
199instance Floating (Vector (Complex Double)) where
200 sin = vectorMapC Sin
201 cos = vectorMapC Cos
202 tan = vectorMapC Tan
203 asin = vectorMapC ASin
204 acos = vectorMapC ACos
205 atan = vectorMapC ATan
206 sinh = vectorMapC Sinh
207 cosh = vectorMapC Cosh
208 tanh = vectorMapC Tanh
209 asinh = vectorMapC ASinh
210 acosh = vectorMapC ACosh
211 atanh = vectorMapC ATanh
212 exp = vectorMapC Exp
213 log = vectorMapC Log
214 sqrt = vectorMapC Sqrt
215 (**) = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS))
216 pi = fromList [pi]
217
218-----------------------------------------------------------
219
220instance Floating (Vector (Complex Float)) where
221 sin = vectorMapQ Sin
222 cos = vectorMapQ Cos
223 tan = vectorMapQ Tan
224 asin = vectorMapQ ASin
225 acos = vectorMapQ ACos
226 atan = vectorMapQ ATan
227 sinh = vectorMapQ Sinh
228 cosh = vectorMapQ Cosh
229 tanh = vectorMapQ Tanh
230 asinh = vectorMapQ ASinh
231 acosh = vectorMapQ ACosh
232 atanh = vectorMapQ ATanh
233 exp = vectorMapQ Exp
234 log = vectorMapQ Log
235 sqrt = vectorMapQ Sqrt
236 (**) = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS))
237 pi = fromList [pi]
238
239-----------------------------------------------------------
240
241instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where
242 sin = liftMatrix sin
243 cos = liftMatrix cos
244 tan = liftMatrix tan
245 asin = liftMatrix asin
246 acos = liftMatrix acos
247 atan = liftMatrix atan
248 sinh = liftMatrix sinh
249 cosh = liftMatrix cosh
250 tanh = liftMatrix tanh
251 asinh = liftMatrix asinh
252 acosh = liftMatrix acosh
253 atanh = liftMatrix atanh
254 exp = liftMatrix exp
255 log = liftMatrix log
256 (**) = liftMatrix2Auto (**)
257 sqrt = liftMatrix sqrt
258 pi = (1><1) [pi]
259
260---------------------------------------------------------------
261
262-- instance (Storable a, Num (Vector a)) => Monoid (Vector a) where
263-- mempty = 0 { idim = 0 }
264-- mappend a b = mconcat [a,b]
265-- mconcat = j . filter ((>0).dim)
266-- where j [] = mempty
267-- j l = join l
268
269---------------------------------------------------------------
270
271-- instance (NFData a, Storable a) => NFData (Vector a) where
272-- rnf = rnf . (@>0)
273--
274-- instance (NFData a, Element a) => NFData (Matrix a) where
275-- rnf = rnf . flatten
276
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs
deleted file mode 100644
index fa3e209..0000000
--- a/lib/Numeric/LinearAlgebra/Interface.hs
+++ /dev/null
@@ -1,118 +0,0 @@
1{-# OPTIONS_GHC -fglasgow-exts #-}
2{-# LANGUAGE UndecidableInstances #-}
3-----------------------------------------------------------------------------
4{- |
5Module : Numeric.LinearAlgebra.Interface
6Copyright : (c) Alberto Ruiz 2007
7License : GPL-style
8
9Maintainer : Alberto Ruiz (aruiz at um dot es)
10Stability : provisional
11Portability : portable
12
13Some useful operators, and Show, Read, Eq, Num, Fractional, and Floating instances for Vector and Matrix.
14
15In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand.
16
17
18-}
19-----------------------------------------------------------------------------
20
21module Numeric.LinearAlgebra.Interface(
22 (<>),(<.>),
23 (<\>),
24 (.*),(*/),
25 (<|>),(<->),
26) where
27
28import Numeric.Vector
29import Numeric.Matrix
30import Numeric.LinearAlgebra.Algorithms
31import Numeric.LinearAlgebra.Linear
32
33class Mul a b c | a b -> c where
34 infixl 7 <>
35 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
36 (<>) :: Product t => a t -> b t -> c t
37
38instance Mul Matrix Matrix Matrix where
39 (<>) = mXm
40
41instance Mul Matrix Vector Vector where
42 (<>) m v = flatten $ m <> (asColumn v)
43
44instance Mul Vector Matrix Vector where
45 (<>) v m = flatten $ (asRow v) <> m
46
47---------------------------------------------------
48
49-- | Dot product: @u \<.\> v = dot u v@
50--(<.>) :: (Field t) => Vector t -> Vector t -> t
51(<.>) :: Vectors Vector t => Vector t -> Vector t -> t
52infixl 7 <.>
53(<.>) = dot
54
55----------------------------------------------------
56
57{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}
58
59-- -- | @x .* a = scale x a@
60-- (.*) :: (Linear c a) => a -> c a -> c a
61infixl 7 .*
62a .* x = scale a x
63
64----------------------------------------------------
65
66{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}
67
68-- -- | @a *\/ x = scale (recip x) a@
69-- (*/) :: (Linear c a) => c a -> a -> c a
70infixl 7 */
71v */ x = scale (recip x) v
72
73-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).
74(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a
75infixl 7 <\>
76m <\> v = flatten (linearSolveSVD m (reshape 1 v))
77
78------------------------------------------------
79
80{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}
81{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}
82
83class Joinable a b where
84 joinH :: Element t => a t -> b t -> Matrix t
85 joinV :: Element t => a t -> b t -> Matrix t
86
87instance Joinable Matrix Matrix where
88 joinH m1 m2 = fromBlocks [[m1,m2]]
89 joinV m1 m2 = fromBlocks [[m1],[m2]]
90
91instance Joinable Matrix Vector where
92 joinH m v = joinH m (asColumn v)
93 joinV m v = joinV m (asRow v)
94
95instance Joinable Vector Matrix where
96 joinH v m = joinH (asColumn v) m
97 joinV v m = joinV (asRow v) m
98
99infixl 4 <|>
100infixl 3 <->
101
102{-- - | Horizontal concatenation of matrices and vectors:
103
104@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]
105(6><4)
106 [ 1.0, 0.0, 0.0, 1.0
107 , 0.0, 1.0, 0.0, 2.0
108 , 0.0, 0.0, 1.0, 3.0
109 , 3.0, 0.0, 0.0, 4.0
110 , 0.0, 3.0, 0.0, 5.0
111 , 0.0, 0.0, 3.0, 6.0 ]@
112-}
113-- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
114a <|> b = joinH a b
115
116-- -- | Vertical concatenation of matrices and vectors.
117-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
118a <-> b = joinV a b
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index 5d0154d..288439f 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -44,8 +44,7 @@ module Numeric.LinearAlgebra.LAPACK (
44import Data.Packed.Internal 44import Data.Packed.Internal
45import Data.Packed.Matrix 45import Data.Packed.Matrix
46import Data.Complex 46import Data.Complex
47import Numeric.Vector() 47import Numeric.Conversion
48import Numeric.Container
49import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) 48import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale))
50import Foreign 49import Foreign
51import Foreign.C.Types (CInt) 50import Foreign.C.Types (CInt)
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs
deleted file mode 100644
index 775060e..0000000
--- a/lib/Numeric/LinearAlgebra/Linear.hs
+++ /dev/null
@@ -1,160 +0,0 @@
1{-# LANGUAGE UndecidableInstances, MultiParamTypeClasses, FlexibleInstances #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE TypeFamilies #-}
4-----------------------------------------------------------------------------
5{- |
6Module : Numeric.LinearAlgebra.Linear
7Copyright : (c) Alberto Ruiz 2006-7
8License : GPL-style
9
10Maintainer : Alberto Ruiz (aruiz at um dot es)
11Stability : provisional
12Portability : uses ffi
13
14Basic optimized operations on vectors and matrices.
15
16-}
17-----------------------------------------------------------------------------
18
19module Numeric.LinearAlgebra.Linear (
20 -- * Linear Algebra Typeclasses
21 Vectors(..),
22 -- * Products
23 Product(..),
24 mXm,mXv,vXm,
25 outer, kronecker,
26 -- * Modules
27 --module Numeric.Vector,
28 --module Numeric.Matrix,
29 module Numeric.Container
30) where
31
32import Data.Packed.Internal.Common
33import Data.Packed.Matrix
34import Data.Packed.Vector
35import Data.Complex
36import Numeric.Container
37import Numeric.Vector()
38import Numeric.Matrix()
39import Numeric.GSL.Vector
40import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
41
42-- | Linear algebraic properties of objects
43class Num e => Vectors a e where
44 -- | dot (inner) product
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
50 -- | euclidean norm
51 norm2 :: a e -> e
52 -- | element of maximum magnitude
53 normInf :: a e -> e
54
55instance 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
62instance 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
69instance 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
76instance 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
83----------------------------------------------------
84
85class Element t => Product t where
86 multiply :: Matrix t -> Matrix t -> Matrix t
87 ctrans :: Matrix t -> Matrix t
88
89instance Product Double where
90 multiply = multiplyR
91 ctrans = trans
92
93instance Product (Complex Double) where
94 multiply = multiplyC
95 ctrans = conj . trans
96
97instance Product Float where
98 multiply = multiplyF
99 ctrans = trans
100
101instance Product (Complex Float) where
102 multiply = multiplyQ
103 ctrans = conj . trans
104
105----------------------------------------------------------
106
107-- synonym for matrix product
108mXm :: Product t => Matrix t -> Matrix t -> Matrix t
109mXm = multiply
110
111-- matrix - vector product
112mXv :: Product t => Matrix t -> Vector t -> Vector t
113mXv m v = flatten $ m `mXm` (asColumn v)
114
115-- vector - matrix product
116vXm :: Product t => Vector t -> Matrix t -> Vector t
117vXm v m = flatten $ (asRow v) `mXm` m
118
119{- | Outer product of two vectors.
120
121@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
122(3><3)
123 [ 5.0, 2.0, 3.0
124 , 10.0, 4.0, 6.0
125 , 15.0, 6.0, 9.0 ]@
126-}
127outer :: (Product t) => Vector t -> Vector t -> Matrix t
128outer u v = asColumn u `multiply` asRow v
129
130{- | Kronecker product of two matrices.
131
132@m1=(2><3)
133 [ 1.0, 2.0, 0.0
134 , 0.0, -1.0, 3.0 ]
135m2=(4><3)
136 [ 1.0, 2.0, 3.0
137 , 4.0, 5.0, 6.0
138 , 7.0, 8.0, 9.0
139 , 10.0, 11.0, 12.0 ]@
140
141@\> kronecker m1 m2
142(8><9)
143 [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0
144 , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0
145 , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0
146 , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0
147 , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0
148 , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0
149 , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0
150 , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
151-}
152kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
153kronecker a b = fromBlocks
154 . splitEvery (cols a)
155 . map (reshape (cols b))
156 . toRows
157 $ flatten a `outer` flatten b
158
159
160-------------------------------------------------------------------
diff --git a/lib/Numeric/Matrix.hs b/lib/Numeric/Matrix.hs
index fa3f94a..6ba870f 100644
--- a/lib/Numeric/Matrix.hs
+++ b/lib/Numeric/Matrix.hs
@@ -3,6 +3,7 @@
3{-# LANGUAGE FlexibleInstances #-} 3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE UndecidableInstances #-} 4{-# LANGUAGE UndecidableInstances #-}
5{-# LANGUAGE MultiParamTypeClasses #-} 5{-# LANGUAGE MultiParamTypeClasses #-}
6{-# LANGUAGE FunctionalDependencies #-}
6 7
7----------------------------------------------------------------------------- 8-----------------------------------------------------------------------------
8-- | 9-- |
@@ -14,12 +15,22 @@
14-- Stability : provisional 15-- Stability : provisional
15-- Portability : portable 16-- Portability : portable
16-- 17--
17-- Numeric instances and functions for 'Data.Packed.Matrix's 18-- Numeric instances and functions for 'Matrix'.
19-- In the context of the standard numeric operators, one-component
20-- vectors and matrices automatically expand to match the dimensions of the other operand.
18-- 21--
19----------------------------------------------------------------------------- 22-----------------------------------------------------------------------------
20 23
21module Numeric.Matrix ( 24module Numeric.Matrix (
22 module Data.Packed.Matrix, 25 -- * Basic functions
26 module Data.Packed.Matrix,
27 module Numeric.Vector,
28 --module Numeric.Container,
29 -- * Operators
30 (<>), (<\>),
31 -- * Deprecated
32 (.*),(*/),(<|>),(<->),
33 vectorMax,vectorMin
23 ) where 34 ) where
24 35
25------------------------------------------------------------------- 36-------------------------------------------------------------------
@@ -28,18 +39,18 @@ import Data.Packed.Vector
28import Data.Packed.Matrix 39import Data.Packed.Matrix
29import Numeric.Container 40import Numeric.Container
30--import Numeric.LinearAlgebra.Linear 41--import Numeric.LinearAlgebra.Linear
31import Numeric.Vector() 42import Numeric.Vector
43import Numeric.LinearAlgebra.Algorithms
44--import Control.Monad(ap)
32 45
33import Control.Monad(ap) 46--import Control.Arrow((***))
34
35import Control.Arrow((***))
36 47
37------------------------------------------------------------------- 48-------------------------------------------------------------------
38 49
39instance Linear Matrix a => Eq (Matrix a) where 50instance Container Matrix a => Eq (Matrix a) where
40 (==) = equal 51 (==) = equal
41 52
42instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where 53instance (Container Matrix a, Num (Vector a)) => Num (Matrix a) where
43 (+) = liftMatrix2Auto (+) 54 (+) = liftMatrix2Auto (+)
44 (-) = liftMatrix2Auto (-) 55 (-) = liftMatrix2Auto (-)
45 negate = liftMatrix negate 56 negate = liftMatrix negate
@@ -50,13 +61,13 @@ instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where
50 61
51--------------------------------------------------- 62---------------------------------------------------
52 63
53instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where 64instance (Container Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where
54 fromRational n = (1><1) [fromRational n] 65 fromRational n = (1><1) [fromRational n]
55 (/) = liftMatrix2Auto (/) 66 (/) = liftMatrix2Auto (/)
56 67
57--------------------------------------------------------- 68---------------------------------------------------------
58 69
59instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where 70instance (Container Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where
60 sin = liftMatrix sin 71 sin = liftMatrix sin
61 cos = liftMatrix cos 72 cos = liftMatrix cos
62 tan = liftMatrix tan 73 tan = liftMatrix tan
@@ -75,43 +86,93 @@ instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floati
75 sqrt = liftMatrix sqrt 86 sqrt = liftMatrix sqrt
76 pi = (1><1) [pi] 87 pi = (1><1) [pi]
77 88
78--------------------------------------------------------------- 89--------------------------------------------------------
79 90
80instance NumericContainer Matrix where 91class Mul a b c | a b -> c where
81 toComplex = uncurry $ liftMatrix2 $ curry toComplex 92 infixl 7 <>
82 fromComplex z = (reshape c *** reshape c) . fromComplex . flatten $ z 93 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
83 where c = cols z 94 (<>) :: Product t => a t -> b t -> c t
84 complex' = liftMatrix complex' 95
85 conj = liftMatrix conj 96instance Mul Matrix Matrix Matrix where
86-- cmap f = liftMatrix (cmap f) 97 (<>) = mXm
87 single' = liftMatrix single' 98
88 double' = liftMatrix double' 99instance Mul Matrix Vector Vector where
89 100 (<>) m v = flatten $ m <> (asColumn v)
90--------------------------------------------------------------- 101
91 102instance Mul Vector Matrix Vector where
92instance (Linear Vector a, Container Matrix a) => Linear Matrix a where 103 (<>) v m = flatten $ (asRow v) <> m
93 scale x = liftMatrix (scale x)
94 scaleRecip x = liftMatrix (scaleRecip x)
95 addConstant x = liftMatrix (addConstant x)
96 add = liftMatrix2 add
97 sub = liftMatrix2 sub
98 mul = liftMatrix2 mul
99 divide = liftMatrix2 divide
100 equal a b = cols a == cols b && flatten a `equal` flatten b
101 scalar x = (1><1) [x]
102 --
103instance (Container Vector a) => Container Matrix a where
104 cmap f = liftMatrix (mapVector f)
105 atIndex = (@@>)
106 minIndex m = let (r,c) = (rows m,cols m)
107 i = (minIndex $ flatten m)
108 in (i `div` c,(i `mod` c) + 1)
109 maxIndex m = let (r,c) = (rows m,cols m)
110 i = (maxIndex $ flatten m)
111 in (i `div` c,(i `mod` c) + 1)
112 minElement = ap (@@>) minIndex
113 maxElement = ap (@@>) maxIndex
114 sumElements = sumElements . flatten
115 prodElements = prodElements . flatten
116 104
117---------------------------------------------------- 105----------------------------------------------------
106
107{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}
108
109-- -- | @x .* a = scale x a@
110-- (.*) :: (Linear c a) => a -> c a -> c a
111infixl 7 .*
112a .* x = scale a x
113
114----------------------------------------------------
115
116{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}
117
118-- -- | @a *\/ x = scale (recip x) a@
119-- (*/) :: (Linear c a) => c a -> a -> c a
120infixl 7 */
121v */ x = scale (recip x) v
122
123-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).
124(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a
125infixl 7 <\>
126m <\> v = flatten (linearSolveSVD m (reshape 1 v))
127
128------------------------------------------------
129
130{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}
131{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}
132
133class Joinable a b where
134 joinH :: Element t => a t -> b t -> Matrix t
135 joinV :: Element t => a t -> b t -> Matrix t
136
137instance Joinable Matrix Matrix where
138 joinH m1 m2 = fromBlocks [[m1,m2]]
139 joinV m1 m2 = fromBlocks [[m1],[m2]]
140
141instance Joinable Matrix Vector where
142 joinH m v = joinH m (asColumn v)
143 joinV m v = joinV m (asRow v)
144
145instance Joinable Vector Matrix where
146 joinH v m = joinH (asColumn v) m
147 joinV v m = joinV (asRow v) m
148
149infixl 4 <|>
150infixl 3 <->
151
152{-- - | Horizontal concatenation of matrices and vectors:
153
154@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]
155(6><4)
156 [ 1.0, 0.0, 0.0, 1.0
157 , 0.0, 1.0, 0.0, 2.0
158 , 0.0, 0.0, 1.0, 3.0
159 , 3.0, 0.0, 0.0, 4.0
160 , 0.0, 3.0, 0.0, 5.0
161 , 0.0, 0.0, 3.0, 6.0 ]@
162-}
163-- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
164a <|> b = joinH a b
165
166-- -- | Vertical concatenation of matrices and vectors.
167-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
168a <-> b = joinV a b
169
170-------------------------------------------------------------------
171
172{-# DEPRECATED vectorMin "use minElement" #-}
173vectorMin :: (Container Vector t, Element t) => Vector t -> t
174vectorMin = minElement
175
176{-# DEPRECATED vectorMax "use maxElement" #-}
177vectorMax :: (Container Vector t, Element t) => Vector t -> t
178vectorMax = maxElement
diff --git a/lib/Numeric/Vector.hs b/lib/Numeric/Vector.hs
index 55d645a..9427243 100644
--- a/lib/Numeric/Vector.hs
+++ b/lib/Numeric/Vector.hs
@@ -14,24 +14,26 @@
14-- Stability : provisional 14-- Stability : provisional
15-- Portability : portable 15-- Portability : portable
16-- 16--
17-- Numeric instances and functions for 'Data.Packed.Vector's 17-- Numeric instances and functions for 'Vector'.
18-- 18--
19----------------------------------------------------------------------------- 19-----------------------------------------------------------------------------
20 20
21module Numeric.Vector ( 21module Numeric.Vector (
22 -- * Vector creation 22 -- * Basic vector functions
23 constant, linspace, 23 module Data.Packed.Vector,
24 module Data.Packed.Vector 24 module Numeric.Container,
25 -- * Vector creation
26 constant, linspace,
27 -- * Operators
28 (<.>)
25 ) where 29 ) where
26 30
27import Data.Complex 31import Data.Complex
28 32
29import Control.Monad(ap)
30
31import Data.Packed.Vector 33import Data.Packed.Vector
32import Data.Packed.Internal.Matrix(Element(..)) 34import Data.Packed.Internal.Matrix
33import Data.Packed.Internal.Vector(asComplex,asReal) 35--import Data.Packed.Internal.Vector(asComplex,asReal)
34import Data.Packed.Matrix(toColumns,fromColumns,flatten,reshape) 36--import Data.Packed.Matrix(toColumns,fromColumns,flatten,reshape)
35import Numeric.GSL.Vector 37import Numeric.GSL.Vector
36 38
37import Numeric.Container 39import Numeric.Container
@@ -92,10 +94,15 @@ Logarithmic spacing can be defined as follows:
92 94
93@logspace n (a,b) = 10 ** linspace n (a,b)@ 95@logspace n (a,b) = 10 ** linspace n (a,b)@
94-} 96-}
95linspace :: (Enum e, Linear Vector e) => Int -> (e, e) -> Vector e 97linspace :: (Enum e, Container Vector e) => Int -> (e, e) -> Vector e
96linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1] 98linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1]
97 where s = (b-a)/fromIntegral (n-1) 99 where s = (b-a)/fromIntegral (n-1)
98 100
101-- | Dot product: @u \<.\> v = dot u v@
102(<.>) :: Vectors Vector t => Vector t -> Vector t -> t
103infixl 7 <.>
104(<.>) = dot
105
99------------------------------------------------------------------ 106------------------------------------------------------------------
100 107
101adaptScalar f1 f2 f3 x y 108adaptScalar f1 f2 f3 x y
@@ -107,7 +114,7 @@ adaptScalar f1 f2 f3 x y
107 114
108#ifndef VECTOR 115#ifndef VECTOR
109 116
110instance Linear Vector a => Eq (Vector a) where 117instance Container Vector a => Eq (Vector a) where
111 (==) = equal 118 (==) = equal
112 119
113#endif 120#endif
@@ -146,7 +153,7 @@ instance Num (Vector (Complex Float)) where
146 153
147--------------------------------------------------- 154---------------------------------------------------
148 155
149instance (Linear Vector a, Num (Vector a)) => Fractional (Vector a) where 156instance (Container Vector a, Num (Vector a)) => Fractional (Vector a) where
150 fromRational n = fromList [fromRational n] 157 fromRational n = fromList [fromRational n]
151 (/) = adaptScalar f divide g where 158 (/) = adaptScalar f divide g where
152 r `f` v = scaleRecip r v 159 r `f` v = scaleRecip r v
@@ -254,116 +261,6 @@ instance Floating (Vector (Complex Float)) where
254-- instance (NFData a, Element a) => NFData (Matrix a) where 261-- instance (NFData a, Element a) => NFData (Matrix a) where
255-- rnf = rnf . flatten 262-- rnf = rnf . flatten
256 263
257---------------------------------------------------------------
258
259-- | obtains the complex conjugate of a complex vector
260conjV :: (RealElement a) => Vector (Complex a) -> Vector (Complex a)
261conjV = mapVector conjugate
262
263-- | creates a complex vector from vectors with real and imaginary parts
264toComplexV :: (RealElement a) => (Vector a, Vector a) -> Vector (Complex a)
265toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i]
266
267-- | the inverse of 'toComplex'
268fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a)
269fromComplexV z = (r,i) where
270 [r,i] = toColumns $ reshape 2 $ asReal z
271 264
272-------------------------------------------------------------------------- 265--------------------------------------------------------------------------
273 266
274instance 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
285instance Linear Vector Float where
286 scale = vectorMapValF Scale
287 scaleRecip = vectorMapValF Recip
288 addConstant = vectorMapValF AddConstant
289 add = vectorZipF Add
290 sub = vectorZipF Sub
291 mul = vectorZipF Mul
292 divide = vectorZipF Div
293 equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
294 scalar x = fromList [x]
295 --
296instance Container Vector Float where
297 cmap = mapVector
298 atIndex = (@>)
299 minIndex = round . toScalarF MinIdx
300 maxIndex = round . toScalarF MaxIdx
301 minElement = toScalarF Min
302 maxElement = toScalarF Max
303 sumElements = sumF
304 prodElements = prodF
305
306instance Linear Vector Double where
307 scale = vectorMapValR Scale
308 scaleRecip = vectorMapValR Recip
309 addConstant = vectorMapValR AddConstant
310 add = vectorZipR Add
311 sub = vectorZipR Sub
312 mul = vectorZipR Mul
313 divide = vectorZipR Div
314 equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
315 scalar x = fromList [x]
316 --
317instance Container Vector Double where
318 cmap = mapVector
319 atIndex = (@>)
320 minIndex = round . toScalarR MinIdx
321 maxIndex = round . toScalarR MaxIdx
322 minElement = toScalarR Min
323 maxElement = toScalarR Max
324 sumElements = sumR
325 prodElements = prodR
326
327instance Linear Vector (Complex Double) where
328 scale = vectorMapValC Scale
329 scaleRecip = vectorMapValC Recip
330 addConstant = vectorMapValC AddConstant
331 add = vectorZipC Add
332 sub = vectorZipC Sub
333 mul = vectorZipC Mul
334 divide = vectorZipC Div
335 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
336 scalar x = fromList [x]
337 --
338instance Container Vector (Complex Double) where
339 cmap = mapVector
340 atIndex = (@>)
341 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
342 maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
343 minElement = ap (@>) minIndex
344 maxElement = ap (@>) maxIndex
345 sumElements = sumC
346 prodElements = prodC
347
348instance Linear Vector (Complex Float) where
349 scale = vectorMapValQ Scale
350 scaleRecip = vectorMapValQ Recip
351 addConstant = vectorMapValQ AddConstant
352 add = vectorZipQ Add
353 sub = vectorZipQ Sub
354 mul = vectorZipQ Mul
355 divide = vectorZipQ Div
356 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
357 scalar x = fromList [x]
358 --
359instance Container Vector (Complex Float) where
360 cmap = mapVector
361 atIndex = (@>)
362 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
363 maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
364 minElement = ap (@>) minIndex
365 maxElement = ap (@>) maxIndex
366 sumElements = sumQ
367 prodElements = prodQ
368
369---------------------------------------------------------------