summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--hmatrix.cabal21
-rw-r--r--lib/Data/Packed/Random.hs2
-rw-r--r--lib/Graphics/Plot.hs2
-rw-r--r--lib/Numeric/Chain.hs2
-rw-r--r--lib/Numeric/Container.hs557
-rw-r--r--lib/Numeric/ContainerBoot.hs575
-rw-r--r--lib/Numeric/LinearAlgebra.hs19
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs15
-rw-r--r--lib/Numeric/Matrix.hs53
-rw-r--r--lib/Numeric/MatrixBoot.hs41
-rw-r--r--lib/Numeric/Vector.hs42
11 files changed, 674 insertions, 655 deletions
diff --git a/hmatrix.cabal b/hmatrix.cabal
index ae03b0f..fb5ed05 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -93,31 +93,28 @@ library
93 Numeric.GSL.Root, 93 Numeric.GSL.Root,
94 Numeric.GSL.Fitting, 94 Numeric.GSL.Fitting,
95 Numeric.GSL.ODE, 95 Numeric.GSL.ODE,
96 Numeric.GSL.Vector,
97 Numeric.GSL, 96 Numeric.GSL,
98 Numeric.Vector,
99 Numeric.Matrix,
100 Numeric.Container, 97 Numeric.Container,
101 Numeric.LinearAlgebra, 98 Numeric.LinearAlgebra,
102 Numeric.LinearAlgebra.LAPACK, 99 Numeric.LinearAlgebra.LAPACK,
103 --Numeric.LinearAlgebra.Interface,
104 --Numeric.LinearAlgebra.Linear,
105 Numeric.LinearAlgebra.Algorithms, 100 Numeric.LinearAlgebra.Algorithms,
106 Graphics.Plot, 101 Graphics.Plot,
107 -- Data.Packed.Convert,
108 Data.Packed.ST, 102 Data.Packed.ST,
109 Data.Packed.Development, 103 Data.Packed.Development
110 Data.Packed.Random
111 other-modules: Data.Packed.Internal, 104 other-modules: Data.Packed.Internal,
112 Data.Packed.Internal.Common, 105 Data.Packed.Internal.Common,
113 Data.Packed.Internal.Signatures, 106 Data.Packed.Internal.Signatures,
114 Data.Packed.Internal.Vector, 107 Data.Packed.Internal.Vector,
115 Data.Packed.Internal.Matrix, 108 Data.Packed.Internal.Matrix,
109 Data.Packed.Random,
116 Numeric.GSL.Internal, 110 Numeric.GSL.Internal,
117 Numeric.Conversion 111 Numeric.GSL.Vector,
118 Numeric.MatrixBoot 112 Numeric.Conversion,
119 Numeric.IO 113 Numeric.ContainerBoot,
120 Numeric.Chain 114 Numeric.IO,
115 Numeric.Chain,
116 Numeric.Vector,
117 Numeric.Matrix
121 118
122 C-sources: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c, 119 C-sources: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c,
123 lib/Numeric/GSL/gsl-aux.c 120 lib/Numeric/GSL/gsl-aux.c
diff --git a/lib/Data/Packed/Random.hs b/lib/Data/Packed/Random.hs
index c34f212..14d91ea 100644
--- a/lib/Data/Packed/Random.hs
+++ b/lib/Data/Packed/Random.hs
@@ -21,7 +21,7 @@ module Data.Packed.Random (
21 21
22import Numeric.GSL.Vector 22import Numeric.GSL.Vector
23import Data.Packed 23import Data.Packed
24import Numeric.Container 24import Numeric.ContainerBoot
25import Numeric.LinearAlgebra.Algorithms 25import Numeric.LinearAlgebra.Algorithms
26 26
27 27
diff --git a/lib/Graphics/Plot.hs b/lib/Graphics/Plot.hs
index 0bdd803..32106df 100644
--- a/lib/Graphics/Plot.hs
+++ b/lib/Graphics/Plot.hs
@@ -26,7 +26,7 @@ module Graphics.Plot(
26 26
27) where 27) where
28 28
29import Numeric.Matrix 29import Numeric.Container
30import Data.List(intersperse) 30import Data.List(intersperse)
31import System.Process (system) 31import System.Process (system)
32 32
diff --git a/lib/Numeric/Chain.hs b/lib/Numeric/Chain.hs
index 299d8fa..03ca88d 100644
--- a/lib/Numeric/Chain.hs
+++ b/lib/Numeric/Chain.hs
@@ -19,7 +19,7 @@ module Numeric.Chain (
19import Data.Maybe 19import Data.Maybe
20 20
21import Data.Packed.Matrix 21import Data.Packed.Matrix
22import Numeric.Container 22import Numeric.ContainerBoot
23 23
24import qualified Data.Array.IArray as A 24import qualified Data.Array.IArray as A
25 25
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs
index 6b73a4e..b3f46de 100644
--- a/lib/Numeric/Container.hs
+++ b/lib/Numeric/Container.hs
@@ -26,13 +26,24 @@
26----------------------------------------------------------------------------- 26-----------------------------------------------------------------------------
27 27
28module Numeric.Container ( 28module Numeric.Container (
29 -- * Basic functions
30 module Data.Packed,
31 constant, linspace,
32 diag, ident,
33 ctrans,
29 -- * Generic operations 34 -- * Generic operations
30 Container(..), 35 Container(..),
31 -- * Matrix product and related functions 36 -- * Matrix product
32 Product(..), 37 Product(..),
33 mXm,mXv,vXm, 38 optimiseMult,
39 mXm,mXv,vXm,(<.>),(<>),(<\>),
34 outer, kronecker, 40 outer, kronecker,
35 41 -- * Random numbers
42 RandDist(..),
43 randomVector,
44 gaussianSample,
45 uniformSample,
46 meanCov,
36 -- * Element conversion 47 -- * Element conversion
37 Convert(..), 48 Convert(..),
38 Complexable(), 49 Complexable(),
@@ -42,6 +53,11 @@ module Numeric.Container (
42 53
43 IndexOf, 54 IndexOf,
44 module Data.Complex, 55 module Data.Complex,
56 -- * Input / Output
57 dispf, disps, dispcf, vecdisp, latexFormat, format,
58 loadMatrix, saveMatrix, fromFile, fileDimensions,
59 readMatrix,
60 fscanfVector, fprintfVector, freadVector, fwriteVector,
45 -- * Experimental 61 -- * Experimental
46 build', konst', 62 build', konst',
47 -- * Deprecated 63 -- * Deprecated
@@ -51,517 +67,64 @@ module Numeric.Container (
51) where 67) where
52 68
53import Data.Packed 69import Data.Packed
54import Numeric.Conversion 70import Data.Packed.Internal(constantD)
55import Data.Packed.Internal 71import Numeric.ContainerBoot
56import Numeric.GSL.Vector 72import Numeric.Chain
57 73import Numeric.IO
58import Data.Complex 74import Data.Complex
59import Control.Monad(ap) 75import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
60 76import Data.Packed.Random
61import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
62
63import System.IO.Unsafe
64
65-------------------------------------------------------------------
66
67type family IndexOf c
68
69type instance IndexOf Vector = Int
70type instance IndexOf Matrix = (Int,Int)
71
72type family ArgOf c a
73
74type instance ArgOf Vector a = a -> a
75type instance ArgOf Matrix a = a -> a -> a
76
77-------------------------------------------------------------------
78
79-- | Basic element-by-element functions for numeric containers
80class (Complexable c, Fractional e, Element e) => Container c e where
81 -- | create a structure with a single element
82 scalar :: e -> c e
83 -- | complex conjugate
84 conj :: c e -> c e
85 scale :: e -> c e -> c e
86 -- | scale the element by element reciprocal of the object:
87 --
88 -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
89 scaleRecip :: e -> c e -> c e
90 addConstant :: e -> c e -> c e
91 add :: c e -> c e -> c e
92 sub :: c e -> c e -> c e
93 -- | element by element multiplication
94 mul :: c e -> c e -> c e
95 -- | element by element division
96 divide :: c e -> c e -> c e
97 equal :: c e -> c e -> Bool
98 --
99 -- | cannot implement instance Functor because of Element class constraint
100 cmap :: (Element a, Element b) => (a -> b) -> c a -> c b
101 -- | constant structure of given size
102 konst :: e -> IndexOf c -> c e
103 -- | create a structure using a function
104 --
105 -- Hilbert matrix of order N:
106 --
107 -- @hilb n = build (n,n) (\\i j -> 1/(i+j+1))@
108 build :: IndexOf c -> (ArgOf c e) -> c e
109 --build :: BoundsOf f -> f -> (ContainerOf f) e
110 --
111 -- | indexing function
112 atIndex :: c e -> IndexOf c -> e
113 -- | index of min element
114 minIndex :: c e -> IndexOf c
115 -- | index of max element
116 maxIndex :: c e -> IndexOf c
117 -- | value of min element
118 minElement :: c e -> e
119 -- | value of max element
120 maxElement :: c e -> e
121 -- the C functions sumX/prodX are twice as fast as using foldVector
122 -- | the sum of elements (faster than using @fold@)
123 sumElements :: c e -> e
124 -- | the product of elements (faster than using @fold@)
125 prodElements :: c e -> e
126
127--------------------------------------------------------------------------
128
129instance Container Vector Float where
130 scale = vectorMapValF Scale
131 scaleRecip = vectorMapValF Recip
132 addConstant = vectorMapValF AddConstant
133 add = vectorZipF Add
134 sub = vectorZipF Sub
135 mul = vectorZipF Mul
136 divide = vectorZipF Div
137 equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
138 scalar x = fromList [x]
139 konst = constantD
140 build = buildV
141 conj = id
142 cmap = mapVector
143 atIndex = (@>)
144 minIndex = round . toScalarF MinIdx
145 maxIndex = round . toScalarF MaxIdx
146 minElement = toScalarF Min
147 maxElement = toScalarF Max
148 sumElements = sumF
149 prodElements = prodF
150
151instance Container Vector Double where
152 scale = vectorMapValR Scale
153 scaleRecip = vectorMapValR Recip
154 addConstant = vectorMapValR AddConstant
155 add = vectorZipR Add
156 sub = vectorZipR Sub
157 mul = vectorZipR Mul
158 divide = vectorZipR Div
159 equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
160 scalar x = fromList [x]
161 konst = constantD
162 build = buildV
163 conj = id
164 cmap = mapVector
165 atIndex = (@>)
166 minIndex = round . toScalarR MinIdx
167 maxIndex = round . toScalarR MaxIdx
168 minElement = toScalarR Min
169 maxElement = toScalarR Max
170 sumElements = sumR
171 prodElements = prodR
172
173instance Container Vector (Complex Double) where
174 scale = vectorMapValC Scale
175 scaleRecip = vectorMapValC Recip
176 addConstant = vectorMapValC AddConstant
177 add = vectorZipC Add
178 sub = vectorZipC Sub
179 mul = vectorZipC Mul
180 divide = vectorZipC Div
181 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
182 scalar x = fromList [x]
183 konst = constantD
184 build = buildV
185 conj = conjugateC
186 cmap = mapVector
187 atIndex = (@>)
188 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
189 maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
190 minElement = ap (@>) minIndex
191 maxElement = ap (@>) maxIndex
192 sumElements = sumC
193 prodElements = prodC
194
195instance Container Vector (Complex Float) where
196 scale = vectorMapValQ Scale
197 scaleRecip = vectorMapValQ Recip
198 addConstant = vectorMapValQ AddConstant
199 add = vectorZipQ Add
200 sub = vectorZipQ Sub
201 mul = vectorZipQ Mul
202 divide = vectorZipQ Div
203 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
204 scalar x = fromList [x]
205 konst = constantD
206 build = buildV
207 conj = conjugateQ
208 cmap = mapVector
209 atIndex = (@>)
210 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
211 maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
212 minElement = ap (@>) minIndex
213 maxElement = ap (@>) maxIndex
214 sumElements = sumQ
215 prodElements = prodQ
216
217---------------------------------------------------------------
218
219instance (Container Vector a) => Container Matrix a where
220 scale x = liftMatrix (scale x)
221 scaleRecip x = liftMatrix (scaleRecip x)
222 addConstant x = liftMatrix (addConstant x)
223 add = liftMatrix2 add
224 sub = liftMatrix2 sub
225 mul = liftMatrix2 mul
226 divide = liftMatrix2 divide
227 equal a b = cols a == cols b && flatten a `equal` flatten b
228 scalar x = (1><1) [x]
229 konst v (r,c) = reshape c (konst v (r*c))
230 build = buildM
231 conj = liftMatrix conj
232 cmap f = liftMatrix (mapVector f)
233 atIndex = (@@>)
234 minIndex m = let (r,c) = (rows m,cols m)
235 i = (minIndex $ flatten m)
236 in (i `div` c,i `mod` c)
237 maxIndex m = let (r,c) = (rows m,cols m)
238 i = (maxIndex $ flatten m)
239 in (i `div` c,i `mod` c)
240 minElement = ap (@@>) minIndex
241 maxElement = ap (@@>) maxIndex
242 sumElements = sumElements . flatten
243 prodElements = prodElements . flatten
244
245----------------------------------------------------
246
247-- | Matrix product and related functions
248class Element e => Product e where
249 -- | matrix product
250 multiply :: Matrix e -> Matrix e -> Matrix e
251 -- | dot (inner) product
252 dot :: Vector e -> Vector e -> e
253 -- | sum of absolute value of elements (differs in complex case from @norm1@)
254 absSum :: Vector e -> RealOf e
255 -- | sum of absolute value of elements
256 norm1 :: Vector e -> RealOf e
257 -- | euclidean norm
258 norm2 :: Vector e -> RealOf e
259 -- | element of maximum magnitude
260 normInf :: Vector e -> RealOf e
261
262instance Product Float where
263 norm2 = toScalarF Norm2
264 absSum = toScalarF AbsSum
265 dot = dotF
266 norm1 = toScalarF AbsSum
267 normInf = maxElement . vectorMapF Abs
268 multiply = multiplyF
269 77
270instance Product Double where 78------------------------------------------------------------------
271 norm2 = toScalarR Norm2
272 absSum = toScalarR AbsSum
273 dot = dotR
274 norm1 = toScalarR AbsSum
275 normInf = maxElement . vectorMapR Abs
276 multiply = multiplyR
277 79
278instance Product (Complex Float) where 80{- | creates a vector with a given number of equal components:
279 norm2 = toScalarQ Norm2
280 absSum = toScalarQ AbsSum
281 dot = dotQ
282 norm1 = sumElements . fst . fromComplex . vectorMapQ Abs
283 normInf = maxElement . fst . fromComplex . vectorMapQ Abs
284 multiply = multiplyQ
285 81
286instance Product (Complex Double) where 82@> constant 2 7
287 norm2 = toScalarC Norm2 837 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
288 absSum = toScalarC AbsSum
289 dot = dotC
290 norm1 = sumElements . fst . fromComplex . vectorMapC Abs
291 normInf = maxElement . fst . fromComplex . vectorMapC Abs
292 multiply = multiplyC
293
294----------------------------------------------------------
295
296-- synonym for matrix product
297mXm :: Product t => Matrix t -> Matrix t -> Matrix t
298mXm = multiply
299
300-- matrix - vector product
301mXv :: Product t => Matrix t -> Vector t -> Vector t
302mXv m v = flatten $ m `mXm` (asColumn v)
303
304-- vector - matrix product
305vXm :: Product t => Vector t -> Matrix t -> Vector t
306vXm v m = flatten $ (asRow v) `mXm` m
307
308{- | Outer product of two vectors.
309
310@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
311(3><3)
312 [ 5.0, 2.0, 3.0
313 , 10.0, 4.0, 6.0
314 , 15.0, 6.0, 9.0 ]@
315-} 84-}
316outer :: (Product t) => Vector t -> Vector t -> Matrix t 85constant :: Element a => a -> Int -> Vector a
317outer u v = asColumn u `multiply` asRow v 86-- constant x n = runSTVector (newVector x n)
318 87constant = constantD-- about 2x faster
319{- | Kronecker product of two matrices.
320
321@m1=(2><3)
322 [ 1.0, 2.0, 0.0
323 , 0.0, -1.0, 3.0 ]
324m2=(4><3)
325 [ 1.0, 2.0, 3.0
326 , 4.0, 5.0, 6.0
327 , 7.0, 8.0, 9.0
328 , 10.0, 11.0, 12.0 ]@
329
330@\> kronecker m1 m2
331(8><9)
332 [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0
333 , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0
334 , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0
335 , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0
336 , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0
337 , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0
338 , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0
339 , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
340-}
341kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
342kronecker a b = fromBlocks
343 . splitEvery (cols a)
344 . map (reshape (cols b))
345 . toRows
346 $ flatten a `outer` flatten b
347
348-------------------------------------------------------------------
349
350
351class Convert t where
352 real :: Container c t => c (RealOf t) -> c t
353 complex :: Container c t => c t -> c (ComplexOf t)
354 single :: Container c t => c t -> c (SingleOf t)
355 double :: Container c t => c t -> c (DoubleOf t)
356 toComplex :: (Container c t, RealElement t) => (c t, c t) -> c (Complex t)
357 fromComplex :: (Container c t, RealElement t) => c (Complex t) -> (c t, c t)
358
359
360instance Convert Double where
361 real = id
362 complex = comp'
363 single = single'
364 double = id
365 toComplex = toComplex'
366 fromComplex = fromComplex'
367
368instance Convert Float where
369 real = id
370 complex = comp'
371 single = id
372 double = double'
373 toComplex = toComplex'
374 fromComplex = fromComplex'
375
376instance Convert (Complex Double) where
377 real = comp'
378 complex = id
379 single = single'
380 double = id
381 toComplex = toComplex'
382 fromComplex = fromComplex'
383
384instance Convert (Complex Float) where
385 real = comp'
386 complex = id
387 single = id
388 double = double'
389 toComplex = toComplex'
390 fromComplex = fromComplex'
391
392-------------------------------------------------------------------
393
394type family RealOf x
395
396type instance RealOf Double = Double
397type instance RealOf (Complex Double) = Double
398 88
399type instance RealOf Float = Float 89{- | Creates a real vector containing a range of values:
400type instance RealOf (Complex Float) = Float
401 90
402type family ComplexOf x 91@\> linspace 5 (-3,7)
925 |> [-3.0,-0.5,2.0,4.5,7.0]@
403 93
404type instance ComplexOf Double = Complex Double 94Logarithmic spacing can be defined as follows:
405type instance ComplexOf (Complex Double) = Complex Double
406 95
407type instance ComplexOf Float = Complex Float 96@logspace n (a,b) = 10 ** linspace n (a,b)@
408type instance ComplexOf (Complex Float) = Complex Float
409
410type family SingleOf x
411
412type instance SingleOf Double = Float
413type instance SingleOf Float = Float
414
415type instance SingleOf (Complex a) = Complex (SingleOf a)
416
417type family DoubleOf x
418
419type instance DoubleOf Double = Double
420type instance DoubleOf Float = Double
421
422type instance DoubleOf (Complex a) = Complex (DoubleOf a)
423
424type family ElementOf c
425
426type instance ElementOf (Vector a) = a
427type instance ElementOf (Matrix a) = a
428
429------------------------------------------------------------
430
431conjugateAux fun x = unsafePerformIO $ do
432 v <- createVector (dim x)
433 app2 fun vec x vec v "conjugateAux"
434 return v
435
436conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
437conjugateQ = conjugateAux c_conjugateQ
438foreign import ccall "conjugateQ" c_conjugateQ :: TQVQV
439
440conjugateC :: Vector (Complex Double) -> Vector (Complex Double)
441conjugateC = conjugateAux c_conjugateC
442foreign import ccall "conjugateC" c_conjugateC :: TCVCV
443
444----------------------------------------------------
445
446{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}
447
448-- -- | @x .* a = scale x a@
449-- (.*) :: (Linear c a) => a -> c a -> c a
450infixl 7 .*
451a .* x = scale a x
452
453----------------------------------------------------
454
455{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}
456
457-- -- | @a *\/ x = scale (recip x) a@
458-- (*/) :: (Linear c a) => c a -> a -> c a
459infixl 7 */
460v */ x = scale (recip x) v
461
462
463------------------------------------------------
464
465{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}
466{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}
467
468class Joinable a b where
469 joinH :: Element t => a t -> b t -> Matrix t
470 joinV :: Element t => a t -> b t -> Matrix t
471
472instance Joinable Matrix Matrix where
473 joinH m1 m2 = fromBlocks [[m1,m2]]
474 joinV m1 m2 = fromBlocks [[m1],[m2]]
475
476instance Joinable Matrix Vector where
477 joinH m v = joinH m (asColumn v)
478 joinV m v = joinV m (asRow v)
479
480instance Joinable Vector Matrix where
481 joinH v m = joinH (asColumn v) m
482 joinV v m = joinV (asRow v) m
483
484infixl 4 <|>
485infixl 3 <->
486
487{-- - | Horizontal concatenation of matrices and vectors:
488
489@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]
490(6><4)
491 [ 1.0, 0.0, 0.0, 1.0
492 , 0.0, 1.0, 0.0, 2.0
493 , 0.0, 0.0, 1.0, 3.0
494 , 3.0, 0.0, 0.0, 4.0
495 , 0.0, 3.0, 0.0, 5.0
496 , 0.0, 0.0, 3.0, 6.0 ]@
497-} 97-}
498-- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t 98linspace :: (Enum e, Container Vector e) => Int -> (e, e) -> Vector e
499a <|> b = joinH a b 99linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1]
500 100 where s = (b-a)/fromIntegral (n-1)
501-- -- | Vertical concatenation of matrices and vectors.
502-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
503a <-> b = joinV a b
504
505-------------------------------------------------------------------
506
507{-# DEPRECATED vectorMin "use minElement" #-}
508vectorMin :: (Container Vector t, Element t) => Vector t -> t
509vectorMin = minElement
510
511{-# DEPRECATED vectorMax "use maxElement" #-}
512vectorMax :: (Container Vector t, Element t) => Vector t -> t
513vectorMax = maxElement
514
515
516{-# DEPRECATED vectorMaxIndex "use minIndex" #-}
517vectorMaxIndex :: Vector Double -> Int
518vectorMaxIndex = round . toScalarR MaxIdx
519
520{-# DEPRECATED vectorMinIndex "use maxIndex" #-}
521vectorMinIndex :: Vector Double -> Int
522vectorMinIndex = round . toScalarR MinIdx
523
524-----------------------------------------------------
525
526class Build f where
527 build' :: BoundsOf f -> f -> ContainerOf f
528
529type family BoundsOf x
530
531type instance BoundsOf (a->a) = Int
532type instance BoundsOf (a->a->a) = (Int,Int)
533
534type family ContainerOf x
535
536type instance ContainerOf (a->a) = Vector a
537type instance ContainerOf (a->a->a) = Matrix a
538 101
539instance (Element a, Num a) => Build (a->a) where 102-- | Dot product: @u \<.\> v = dot u v@
540 build' = buildV 103(<.>) :: Product t => Vector t -> Vector t -> t
104infixl 7 <.>
105(<.>) = dot
541 106
542instance (Element a, Num a) => Build (a->a->a) where
543 build' = buildM
544 107
545buildM (rc,cc) f = fromLists [ [f r c | c <- cs] | r <- rs ]
546 where rs = map fromIntegral [0 .. (rc-1)]
547 cs = map fromIntegral [0 .. (cc-1)]
548 108
549buildV n f = fromList [f k | k <- ks] 109--------------------------------------------------------
550 where ks = map fromIntegral [0 .. (n-1)]
551 110
552---------------------------------------------------- 111class Mul a b c | a b -> c where
553-- experimental 112 infixl 7 <>
113 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
114 (<>) :: Product t => a t -> b t -> c t
554 115
555class Konst s where 116instance Mul Matrix Matrix Matrix where
556 konst' :: Element e => e -> s -> ContainerOf' s e 117 (<>) = mXm
557 118
558type family ContainerOf' x y 119instance Mul Matrix Vector Vector where
120 (<>) m v = flatten $ m <> (asColumn v)
559 121
560type instance ContainerOf' Int a = Vector a 122instance Mul Vector Matrix Vector where
561type instance ContainerOf' (Int,Int) a = Matrix a 123 (<>) v m = flatten $ (asRow v) <> m
562 124
563instance Konst Int where 125--------------------------------------------------------
564 konst' = constantD
565 126
566instance Konst (Int,Int) where 127-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).
567 konst' k (r,c) = reshape c $ konst' k (r*c) 128(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a
129infixl 7 <\>
130m <\> v = flatten (linearSolveSVD m (reshape 1 v))
diff --git a/lib/Numeric/ContainerBoot.hs b/lib/Numeric/ContainerBoot.hs
new file mode 100644
index 0000000..1284639
--- /dev/null
+++ b/lib/Numeric/ContainerBoot.hs
@@ -0,0 +1,575 @@
1{-# LANGUAGE TypeFamilies #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE UndecidableInstances #-}
7
8-----------------------------------------------------------------------------
9-- |
10-- Module : Numeric.ContainerBoot
11-- Copyright : (c) Alberto Ruiz 2010
12-- License : GPL-style
13--
14-- Maintainer : Alberto Ruiz <aruiz@um.es>
15-- Stability : provisional
16-- Portability : portable
17--
18-- Module to avoid cyclyc dependencies.
19--
20-----------------------------------------------------------------------------
21
22module Numeric.ContainerBoot (
23 -- * Basic functions
24 ident, diag, ctrans,
25 -- * Generic operations
26 Container(..),
27 -- * Matrix product and related functions
28 Product(..),
29 mXm,mXv,vXm,
30 outer, kronecker,
31 -- * Element conversion
32 Convert(..),
33 Complexable(),
34 RealElement(),
35
36 RealOf, ComplexOf, SingleOf, DoubleOf,
37
38 IndexOf,
39 module Data.Complex,
40 -- * Experimental
41 build', konst',
42 -- * Deprecated
43 (.*),(*/),(<|>),(<->),
44 vectorMax,vectorMin,
45 vectorMaxIndex, vectorMinIndex
46) where
47
48import Data.Packed
49import Numeric.Conversion
50import Data.Packed.Internal
51import Numeric.GSL.Vector
52
53import Data.Complex
54import Control.Monad(ap)
55
56import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
57
58import System.IO.Unsafe
59
60-------------------------------------------------------------------
61
62type family IndexOf c
63
64type instance IndexOf Vector = Int
65type instance IndexOf Matrix = (Int,Int)
66
67type family ArgOf c a
68
69type instance ArgOf Vector a = a -> a
70type instance ArgOf Matrix a = a -> a -> a
71
72-------------------------------------------------------------------
73
74-- | Basic element-by-element functions for numeric containers
75class (Complexable c, Fractional e, Element e) => Container c e where
76 -- | create a structure with a single element
77 scalar :: e -> c e
78 -- | complex conjugate
79 conj :: c e -> c e
80 scale :: e -> c e -> c e
81 -- | scale the element by element reciprocal of the object:
82 --
83 -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@
84 scaleRecip :: e -> c e -> c e
85 addConstant :: e -> c e -> c e
86 add :: c e -> c e -> c e
87 sub :: c e -> c e -> c e
88 -- | element by element multiplication
89 mul :: c e -> c e -> c e
90 -- | element by element division
91 divide :: c e -> c e -> c e
92 equal :: c e -> c e -> Bool
93 --
94 -- | cannot implement instance Functor because of Element class constraint
95 cmap :: (Element a, Element b) => (a -> b) -> c a -> c b
96 -- | constant structure of given size
97 konst :: e -> IndexOf c -> c e
98 -- | create a structure using a function
99 --
100 -- Hilbert matrix of order N:
101 --
102 -- @hilb n = build (n,n) (\\i j -> 1/(i+j+1))@
103 build :: IndexOf c -> (ArgOf c e) -> c e
104 --build :: BoundsOf f -> f -> (ContainerOf f) e
105 --
106 -- | indexing function
107 atIndex :: c e -> IndexOf c -> e
108 -- | index of min element
109 minIndex :: c e -> IndexOf c
110 -- | index of max element
111 maxIndex :: c e -> IndexOf c
112 -- | value of min element
113 minElement :: c e -> e
114 -- | value of max element
115 maxElement :: c e -> e
116 -- the C functions sumX/prodX are twice as fast as using foldVector
117 -- | the sum of elements (faster than using @fold@)
118 sumElements :: c e -> e
119 -- | the product of elements (faster than using @fold@)
120 prodElements :: c e -> e
121
122--------------------------------------------------------------------------
123
124instance Container Vector Float where
125 scale = vectorMapValF Scale
126 scaleRecip = vectorMapValF Recip
127 addConstant = vectorMapValF AddConstant
128 add = vectorZipF Add
129 sub = vectorZipF Sub
130 mul = vectorZipF Mul
131 divide = vectorZipF Div
132 equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
133 scalar x = fromList [x]
134 konst = constantD
135 build = buildV
136 conj = id
137 cmap = mapVector
138 atIndex = (@>)
139 minIndex = round . toScalarF MinIdx
140 maxIndex = round . toScalarF MaxIdx
141 minElement = toScalarF Min
142 maxElement = toScalarF Max
143 sumElements = sumF
144 prodElements = prodF
145
146instance Container Vector Double where
147 scale = vectorMapValR Scale
148 scaleRecip = vectorMapValR Recip
149 addConstant = vectorMapValR AddConstant
150 add = vectorZipR Add
151 sub = vectorZipR Sub
152 mul = vectorZipR Mul
153 divide = vectorZipR Div
154 equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
155 scalar x = fromList [x]
156 konst = constantD
157 build = buildV
158 conj = id
159 cmap = mapVector
160 atIndex = (@>)
161 minIndex = round . toScalarR MinIdx
162 maxIndex = round . toScalarR MaxIdx
163 minElement = toScalarR Min
164 maxElement = toScalarR Max
165 sumElements = sumR
166 prodElements = prodR
167
168instance Container Vector (Complex Double) where
169 scale = vectorMapValC Scale
170 scaleRecip = vectorMapValC Recip
171 addConstant = vectorMapValC AddConstant
172 add = vectorZipC Add
173 sub = vectorZipC Sub
174 mul = vectorZipC Mul
175 divide = vectorZipC Div
176 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
177 scalar x = fromList [x]
178 konst = constantD
179 build = buildV
180 conj = conjugateC
181 cmap = mapVector
182 atIndex = (@>)
183 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
184 maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
185 minElement = ap (@>) minIndex
186 maxElement = ap (@>) maxIndex
187 sumElements = sumC
188 prodElements = prodC
189
190instance Container Vector (Complex Float) where
191 scale = vectorMapValQ Scale
192 scaleRecip = vectorMapValQ Recip
193 addConstant = vectorMapValQ AddConstant
194 add = vectorZipQ Add
195 sub = vectorZipQ Sub
196 mul = vectorZipQ Mul
197 divide = vectorZipQ Div
198 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
199 scalar x = fromList [x]
200 konst = constantD
201 build = buildV
202 conj = conjugateQ
203 cmap = mapVector
204 atIndex = (@>)
205 minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
206 maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate)
207 minElement = ap (@>) minIndex
208 maxElement = ap (@>) maxIndex
209 sumElements = sumQ
210 prodElements = prodQ
211
212---------------------------------------------------------------
213
214instance (Container Vector a) => Container Matrix a where
215 scale x = liftMatrix (scale x)
216 scaleRecip x = liftMatrix (scaleRecip x)
217 addConstant x = liftMatrix (addConstant x)
218 add = liftMatrix2 add
219 sub = liftMatrix2 sub
220 mul = liftMatrix2 mul
221 divide = liftMatrix2 divide
222 equal a b = cols a == cols b && flatten a `equal` flatten b
223 scalar x = (1><1) [x]
224 konst v (r,c) = reshape c (konst v (r*c))
225 build = buildM
226 conj = liftMatrix conj
227 cmap f = liftMatrix (mapVector f)
228 atIndex = (@@>)
229 minIndex m = let (r,c) = (rows m,cols m)
230 i = (minIndex $ flatten m)
231 in (i `div` c,i `mod` c)
232 maxIndex m = let (r,c) = (rows m,cols m)
233 i = (maxIndex $ flatten m)
234 in (i `div` c,i `mod` c)
235 minElement = ap (@@>) minIndex
236 maxElement = ap (@@>) maxIndex
237 sumElements = sumElements . flatten
238 prodElements = prodElements . flatten
239
240----------------------------------------------------
241
242-- | Matrix product and related functions
243class Element e => Product e where
244 -- | matrix product
245 multiply :: Matrix e -> Matrix e -> Matrix e
246 -- | dot (inner) product
247 dot :: Vector e -> Vector e -> e
248 -- | sum of absolute value of elements (differs in complex case from @norm1@)
249 absSum :: Vector e -> RealOf e
250 -- | sum of absolute value of elements
251 norm1 :: Vector e -> RealOf e
252 -- | euclidean norm
253 norm2 :: Vector e -> RealOf e
254 -- | element of maximum magnitude
255 normInf :: Vector e -> RealOf e
256
257instance Product Float where
258 norm2 = toScalarF Norm2
259 absSum = toScalarF AbsSum
260 dot = dotF
261 norm1 = toScalarF AbsSum
262 normInf = maxElement . vectorMapF Abs
263 multiply = multiplyF
264
265instance Product Double where
266 norm2 = toScalarR Norm2
267 absSum = toScalarR AbsSum
268 dot = dotR
269 norm1 = toScalarR AbsSum
270 normInf = maxElement . vectorMapR Abs
271 multiply = multiplyR
272
273instance Product (Complex Float) where
274 norm2 = toScalarQ Norm2
275 absSum = toScalarQ AbsSum
276 dot = dotQ
277 norm1 = sumElements . fst . fromComplex . vectorMapQ Abs
278 normInf = maxElement . fst . fromComplex . vectorMapQ Abs
279 multiply = multiplyQ
280
281instance Product (Complex Double) where
282 norm2 = toScalarC Norm2
283 absSum = toScalarC AbsSum
284 dot = dotC
285 norm1 = sumElements . fst . fromComplex . vectorMapC Abs
286 normInf = maxElement . fst . fromComplex . vectorMapC Abs
287 multiply = multiplyC
288
289----------------------------------------------------------
290
291-- synonym for matrix product
292mXm :: Product t => Matrix t -> Matrix t -> Matrix t
293mXm = multiply
294
295-- matrix - vector product
296mXv :: Product t => Matrix t -> Vector t -> Vector t
297mXv m v = flatten $ m `mXm` (asColumn v)
298
299-- vector - matrix product
300vXm :: Product t => Vector t -> Matrix t -> Vector t
301vXm v m = flatten $ (asRow v) `mXm` m
302
303{- | Outer product of two vectors.
304
305@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3]
306(3><3)
307 [ 5.0, 2.0, 3.0
308 , 10.0, 4.0, 6.0
309 , 15.0, 6.0, 9.0 ]@
310-}
311outer :: (Product t) => Vector t -> Vector t -> Matrix t
312outer u v = asColumn u `multiply` asRow v
313
314{- | Kronecker product of two matrices.
315
316@m1=(2><3)
317 [ 1.0, 2.0, 0.0
318 , 0.0, -1.0, 3.0 ]
319m2=(4><3)
320 [ 1.0, 2.0, 3.0
321 , 4.0, 5.0, 6.0
322 , 7.0, 8.0, 9.0
323 , 10.0, 11.0, 12.0 ]@
324
325@\> kronecker m1 m2
326(8><9)
327 [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0
328 , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0
329 , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0
330 , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0
331 , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0
332 , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0
333 , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0
334 , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@
335-}
336kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
337kronecker a b = fromBlocks
338 . splitEvery (cols a)
339 . map (reshape (cols b))
340 . toRows
341 $ flatten a `outer` flatten b
342
343-------------------------------------------------------------------
344
345
346class Convert t where
347 real :: Container c t => c (RealOf t) -> c t
348 complex :: Container c t => c t -> c (ComplexOf t)
349 single :: Container c t => c t -> c (SingleOf t)
350 double :: Container c t => c t -> c (DoubleOf t)
351 toComplex :: (Container c t, RealElement t) => (c t, c t) -> c (Complex t)
352 fromComplex :: (Container c t, RealElement t) => c (Complex t) -> (c t, c t)
353
354
355instance Convert Double where
356 real = id
357 complex = comp'
358 single = single'
359 double = id
360 toComplex = toComplex'
361 fromComplex = fromComplex'
362
363instance Convert Float where
364 real = id
365 complex = comp'
366 single = id
367 double = double'
368 toComplex = toComplex'
369 fromComplex = fromComplex'
370
371instance Convert (Complex Double) where
372 real = comp'
373 complex = id
374 single = single'
375 double = id
376 toComplex = toComplex'
377 fromComplex = fromComplex'
378
379instance Convert (Complex Float) where
380 real = comp'
381 complex = id
382 single = id
383 double = double'
384 toComplex = toComplex'
385 fromComplex = fromComplex'
386
387-------------------------------------------------------------------
388
389type family RealOf x
390
391type instance RealOf Double = Double
392type instance RealOf (Complex Double) = Double
393
394type instance RealOf Float = Float
395type instance RealOf (Complex Float) = Float
396
397type family ComplexOf x
398
399type instance ComplexOf Double = Complex Double
400type instance ComplexOf (Complex Double) = Complex Double
401
402type instance ComplexOf Float = Complex Float
403type instance ComplexOf (Complex Float) = Complex Float
404
405type family SingleOf x
406
407type instance SingleOf Double = Float
408type instance SingleOf Float = Float
409
410type instance SingleOf (Complex a) = Complex (SingleOf a)
411
412type family DoubleOf x
413
414type instance DoubleOf Double = Double
415type instance DoubleOf Float = Double
416
417type instance DoubleOf (Complex a) = Complex (DoubleOf a)
418
419type family ElementOf c
420
421type instance ElementOf (Vector a) = a
422type instance ElementOf (Matrix a) = a
423
424------------------------------------------------------------
425
426conjugateAux fun x = unsafePerformIO $ do
427 v <- createVector (dim x)
428 app2 fun vec x vec v "conjugateAux"
429 return v
430
431conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
432conjugateQ = conjugateAux c_conjugateQ
433foreign import ccall "conjugateQ" c_conjugateQ :: TQVQV
434
435conjugateC :: Vector (Complex Double) -> Vector (Complex Double)
436conjugateC = conjugateAux c_conjugateC
437foreign import ccall "conjugateC" c_conjugateC :: TCVCV
438
439----------------------------------------------------
440
441{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}
442
443-- -- | @x .* a = scale x a@
444-- (.*) :: (Linear c a) => a -> c a -> c a
445infixl 7 .*
446a .* x = scale a x
447
448----------------------------------------------------
449
450{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}
451
452-- -- | @a *\/ x = scale (recip x) a@
453-- (*/) :: (Linear c a) => c a -> a -> c a
454infixl 7 */
455v */ x = scale (recip x) v
456
457
458------------------------------------------------
459
460{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}
461{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}
462
463class Joinable a b where
464 joinH :: Element t => a t -> b t -> Matrix t
465 joinV :: Element t => a t -> b t -> Matrix t
466
467instance Joinable Matrix Matrix where
468 joinH m1 m2 = fromBlocks [[m1,m2]]
469 joinV m1 m2 = fromBlocks [[m1],[m2]]
470
471instance Joinable Matrix Vector where
472 joinH m v = joinH m (asColumn v)
473 joinV m v = joinV m (asRow v)
474
475instance Joinable Vector Matrix where
476 joinH v m = joinH (asColumn v) m
477 joinV v m = joinV (asRow v) m
478
479infixl 4 <|>
480infixl 3 <->
481
482{-- - | Horizontal concatenation of matrices and vectors:
483
484@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]
485(6><4)
486 [ 1.0, 0.0, 0.0, 1.0
487 , 0.0, 1.0, 0.0, 2.0
488 , 0.0, 0.0, 1.0, 3.0
489 , 3.0, 0.0, 0.0, 4.0
490 , 0.0, 3.0, 0.0, 5.0
491 , 0.0, 0.0, 3.0, 6.0 ]@
492-}
493-- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
494a <|> b = joinH a b
495
496-- -- | Vertical concatenation of matrices and vectors.
497-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
498a <-> b = joinV a b
499
500-------------------------------------------------------------------
501
502{-# DEPRECATED vectorMin "use minElement" #-}
503vectorMin :: (Container Vector t, Element t) => Vector t -> t
504vectorMin = minElement
505
506{-# DEPRECATED vectorMax "use maxElement" #-}
507vectorMax :: (Container Vector t, Element t) => Vector t -> t
508vectorMax = maxElement
509
510
511{-# DEPRECATED vectorMaxIndex "use minIndex" #-}
512vectorMaxIndex :: Vector Double -> Int
513vectorMaxIndex = round . toScalarR MaxIdx
514
515{-# DEPRECATED vectorMinIndex "use maxIndex" #-}
516vectorMinIndex :: Vector Double -> Int
517vectorMinIndex = round . toScalarR MinIdx
518
519-----------------------------------------------------
520
521class Build f where
522 build' :: BoundsOf f -> f -> ContainerOf f
523
524type family BoundsOf x
525
526type instance BoundsOf (a->a) = Int
527type instance BoundsOf (a->a->a) = (Int,Int)
528
529type family ContainerOf x
530
531type instance ContainerOf (a->a) = Vector a
532type instance ContainerOf (a->a->a) = Matrix a
533
534instance (Element a, Num a) => Build (a->a) where
535 build' = buildV
536
537instance (Element a, Num a) => Build (a->a->a) where
538 build' = buildM
539
540buildM (rc,cc) f = fromLists [ [f r c | c <- cs] | r <- rs ]
541 where rs = map fromIntegral [0 .. (rc-1)]
542 cs = map fromIntegral [0 .. (cc-1)]
543
544buildV n f = fromList [f k | k <- ks]
545 where ks = map fromIntegral [0 .. (n-1)]
546
547----------------------------------------------------
548-- experimental
549
550class Konst s where
551 konst' :: Element e => e -> s -> ContainerOf' s e
552
553type family ContainerOf' x y
554
555type instance ContainerOf' Int a = Vector a
556type instance ContainerOf' (Int,Int) a = Matrix a
557
558instance Konst Int where
559 konst' = constantD
560
561instance Konst (Int,Int) where
562 konst' k (r,c) = reshape c $ konst' k (r*c)
563
564--------------------------------------------------------
565-- | conjugate transpose
566ctrans :: (Container Vector e, Element e) => Matrix e -> Matrix e
567ctrans = liftMatrix conj . trans
568
569-- | Creates a square matrix with a given diagonal.
570diag :: (Num a, Element a) => Vector a -> Matrix a
571diag v = diagRect 0 v n n where n = dim v
572
573-- | creates the identity matrix of given dimension
574ident :: (Num a, Element a) => Int -> Matrix a
575ident n = diag (constantD 1 n)
diff --git a/lib/Numeric/LinearAlgebra.hs b/lib/Numeric/LinearAlgebra.hs
index edfd87e..6b1a0bd 100644
--- a/lib/Numeric/LinearAlgebra.hs
+++ b/lib/Numeric/LinearAlgebra.hs
@@ -1,7 +1,7 @@
1----------------------------------------------------------------------------- 1-----------------------------------------------------------------------------
2{- | 2{- |
3Module : Numeric.LinearAlgebra 3Module : Numeric.LinearAlgebra
4Copyright : (c) Alberto Ruiz 2006-9 4Copyright : (c) Alberto Ruiz 2006-10
5License : GPL-style 5License : GPL-style
6 6
7Maintainer : Alberto Ruiz (aruiz at um dot es) 7Maintainer : Alberto Ruiz (aruiz at um dot es)
@@ -10,16 +10,19 @@ Portability : uses ffi
10 10
11This module reexports all normally required functions for Linear Algebra applications. 11This module reexports all normally required functions for Linear Algebra applications.
12 12
13It also provides instances of standard classes 'Show', 'Read', 'Eq',
14'Num', 'Fractional', and 'Floating' for 'Vector' and 'Matrix'.
15In arithmetic operations one-component vectors and matrices automatically
16expand to match the dimensions of the other operand.
17
13-} 18-}
14----------------------------------------------------------------------------- 19-----------------------------------------------------------------------------
15module Numeric.LinearAlgebra ( 20module Numeric.LinearAlgebra (
16 module Numeric.LinearAlgebra.Algorithms, 21 module Numeric.Container,
17-- module Numeric.LinearAlgebra.Interface, 22 module Numeric.LinearAlgebra.Algorithms
18 module Numeric.Matrix,
19 module Data.Packed.Random
20) where 23) where
21 24
25import Numeric.Container
22import Numeric.LinearAlgebra.Algorithms 26import Numeric.LinearAlgebra.Algorithms
23--import Numeric.LinearAlgebra.Interface 27import Numeric.Matrix()
24import Numeric.Matrix 28import Numeric.Vector() \ No newline at end of file
25import Data.Packed.Random
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index c49bec7..f4f8bca 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -13,7 +13,7 @@ Maintainer : Alberto Ruiz (aruiz at um dot es)
13Stability : provisional 13Stability : provisional
14Portability : uses ffi 14Portability : uses ffi
15 15
16Generic interface for the most common functions. Using it we can write higher level algorithms and testing properties for both real and complex matrices. 16High level generic interface to common matrix computations.
17 17
18Specific functions for particular base types can also be explicitly 18Specific functions for particular base types can also be explicitly
19imported from "Numeric.LinearAlgebra.LAPACK". 19imported from "Numeric.LinearAlgebra.LAPACK".
@@ -63,6 +63,7 @@ module Numeric.LinearAlgebra.Algorithms (
63 nullspaceSVD, 63 nullspaceSVD,
64-- * Norms 64-- * Norms
65 Normed(..), NormType(..), 65 Normed(..), NormType(..),
66 relativeError,
66-- * Misc 67-- * Misc
67 eps, peps, i, 68 eps, peps, i,
68-- * Util 69-- * Util
@@ -79,10 +80,10 @@ import Data.Packed.Matrix
79import Numeric.LinearAlgebra.LAPACK as LAPACK 80import Numeric.LinearAlgebra.LAPACK as LAPACK
80import Data.List(foldl1') 81import Data.List(foldl1')
81import Data.Array 82import Data.Array
82import Numeric.Container hiding ((.*),(*/)) 83import Numeric.ContainerBoot hiding ((.*),(*/))
83import Numeric.MatrixBoot
84 84
85{- | Auxiliary typeclass used to define generic linear algebra computations for both real and complex matrices. Only double precision is supported in this version (we can 85
86{- | Class used to define generic linear algebra computations for both real and complex matrices. Only double precision is supported in this version (we can
86transform single precision objects using 'single' and 'double'). 87transform single precision objects using 'single' and 'double').
87 88
88-} 89-}
@@ -691,3 +692,9 @@ instance Normed Matrix (Complex Float) where
691 pnorm PNorm2 = realToFrac . (@>0) . singularValues . double 692 pnorm PNorm2 = realToFrac . (@>0) . singularValues . double
692 pnorm Infinity = pnorm PNorm1 . trans 693 pnorm Infinity = pnorm PNorm1 . trans
693 pnorm Frobenius = pnorm PNorm2 . flatten 694 pnorm Frobenius = pnorm PNorm2 . flatten
695
696-- | Approximate number of common digits in the maximum element.
697relativeError :: (Normed c t, Container c t) => c t -> c t -> Int
698relativeError x y = dig (norm (x `sub` y) / norm x)
699 where norm = pnorm Infinity
700 dig r = round $ -logBase 10 (realToFrac r :: Double)
diff --git a/lib/Numeric/Matrix.hs b/lib/Numeric/Matrix.hs
index 7bb79f3..bab5a27 100644
--- a/lib/Numeric/Matrix.hs
+++ b/lib/Numeric/Matrix.hs
@@ -15,40 +15,20 @@
15-- Stability : provisional 15-- Stability : provisional
16-- Portability : portable 16-- Portability : portable
17-- 17--
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.
21--
22-- Provides instances of standard classes 'Show', 'Read', 'Eq', 18-- Provides instances of standard classes 'Show', 'Read', 'Eq',
23-- 'Num', 'Fractional', and 'Floating' for 'Matrix'. 19-- 'Num', 'Fractional', and 'Floating' for 'Matrix'.
24-- 20--
21-- In arithmetic operations one-component
22-- vectors and matrices automatically expand to match the dimensions of the other operand.
23
25----------------------------------------------------------------------------- 24-----------------------------------------------------------------------------
26 25
27module Numeric.Matrix ( 26module Numeric.Matrix (
28 -- * Basic functions
29 module Data.Packed.Matrix,
30 module Numeric.Vector,
31 -- * Matrix creation
32 diag, ident,
33 -- * matrix operations
34 ctrans,
35 optimiseMult,
36 -- * Operators
37 (<>), (<\>),
38 -- * IO
39 dispf, disps, dispcf, vecdisp, latexFormat, format,
40 loadMatrix, saveMatrix, fromFile, fileDimensions,
41 readMatrix
42 ) where 27 ) where
43 28
44------------------------------------------------------------------- 29-------------------------------------------------------------------
45 30
46import Data.Packed.Matrix 31import Numeric.Container
47import Numeric.Vector
48import Numeric.Chain
49import Numeric.MatrixBoot
50import Numeric.IO
51import Numeric.LinearAlgebra.Algorithms
52 32
53------------------------------------------------------------------- 33-------------------------------------------------------------------
54 34
@@ -90,26 +70,3 @@ instance (Floating a, Container Vector a, Floating (Vector a), Fractional (Matri
90 (**) = liftMatrix2Auto (**) 70 (**) = liftMatrix2Auto (**)
91 sqrt = liftMatrix sqrt 71 sqrt = liftMatrix sqrt
92 pi = (1><1) [pi] 72 pi = (1><1) [pi]
93
94--------------------------------------------------------
95
96class Mul a b c | a b -> c where
97 infixl 7 <>
98 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
99 (<>) :: Product t => a t -> b t -> c t
100
101instance Mul Matrix Matrix Matrix where
102 (<>) = mXm
103
104instance Mul Matrix Vector Vector where
105 (<>) m v = flatten $ m <> (asColumn v)
106
107instance Mul Vector Matrix Vector where
108 (<>) v m = flatten $ (asRow v) <> m
109
110--------------------------------------------------------
111
112-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).
113(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a
114infixl 7 <\>
115m <\> v = flatten (linearSolveSVD m (reshape 1 v))
diff --git a/lib/Numeric/MatrixBoot.hs b/lib/Numeric/MatrixBoot.hs
deleted file mode 100644
index 13560b1..0000000
--- a/lib/Numeric/MatrixBoot.hs
+++ /dev/null
@@ -1,41 +0,0 @@
1{-# LANGUAGE FlexibleContexts #-}
2-----------------------------------------------------------------------------
3-- |
4-- Module : Numeric.MatrixBoot
5-- Copyright : (c) Alberto Ruiz 2010
6-- License : GPL-style
7--
8-- Maintainer : Alberto Ruiz <aruiz@um.es>
9-- Stability : provisional
10-- Portability : portable
11--
12-- Module to avoid cyclic dependency
13--
14-----------------------------------------------------------------------------
15
16module Numeric.MatrixBoot (
17 ctrans, diag, ident,
18 ) where
19
20import Data.Packed.Vector
21import Data.Packed.Matrix
22import Data.Packed.Internal.Matrix
23import Numeric.Container
24
25-------------------------------------------------------------------
26
27-- | conjugate transpose
28ctrans :: (Container Vector e, Element e) => Matrix e -> Matrix e
29ctrans = liftMatrix conj . trans
30
31-- | Creates a square matrix with a given diagonal.
32diag :: (Num a, Element a) => Vector a -> Matrix a
33diag v = diagRect 0 v n n where n = dim v
34
35-- | creates the identity matrix of given dimension
36ident :: (Num a, Element a) => Int -> Matrix a
37ident n = diag (constantD 1 n)
38
39----------------------------------------------------
40
41
diff --git a/lib/Numeric/Vector.hs b/lib/Numeric/Vector.hs
index 9acd6f0..3fcd880 100644
--- a/lib/Numeric/Vector.hs
+++ b/lib/Numeric/Vector.hs
@@ -13,30 +13,16 @@
13-- Stability : provisional 13-- Stability : provisional
14-- Portability : portable 14-- Portability : portable
15-- 15--
16-- Numeric instances and functions for 'Vector'.
17--
18-- Provides instances of standard classes 'Show', 'Read', 'Eq', 16-- Provides instances of standard classes 'Show', 'Read', 'Eq',
19-- 'Num', 'Fractional', and 'Floating' for 'Vector'. 17-- 'Num', 'Fractional', and 'Floating' for 'Vector'.
20-- 18--
21----------------------------------------------------------------------------- 19-----------------------------------------------------------------------------
22 20
23module Numeric.Vector ( 21module Numeric.Vector (
24 -- * Basic vector functions
25 module Data.Packed.Vector,
26 module Numeric.Container,
27 -- * Vector creation
28 constant, linspace,
29 -- * Operators
30 (<.>),
31 -- * IO
32 fscanfVector, fprintfVector, freadVector, fwriteVector
33 ) where 22 ) where
34 23
35import Data.Packed.Vector
36import Data.Packed.Internal.Matrix
37import Numeric.GSL.Vector 24import Numeric.GSL.Vector
38import Numeric.Container 25import Numeric.Container
39import Numeric.IO
40 26
41------------------------------------------------------------------- 27-------------------------------------------------------------------
42 28
@@ -73,34 +59,6 @@ instance (Element a, Read a) => Read (Vector a) where
73breakAt c l = (a++[c],tail b) where 59breakAt c l = (a++[c],tail b) where
74 (a,b) = break (==c) l 60 (a,b) = break (==c) l
75 61
76------------------------------------------------------------------
77
78{- | creates a vector with a given number of equal components:
79
80@> constant 2 7
817 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
82-}
83constant :: Element a => a -> Int -> Vector a
84-- constant x n = runSTVector (newVector x n)
85constant = constantD -- about 2x faster
86
87{- | Creates a real vector containing a range of values:
88
89@\> linspace 5 (-3,7)
905 |> [-3.0,-0.5,2.0,4.5,7.0]@
91
92Logarithmic spacing can be defined as follows:
93
94@logspace n (a,b) = 10 ** linspace n (a,b)@
95-}
96linspace :: (Enum e, Container Vector e) => Int -> (e, e) -> Vector e
97linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1]
98 where s = (b-a)/fromIntegral (n-1)
99
100-- | Dot product: @u \<.\> v = dot u v@
101(<.>) :: Product t => Vector t -> Vector t -> t
102infixl 7 <.>
103(<.>) = dot
104 62
105------------------------------------------------------------------ 63------------------------------------------------------------------
106 64
generated by cgit v1.2.3 (git 2.39.1) at 2024-09-15 02:41:19 +0000