summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
Diffstat (limited to 'packages')
-rw-r--r--packages/base/hmatrix.cabal11
-rw-r--r--packages/base/src/Data/Packed/Development.hs3
-rw-r--r--packages/base/src/Data/Packed/IO.hs2
-rw-r--r--packages/base/src/Data/Packed/Internal/Numeric.hs39
-rw-r--r--packages/base/src/Data/Packed/Matrix.hs13
-rw-r--r--packages/base/src/Data/Packed/Numeric.hs288
-rw-r--r--packages/base/src/Data/Packed/Vector.hs16
-rw-r--r--packages/base/src/Numeric/Container.hs273
-rw-r--r--packages/base/src/Numeric/HMatrix.hs158
-rw-r--r--packages/base/src/Numeric/LinearAlgebra.hs149
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Algorithms.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Data.hs18
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Util.hs149
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Util/CG.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs18
-rw-r--r--packages/base/src/Numeric/Sparse.hs6
-rw-r--r--packages/gsl/src/Numeric/GSL/Random.hs4
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests.hs8
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs2
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs2
20 files changed, 691 insertions, 472 deletions
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index a92409b..01e3c26 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -38,18 +38,18 @@ library
38 Data.Packed.ST, 38 Data.Packed.ST,
39 Data.Packed.Development, 39 Data.Packed.Development,
40 40
41 Numeric.LinearAlgebra
41 Numeric.LinearAlgebra.LAPACK 42 Numeric.LinearAlgebra.LAPACK
42 Numeric.LinearAlgebra.Algorithms 43 Numeric.LinearAlgebra.Algorithms
43 Numeric.Container 44 Numeric.Container
44 Numeric.LinearAlgebra.Util 45 Numeric.LinearAlgebra.Util
45 46
46 Numeric.LinearAlgebra 47 Numeric.HMatrix
47 Numeric.LinearAlgebra.Devel 48 Numeric.LinearAlgebra.Devel
48 Numeric.LinearAlgebra.Data 49 Numeric.LinearAlgebra.Data
49 50
50 -- Numeric.LinearAlgebra.Compat 51
51 52
52
53 other-modules: Data.Packed.Internal, 53 other-modules: Data.Packed.Internal,
54 Data.Packed.Internal.Common, 54 Data.Packed.Internal.Common,
55 Data.Packed.Internal.Signatures, 55 Data.Packed.Internal.Signatures,
@@ -61,6 +61,7 @@ library
61 Numeric.Vector 61 Numeric.Vector
62 Numeric.Matrix 62 Numeric.Matrix
63 Data.Packed.Internal.Numeric 63 Data.Packed.Internal.Numeric
64 Data.Packed.Numeric
64 Numeric.LinearAlgebra.Util.Convolution 65 Numeric.LinearAlgebra.Util.Convolution
65 Numeric.LinearAlgebra.Util.CG 66 Numeric.LinearAlgebra.Util.CG
66 Numeric.LinearAlgebra.Random 67 Numeric.LinearAlgebra.Random
@@ -69,7 +70,7 @@ library
69 70
70 C-sources: src/C/lapack-aux.c 71 C-sources: src/C/lapack-aux.c
71 src/C/vector-aux.c 72 src/C/vector-aux.c
72 73
73 74
74 extensions: ForeignFunctionInterface, 75 extensions: ForeignFunctionInterface,
75 CPP 76 CPP
diff --git a/packages/base/src/Data/Packed/Development.hs b/packages/base/src/Data/Packed/Development.hs
index 1efedc9..72eb16b 100644
--- a/packages/base/src/Data/Packed/Development.hs
+++ b/packages/base/src/Data/Packed/Development.hs
@@ -25,8 +25,7 @@ module Data.Packed.Development (
25 unsafeFromForeignPtr, 25 unsafeFromForeignPtr,
26 unsafeToForeignPtr, 26 unsafeToForeignPtr,
27 check, (//), 27 check, (//),
28 at', atM', fi, table, 28 at', atM', fi
29 conformMs, conformVs, shSize, splitEvery
30) where 29) where
31 30
32import Data.Packed.Internal 31import Data.Packed.Internal
diff --git a/packages/base/src/Data/Packed/IO.hs b/packages/base/src/Data/Packed/IO.hs
index 94fb30a..f7afa80 100644
--- a/packages/base/src/Data/Packed/IO.hs
+++ b/packages/base/src/Data/Packed/IO.hs
@@ -18,12 +18,12 @@ module Data.Packed.IO (
18) where 18) where
19 19
20import Data.Packed 20import Data.Packed
21import Data.Packed.Development
22import Text.Printf(printf) 21import Text.Printf(printf)
23import Data.List(intersperse) 22import Data.List(intersperse)
24import Data.Complex 23import Data.Complex
25import Numeric.Vectorized(vectorScan,saveMatrix) 24import Numeric.Vectorized(vectorScan,saveMatrix)
26import Control.Applicative((<$>)) 25import Control.Applicative((<$>))
26import Data.Packed.Internal
27 27
28{- | Creates a string from a matrix given a separator and a function to show each entry. Using 28{- | Creates a string from a matrix given a separator and a function to show each entry. Using
29this function the user can easily define any desired display function: 29this function the user can easily define any desired display function:
diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs
index 9cd18df..3c1c1d0 100644
--- a/packages/base/src/Data/Packed/Internal/Numeric.hs
+++ b/packages/base/src/Data/Packed/Internal/Numeric.hs
@@ -49,6 +49,7 @@ import Data.Complex
49import Control.Applicative((<*>)) 49import Control.Applicative((<*>))
50 50
51import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) 51import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
52import Data.Packed.Internal
52 53
53------------------------------------------------------------------- 54-------------------------------------------------------------------
54 55
@@ -65,7 +66,9 @@ type instance ArgOf Matrix a = a -> a -> a
65------------------------------------------------------------------- 66-------------------------------------------------------------------
66 67
67-- | Basic element-by-element functions for numeric containers 68-- | Basic element-by-element functions for numeric containers
68class (Complexable c, Fractional e, Element e) => Container c e where 69class (Complexable c, Fractional e, Element e) => Container c e
70 where
71 size' :: c e -> IndexOf c
69 scalar' :: e -> c e 72 scalar' :: e -> c e
70 conj' :: c e -> c e 73 conj' :: c e -> c e
71 scale' :: e -> c e -> c e 74 scale' :: e -> c e -> c e
@@ -114,7 +117,9 @@ class (Complexable c, Fractional e, Element e) => Container c e where
114 117
115-------------------------------------------------------------------------- 118--------------------------------------------------------------------------
116 119
117instance Container Vector Float where 120instance Container Vector Float
121 where
122 size' = dim
118 scale' = vectorMapValF Scale 123 scale' = vectorMapValF Scale
119 scaleRecip = vectorMapValF Recip 124 scaleRecip = vectorMapValF Recip
120 addConstant = vectorMapValF AddConstant 125 addConstant = vectorMapValF AddConstant
@@ -125,7 +130,7 @@ instance Container Vector Float where
125 equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 130 equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0
126 arctan2' = vectorZipF ATan2 131 arctan2' = vectorZipF ATan2
127 scalar' x = fromList [x] 132 scalar' x = fromList [x]
128 konst' = constant 133 konst' = constantD
129 build' = buildV 134 build' = buildV
130 conj' = id 135 conj' = id
131 cmap' = mapVector 136 cmap' = mapVector
@@ -142,7 +147,9 @@ instance Container Vector Float where
142 accum' = accumV 147 accum' = accumV
143 cond' = condV condF 148 cond' = condV condF
144 149
145instance Container Vector Double where 150instance Container Vector Double
151 where
152 size' = dim
146 scale' = vectorMapValR Scale 153 scale' = vectorMapValR Scale
147 scaleRecip = vectorMapValR Recip 154 scaleRecip = vectorMapValR Recip
148 addConstant = vectorMapValR AddConstant 155 addConstant = vectorMapValR AddConstant
@@ -153,7 +160,7 @@ instance Container Vector Double where
153 equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 160 equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0
154 arctan2' = vectorZipR ATan2 161 arctan2' = vectorZipR ATan2
155 scalar' x = fromList [x] 162 scalar' x = fromList [x]
156 konst' = constant 163 konst' = constantD
157 build' = buildV 164 build' = buildV
158 conj' = id 165 conj' = id
159 cmap' = mapVector 166 cmap' = mapVector
@@ -170,7 +177,9 @@ instance Container Vector Double where
170 accum' = accumV 177 accum' = accumV
171 cond' = condV condD 178 cond' = condV condD
172 179
173instance Container Vector (Complex Double) where 180instance Container Vector (Complex Double)
181 where
182 size' = dim
174 scale' = vectorMapValC Scale 183 scale' = vectorMapValC Scale
175 scaleRecip = vectorMapValC Recip 184 scaleRecip = vectorMapValC Recip
176 addConstant = vectorMapValC AddConstant 185 addConstant = vectorMapValC AddConstant
@@ -181,7 +190,7 @@ instance Container Vector (Complex Double) where
181 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 190 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
182 arctan2' = vectorZipC ATan2 191 arctan2' = vectorZipC ATan2
183 scalar' x = fromList [x] 192 scalar' x = fromList [x]
184 konst' = constant 193 konst' = constantD
185 build' = buildV 194 build' = buildV
186 conj' = conjugateC 195 conj' = conjugateC
187 cmap' = mapVector 196 cmap' = mapVector
@@ -198,7 +207,9 @@ instance Container Vector (Complex Double) where
198 accum' = accumV 207 accum' = accumV
199 cond' = undefined -- cannot match 208 cond' = undefined -- cannot match
200 209
201instance Container Vector (Complex Float) where 210instance Container Vector (Complex Float)
211 where
212 size' = dim
202 scale' = vectorMapValQ Scale 213 scale' = vectorMapValQ Scale
203 scaleRecip = vectorMapValQ Recip 214 scaleRecip = vectorMapValQ Recip
204 addConstant = vectorMapValQ AddConstant 215 addConstant = vectorMapValQ AddConstant
@@ -209,7 +220,7 @@ instance Container Vector (Complex Float) where
209 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 220 equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0
210 arctan2' = vectorZipQ ATan2 221 arctan2' = vectorZipQ ATan2
211 scalar' x = fromList [x] 222 scalar' x = fromList [x]
212 konst' = constant 223 konst' = constantD
213 build' = buildV 224 build' = buildV
214 conj' = conjugateQ 225 conj' = conjugateQ
215 cmap' = mapVector 226 cmap' = mapVector
@@ -228,7 +239,9 @@ instance Container Vector (Complex Float) where
228 239
229--------------------------------------------------------------- 240---------------------------------------------------------------
230 241
231instance (Container Vector a) => Container Matrix a where 242instance (Container Vector a) => Container Matrix a
243 where
244 size' = size
232 scale' x = liftMatrix (scale' x) 245 scale' x = liftMatrix (scale' x)
233 scaleRecip x = liftMatrix (scaleRecip x) 246 scaleRecip x = liftMatrix (scaleRecip x)
234 addConstant x = liftMatrix (addConstant x) 247 addConstant x = liftMatrix (addConstant x)
@@ -637,7 +650,7 @@ diag v = diagRect 0 v n n where n = dim v
637 650
638-- | creates the identity matrix of given dimension 651-- | creates the identity matrix of given dimension
639ident :: (Num a, Element a) => Int -> Matrix a 652ident :: (Num a, Element a) => Int -> Matrix a
640ident n = diag (constant 1 n) 653ident n = diag (constantD 1 n)
641 654
642-------------------------------------------------------- 655--------------------------------------------------------
643 656
@@ -681,8 +694,12 @@ condV f a b l e t = f a' b' l' e' t'
681 694
682class Transposable t 695class Transposable t
683 where 696 where
697 -- | (conjugate) transpose
684 tr :: t -> t 698 tr :: t -> t
685 699
700instance (Container Vector t) => Transposable (Matrix t)
701 where
702 tr = ctrans
686 703
687class Linear t v 704class Linear t v
688 where 705 where
diff --git a/packages/base/src/Data/Packed/Matrix.hs b/packages/base/src/Data/Packed/Matrix.hs
index b3be823..2acb31a 100644
--- a/packages/base/src/Data/Packed/Matrix.hs
+++ b/packages/base/src/Data/Packed/Matrix.hs
@@ -226,15 +226,12 @@ takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (co
226 226
227------------------------------------------------------------ 227------------------------------------------------------------
228 228
229{- | An easy way to create a matrix: 229{- | create a general matrix
230 230
231>>> (2><3)[2,4,7,-3,11,0] 231>>> (2><3) [2, 4, 7+2*𝑖, -3, 11, 0]
232(2><3) 232(2><3)
233 [ 2.0, 4.0, 7.0 233 [ 2.0 :+ 0.0, 4.0 :+ 0.0, 7.0 :+ 2.0
234 , -3.0, 11.0, 0.0 ] 234 , (-3.0) :+ (-0.0), 11.0 :+ 0.0, 0.0 :+ 0.0 ]
235
236This is the format produced by the instances of Show (Matrix a), which
237can also be used for input.
238 235
239The input list is explicitly truncated, so that it can 236The input list is explicitly truncated, so that it can
240safely be used with lists that are too long (like infinite lists). 237safely be used with lists that are too long (like infinite lists).
@@ -244,6 +241,8 @@ safely be used with lists that are too long (like infinite lists).
244 [ 1.0, 2.0, 3.0 241 [ 1.0, 2.0, 3.0
245 , 4.0, 5.0, 6.0 ] 242 , 4.0, 5.0, 6.0 ]
246 243
244This is the format produced by the instances of Show (Matrix a), which
245can also be used for input.
247 246
248-} 247-}
249(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a 248(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a
diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs
new file mode 100644
index 0000000..01cf6c5
--- /dev/null
+++ b/packages/base/src/Data/Packed/Numeric.hs
@@ -0,0 +1,288 @@
1{-# LANGUAGE TypeFamilies #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE UndecidableInstances #-}
7
8-----------------------------------------------------------------------------
9-- |
10-- Module : Data.Packed.Numeric
11-- Copyright : (c) Alberto Ruiz 2010-14
12-- License : BSD3
13-- Maintainer : Alberto Ruiz
14-- Stability : provisional
15--
16-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines.
17--
18-- The 'Container' class is used to define optimized generic functions which work
19-- on 'Vector' and 'Matrix' with real or complex elements.
20--
21-- Some of these functions are also available in the instances of the standard
22-- numeric Haskell classes provided by "Numeric.LinearAlgebra".
23--
24-----------------------------------------------------------------------------
25{-# OPTIONS_HADDOCK hide #-}
26
27module Data.Packed.Numeric (
28 -- * Basic functions
29 module Data.Packed,
30 Konst(..), Build(..),
31 linspace,
32 diag, ident,
33 ctrans,
34 -- * Generic operations
35 Container(..),
36 -- add, mul, sub, divide, equal, scaleRecip, addConstant,
37 scalar, conj, scale, arctan2, cmap,
38 atIndex, minIndex, maxIndex, minElement, maxElement,
39 sumElements, prodElements,
40 step, cond, find, assoc, accum,
41 Transposable(..), Linear(..),
42 -- * Matrix product
43 Product(..), udot, dot, (◇),
44 Mul(..),
45 Contraction(..),(<.>),
46 optimiseMult,
47 mXm,mXv,vXm,LSDiv,(<\>),
48 outer, kronecker,
49 -- * Random numbers
50 RandDist(..),
51 randomVector,
52 gaussianSample,
53 uniformSample,
54 meanCov,
55 -- * Element conversion
56 Convert(..),
57 Complexable(),
58 RealElement(),
59
60 RealOf, ComplexOf, SingleOf, DoubleOf,
61
62 IndexOf,
63 module Data.Complex,
64 -- * IO
65 module Data.Packed.IO,
66 -- * Misc
67 Testable(..)
68) where
69
70import Data.Packed
71import Data.Packed.Internal.Numeric
72import Data.Complex
73import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
74import Data.Monoid(Monoid(mconcat))
75import Data.Packed.IO
76import Numeric.LinearAlgebra.Random
77
78------------------------------------------------------------------
79
80{- | Creates a real vector containing a range of values:
81
82>>> linspace 5 (-3,7::Double)
83fromList [-3.0,-0.5,2.0,4.5,7.0]@
84
85>>> linspace 5 (8,2+i) :: Vector (Complex Double)
86fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0]
87
88Logarithmic spacing can be defined as follows:
89
90@logspace n (a,b) = 10 ** linspace n (a,b)@
91-}
92linspace :: (Container Vector e) => Int -> (e, e) -> Vector e
93linspace 0 (a,b) = fromList[(a+b)/2]
94linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1]
95 where s = (b-a)/fromIntegral (n-1)
96
97--------------------------------------------------------
98
99{- | Matrix product, matrix - vector product, and dot product (equivalent to 'contraction')
100
101(This operator can also be written using the unicode symbol ◇ (25c7).)
102
103Examples:
104
105>>> let a = (3><4) [1..] :: Matrix Double
106>>> let v = fromList [1,0,2,-1] :: Vector Double
107>>> let u = fromList [1,2,3] :: Vector Double
108
109>>> a
110(3><4)
111 [ 1.0, 2.0, 3.0, 4.0
112 , 5.0, 6.0, 7.0, 8.0
113 , 9.0, 10.0, 11.0, 12.0 ]
114
115matrix × matrix:
116
117>>> disp 2 (a <.> trans a)
1183x3
119 30 70 110
120 70 174 278
121110 278 446
122
123matrix × vector:
124
125>>> a <.> v
126fromList [3.0,11.0,19.0]
127
128dot product:
129
130>>> u <.> fromList[3,2,1::Double]
13110
132
133For complex vectors the first argument is conjugated:
134
135>>> fromList [1,i] <.> fromList[2*i+1,3]
1361.0 :+ (-1.0)
137
138>>> fromList [1,i,1-i] <.> complex a
139fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0]
140-}
141infixl 7 <.>
142(<.>) :: Contraction a b c => a -> b -> c
143(<.>) = contraction
144
145
146class Contraction a b c | a b -> c
147 where
148 -- | Matrix product, matrix - vector product, and dot product
149 contraction :: a -> b -> c
150
151instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where
152 u `contraction` v = conj u `udot` v
153
154instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where
155 contraction = mXv
156
157instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where
158 contraction v m = (conj v) `vXm` m
159
160instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where
161 contraction = mXm
162
163
164--------------------------------------------------------------------------------
165
166class Mul a b c | a b -> c where
167 infixl 7 <>
168 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
169 (<>) :: Product t => a t -> b t -> c t
170
171instance Mul Matrix Matrix Matrix where
172 (<>) = mXm
173
174instance Mul Matrix Vector Vector where
175 (<>) m v = flatten $ m <> asColumn v
176
177instance Mul Vector Matrix Vector where
178 (<>) v m = flatten $ asRow v <> m
179
180--------------------------------------------------------------------------------
181
182-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
183infixl 7 <\>
184(<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t
185(<\>) = linSolve
186
187class LSDiv c
188 where
189 linSolve :: Field t => Matrix t -> c t -> c t
190
191instance LSDiv Vector
192 where
193 linSolve m v = flatten (linearSolveSVD m (reshape 1 v))
194
195instance LSDiv Matrix
196 where
197 linSolve = linearSolveSVD
198
199--------------------------------------------------------------------------------
200
201class Konst e d c | d -> c, c -> d
202 where
203 -- |
204 -- >>> konst 7 3 :: Vector Float
205 -- fromList [7.0,7.0,7.0]
206 --
207 -- >>> konst i (3::Int,4::Int)
208 -- (3><4)
209 -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
210 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
211 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ]
212 --
213 konst :: e -> d -> c e
214
215instance Container Vector e => Konst e Int Vector
216 where
217 konst = konst'
218
219instance Container Vector e => Konst e (Int,Int) Matrix
220 where
221 konst = konst'
222
223--------------------------------------------------------------------------------
224
225class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
226 where
227 -- |
228 -- >>> build 5 (**2) :: Vector Double
229 -- fromList [0.0,1.0,4.0,9.0,16.0]
230 --
231 -- Hilbert matrix of order N:
232 --
233 -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
234 -- >>> putStr . dispf 2 $ hilb 3
235 -- 3x3
236 -- 1.00 0.50 0.33
237 -- 0.50 0.33 0.25
238 -- 0.33 0.25 0.20
239 --
240 build :: d -> f -> c e
241
242instance Container Vector e => Build Int (e -> e) Vector e
243 where
244 build = build'
245
246instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e
247 where
248 build = build'
249
250--------------------------------------------------------------------------------
251
252-- | alternative unicode symbol (25c7) for 'contraction'
253(◇) :: Contraction a b c => a -> b -> c
254infixl 7 ◇
255(◇) = contraction
256
257-- | dot product: @cdot u v = 'udot' ('conj' u) v@
258dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t
259dot u v = udot (conj u) v
260
261--------------------------------------------------------------------------------
262
263optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t
264optimiseMult = mconcat
265
266--------------------------------------------------------------------------------
267
268
269{- | Compute mean vector and covariance matrix of the rows of a matrix.
270
271>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3])
272(fromList [4.010341078059521,5.0197204699640405],
273(2><2)
274 [ 1.9862461923890056, -1.0127225830525157e-2
275 , -1.0127225830525157e-2, 3.0373954915729318 ])
276
277-}
278meanCov :: Matrix Double -> (Vector Double, Matrix Double)
279meanCov x = (med,cov) where
280 r = rows x
281 k = 1 / fromIntegral r
282 med = konst k r `vXm` x
283 meds = konst 1 r `outer` med
284 xc = x `sub` meds
285 cov = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc)
286
287--------------------------------------------------------------------------------
288
diff --git a/packages/base/src/Data/Packed/Vector.hs b/packages/base/src/Data/Packed/Vector.hs
index 53fe563..31dcf47 100644
--- a/packages/base/src/Data/Packed/Vector.hs
+++ b/packages/base/src/Data/Packed/Vector.hs
@@ -17,17 +17,15 @@
17 17
18module Data.Packed.Vector ( 18module Data.Packed.Vector (
19 Vector, 19 Vector,
20 fromList, (|>), toList, buildVector, constant, 20 fromList, (|>), toList, buildVector,
21 dim, (@>), 21 dim, (@>),
22 subVector, takesV, vjoin, join, 22 subVector, takesV, vjoin, join,
23 mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith, 23 mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith,
24 mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_, 24 mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_,
25 foldLoop, foldVector, foldVectorG, foldVectorWithIndex, 25 foldLoop, foldVector, foldVectorG, foldVectorWithIndex
26 stepD, stepF, condD, condF, conjugateC, conjugateQ
27) where 26) where
28 27
29import Data.Packed.Internal.Vector 28import Data.Packed.Internal.Vector
30import Data.Packed.Internal.Matrix
31import Foreign.Storable 29import Foreign.Storable
32 30
33------------------------------------------------------------------- 31-------------------------------------------------------------------
@@ -95,13 +93,3 @@ unzipVector = unzipVectorWith id
95join :: Storable t => [Vector t] -> Vector t 93join :: Storable t => [Vector t] -> Vector t
96join = vjoin 94join = vjoin
97 95
98{- | creates a vector with a given number of equal components:
99
100@> constant 2 7
1017 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
102-}
103constant :: Element a => a -> Int -> Vector a
104-- constant x n = runSTVector (newVector x n)
105constant = constantD-- about 2x faster
106
107
diff --git a/packages/base/src/Numeric/Container.hs b/packages/base/src/Numeric/Container.hs
index 067c5fa..6a841aa 100644
--- a/packages/base/src/Numeric/Container.hs
+++ b/packages/base/src/Numeric/Container.hs
@@ -1,258 +1,49 @@
1{-# LANGUAGE TypeFamilies #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE UndecidableInstances #-}
7
8-----------------------------------------------------------------------------
9-- |
10-- Module : Numeric.Container
11-- Copyright : (c) Alberto Ruiz 2010-14
12-- License : BSD3
13-- Maintainer : Alberto Ruiz
14-- Stability : provisional
15--
16-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines.
17--
18-- The 'Container' class is used to define optimized generic functions which work
19-- on 'Vector' and 'Matrix' with real or complex elements.
20--
21-- Some of these functions are also available in the instances of the standard
22-- numeric Haskell classes provided by "Numeric.LinearAlgebra".
23--
24-----------------------------------------------------------------------------
25{-# OPTIONS_HADDOCK hide #-} 1{-# OPTIONS_HADDOCK hide #-}
26 2
27module Numeric.Container ( 3module Numeric.Container(
28 -- * Basic functions
29 module Data.Packed, 4 module Data.Packed,
30 konst, build, 5 constant,
31 linspace, 6 linspace,
32 diag, ident, 7 diag,
8 ident,
33 ctrans, 9 ctrans,
34 -- * Generic operations 10 Container(scaleRecip, addConstant,add, sub, mul, divide, equal),
35 Container, 11 scalar,
36 add, mul, sub, divide, equal, scaleRecip, addConstant, 12 conj,
37 scalar, conj, scale, arctan2, cmap, 13 scale,
38 atIndex, minIndex, maxIndex, minElement, maxElement, 14 arctan2,
15 cmap,
16 Konst(..),
17 Build(..),
18 atIndex,
19 minIndex, maxIndex, minElement, maxElement,
39 sumElements, prodElements, 20 sumElements, prodElements,
40 step, cond, find, assoc, accum, 21 step, cond, find, assoc, accum,
41 Transposable(..), Linear(..), 22 Element(..),
42 -- * Matrix product 23 Product(..),
43 Product(..), udot, dot, (◇),
44 Mul(..),
45 Contraction(..),(<.>),
46 optimiseMult, 24 optimiseMult,
47 mXm,mXv,vXm,LSDiv(..), 25 mXm, mXv, vXm, (<.>),
26 Mul(..),
27 LSDiv, (<\>),
48 outer, kronecker, 28 outer, kronecker,
49 -- * Random numbers
50 RandDist(..), 29 RandDist(..),
51 randomVector, 30 randomVector, gaussianSample, uniformSample,
52 gaussianSample, 31 meanCov,
53 uniformSample,
54 -- * Element conversion
55 Convert(..), 32 Convert(..),
56 Complexable(), 33 Complexable,
57 RealElement(), 34 RealElement,
58 35 RealOf, ComplexOf, SingleOf, DoubleOf, IndexOf,
59 RealOf, ComplexOf, SingleOf, DoubleOf,
60
61 IndexOf,
62 module Data.Complex, 36 module Data.Complex,
63 -- * IO 37 dispf, disps, dispcf, vecdisp, latexFormat, format,
64 module Data.Packed.IO, 38 loadMatrix, saveMatrix, readMatrix
65 -- * Misc
66 Testable(..)
67) where 39) where
68 40
69import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ)
70import Data.Packed.Internal.Numeric
71import Data.Complex
72import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
73import Data.Monoid(Monoid(mconcat))
74import Data.Packed.IO
75import Numeric.LinearAlgebra.Random
76
77------------------------------------------------------------------
78
79{- | Creates a real vector containing a range of values:
80
81>>> linspace 5 (-3,7::Double)
82fromList [-3.0,-0.5,2.0,4.5,7.0]@
83
84>>> linspace 5 (8,2+i) :: Vector (Complex Double)
85fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0]
86
87Logarithmic spacing can be defined as follows:
88
89@logspace n (a,b) = 10 ** linspace n (a,b)@
90-}
91linspace :: (Container Vector e) => Int -> (e, e) -> Vector e
92linspace 0 (a,b) = fromList[(a+b)/2]
93linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1]
94 where s = (b-a)/fromIntegral (n-1)
95
96--------------------------------------------------------
97
98{- | Matrix product, matrix - vector product, and dot product (equivalent to 'contraction')
99
100(This operator can also be written using the unicode symbol ◇ (25c7).)
101
102Examples:
103
104>>> let a = (3><4) [1..] :: Matrix Double
105>>> let v = fromList [1,0,2,-1] :: Vector Double
106>>> let u = fromList [1,2,3] :: Vector Double
107
108>>> a
109(3><4)
110 [ 1.0, 2.0, 3.0, 4.0
111 , 5.0, 6.0, 7.0, 8.0
112 , 9.0, 10.0, 11.0, 12.0 ]
113
114matrix × matrix:
115
116>>> disp 2 (a <.> trans a)
1173x3
118 30 70 110
119 70 174 278
120110 278 446
121
122matrix × vector:
123
124>>> a <.> v
125fromList [3.0,11.0,19.0]
126
127dot product:
128
129>>> u <.> fromList[3,2,1::Double]
13010
131
132For complex vectors the first argument is conjugated:
133
134>>> fromList [1,i] <.> fromList[2*i+1,3]
1351.0 :+ (-1.0)
136
137>>> fromList [1,i,1-i] <.> complex a
138fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0]
139-}
140infixl 7 <.>
141(<.>) :: Contraction a b c => a -> b -> c
142(<.>) = contraction
143
144 41
145class Contraction a b c | a b -> c 42import Data.Packed.Numeric
146 where 43import Data.Packed
147 -- | Matrix product, matrix - vector product, and dot product 44import Data.Packed.Internal(constantD)
148 contraction :: a -> b -> c 45import Data.Complex
149
150instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where
151 u `contraction` v = conj u `udot` v
152
153instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where
154 contraction = mXv
155
156instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where
157 contraction v m = (conj v) `vXm` m
158
159instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where
160 contraction = mXm
161
162
163--------------------------------------------------------------------------------
164
165class Mul a b c | a b -> c where
166 infixl 7 <>
167 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
168 (<>) :: Product t => a t -> b t -> c t
169
170instance Mul Matrix Matrix Matrix where
171 (<>) = mXm
172
173instance Mul Matrix Vector Vector where
174 (<>) m v = flatten $ m <> asColumn v
175
176instance Mul Vector Matrix Vector where
177 (<>) v m = flatten $ asRow v <> m
178
179--------------------------------------------------------------------------------
180
181class LSDiv c where
182 infixl 7 <\>
183 -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
184 (<\>) :: Field t => Matrix t -> c t -> c t
185
186instance LSDiv Vector where
187 m <\> v = flatten (linearSolveSVD m (reshape 1 v))
188
189instance LSDiv Matrix where
190 (<\>) = linearSolveSVD
191
192--------------------------------------------------------------------------------
193
194class Konst e d c | d -> c, c -> d
195 where
196 -- |
197 -- >>> konst 7 3 :: Vector Float
198 -- fromList [7.0,7.0,7.0]
199 --
200 -- >>> konst i (3::Int,4::Int)
201 -- (3><4)
202 -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
203 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
204 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ]
205 --
206 konst :: e -> d -> c e
207
208instance Container Vector e => Konst e Int Vector
209 where
210 konst = konst'
211
212instance Container Vector e => Konst e (Int,Int) Matrix
213 where
214 konst = konst'
215
216--------------------------------------------------------------------------------
217
218class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
219 where
220 -- |
221 -- >>> build 5 (**2) :: Vector Double
222 -- fromList [0.0,1.0,4.0,9.0,16.0]
223 --
224 -- Hilbert matrix of order N:
225 --
226 -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
227 -- >>> putStr . dispf 2 $ hilb 3
228 -- 3x3
229 -- 1.00 0.50 0.33
230 -- 0.50 0.33 0.25
231 -- 0.33 0.25 0.20
232 --
233 build :: d -> f -> c e
234
235instance Container Vector e => Build Int (e -> e) Vector e
236 where
237 build = build'
238
239instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e
240 where
241 build = build'
242
243--------------------------------------------------------------------------------
244
245-- | alternative unicode symbol (25c7) for 'contraction'
246(◇) :: Contraction a b c => a -> b -> c
247infixl 7 ◇
248(◇) = contraction
249
250-- | dot product: @cdot u v = 'udot' ('conj' u) v@
251dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t
252dot u v = udot (conj u) v
253
254--------------------------------------------------------------------------------
255 46
256optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t 47constant :: Element a => a -> Int -> Vector a
257optimiseMult = mconcat 48constant = constantD
258 49
diff --git a/packages/base/src/Numeric/HMatrix.hs b/packages/base/src/Numeric/HMatrix.hs
new file mode 100644
index 0000000..a56c3d2
--- /dev/null
+++ b/packages/base/src/Numeric/HMatrix.hs
@@ -0,0 +1,158 @@
1-----------------------------------------------------------------------------
2{- |
3Module : Numeric.HMatrix
4Copyright : (c) Alberto Ruiz 2006-14
5License : BSD3
6Maintainer : Alberto Ruiz
7Stability : provisional
8
9-}
10-----------------------------------------------------------------------------
11module Numeric.HMatrix (
12
13 -- * Basic types and data processing
14 module Numeric.LinearAlgebra.Data,
15
16 -- * Arithmetic and numeric classes
17 -- |
18 -- The standard numeric classes are defined elementwise:
19 --
20 -- >>> fromList [1,2,3] * fromList [3,0,-2 :: Double]
21 -- fromList [3.0,0.0,-6.0]
22 --
23 -- >>> (3><3) [1..9] * ident 3 :: Matrix Double
24 -- (3><3)
25 -- [ 1.0, 0.0, 0.0
26 -- , 0.0, 5.0, 0.0
27 -- , 0.0, 0.0, 9.0 ]
28 --
29 -- In arithmetic operations single-element vectors and matrices
30 -- (created from numeric literals or using 'scalar') automatically
31 -- expand to match the dimensions of the other operand:
32 --
33 -- >>> 5 + 2*ident 3 :: Matrix Double
34 -- (3><3)
35 -- [ 7.0, 5.0, 5.0
36 -- , 5.0, 7.0, 5.0
37 -- , 5.0, 5.0, 7.0 ]
38 --
39
40 -- * Matrix product
41 (<.>),
42
43 -- | The overloaded multiplication operators may need type annotations to remove
44 -- ambiguity. In those cases we can also use the specific functions 'mXm', 'mXv', and 'dot'.
45 --
46 -- The matrix x matrix product is also implemented in the "Data.Monoid" instance, where
47 -- single-element matrices (created from numeric literals or using 'scalar')
48 -- are used for scaling.
49 --
50 -- >>> let m = (2><3)[1..] :: Matrix Double
51 -- >>> m <> 2 <> diagl[0.5,1,0]
52 -- (2><3)
53 -- [ 1.0, 4.0, 0.0
54 -- , 4.0, 10.0, 0.0 ]
55 --
56 -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
57
58
59 -- * Other products
60 outer, kronecker, cross,
61 scale,
62 sumElements, prodElements,
63
64 -- * Linear Systems
65 (<\>),
66 linearSolve,
67 linearSolveLS,
68 linearSolveSVD,
69 luSolve,
70 cholSolve,
71 cgSolve,
72 cgSolve',
73
74 -- * Inverse and pseudoinverse
75 inv, pinv, pinvTol,
76
77 -- * Determinant and rank
78 rcond, rank, ranksv,
79 det, invlndet,
80
81 -- * Singular value decomposition
82 svd,
83 fullSVD,
84 thinSVD,
85 compactSVD,
86 singularValues,
87 leftSV, rightSV,
88
89 -- * Eigensystems
90 eig, eigSH, eigSH',
91 eigenvalues, eigenvaluesSH, eigenvaluesSH',
92 geigSH',
93
94 -- * QR
95 qr, rq, qrRaw, qrgr,
96
97 -- * Cholesky
98 chol, cholSH, mbCholSH,
99
100 -- * Hessenberg
101 hess,
102
103 -- * Schur
104 schur,
105
106 -- * LU
107 lu, luPacked,
108
109 -- * Matrix functions
110 expm,
111 sqrtm,
112 matFunc,
113
114 -- * Nullspace
115 nullspacePrec,
116 nullVector,
117 nullspaceSVD,
118 null1, null1sym,
119
120 orth,
121
122 -- * Norms
123 norm_0, norm_1, norm_2, norm_Inf,
124 mnorm_0, mnorm_1, mnorm_2, mnorm_Inf,
125 norm_Frob, norm_nuclear,
126
127 -- * Correlation and convolution
128 corr, conv, corrMin, corr2, conv2,
129
130 -- * Random arrays
131
132 RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
133
134 -- * Misc
135 meanCov, peps, relativeError, haussholder, optimiseMult, dot, udot, mXm, mXv, smXv, (<>), (◇), Seed, checkT,
136 -- * Auxiliary classes
137 Element, Container, Product, Contraction, LSDiv,
138 Complexable, RealElement,
139 RealOf, ComplexOf, SingleOf, DoubleOf,
140 IndexOf,
141 Field,
142 Normed,
143 CGMat, Transposable,
144 ℕ,ℤ,ℝ,ℂ,ℝn,ℂn, 𝑖, i_C --ℍ
145) where
146
147import Numeric.LinearAlgebra.Data
148
149import Numeric.Matrix()
150import Numeric.Vector()
151import Data.Packed.Numeric
152import Numeric.LinearAlgebra.Algorithms
153import Numeric.LinearAlgebra.Util
154import Numeric.LinearAlgebra.Random
155import Numeric.Sparse(smXv)
156import Numeric.LinearAlgebra.Util.CG
157
158
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index 242122f..ad315e4 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -1,4 +1,4 @@
1----------------------------------------------------------------------------- 1--------------------------------------------------------------------------------
2{- | 2{- |
3Module : Numeric.LinearAlgebra 3Module : Numeric.LinearAlgebra
4Copyright : (c) Alberto Ruiz 2006-14 4Copyright : (c) Alberto Ruiz 2006-14
@@ -7,149 +7,16 @@ Maintainer : Alberto Ruiz
7Stability : provisional 7Stability : provisional
8 8
9-} 9-}
10----------------------------------------------------------------------------- 10--------------------------------------------------------------------------------
11module Numeric.LinearAlgebra ( 11{-# OPTIONS_HADDOCK hide #-}
12
13 -- * Basic types and data processing
14 module Numeric.LinearAlgebra.Data,
15
16 -- * Arithmetic and numeric classes
17 -- |
18 -- The standard numeric classes are defined elementwise:
19 --
20 -- >>> fromList [1,2,3] * fromList [3,0,-2 :: Double]
21 -- fromList [3.0,0.0,-6.0]
22 --
23 -- >>> (3><3) [1..9] * ident 3 :: Matrix Double
24 -- (3><3)
25 -- [ 1.0, 0.0, 0.0
26 -- , 0.0, 5.0, 0.0
27 -- , 0.0, 0.0, 9.0 ]
28 --
29 -- In arithmetic operations single-element vectors and matrices
30 -- (created from numeric literals or using 'scalar') automatically
31 -- expand to match the dimensions of the other operand:
32 --
33 -- >>> 5 + 2*ident 3 :: Matrix Double
34 -- (3><3)
35 -- [ 7.0, 5.0, 5.0
36 -- , 5.0, 7.0, 5.0
37 -- , 5.0, 5.0, 7.0 ]
38 --
39
40 -- * Matrix product
41 (<.>),
42
43 -- | The overloaded multiplication operators may need type annotations to remove
44 -- ambiguity. In those cases we can also use the specific functions 'mXm', 'mXv', and 'dot'.
45 --
46 -- The matrix x matrix product is also implemented in the "Data.Monoid" instance, where
47 -- single-element matrices (created from numeric literals or using 'scalar')
48 -- are used for scaling.
49 --
50 -- >>> let m = (2><3)[1..] :: Matrix Double
51 -- >>> m <> 2 <> diagl[0.5,1,0]
52 -- (2><3)
53 -- [ 1.0, 4.0, 0.0
54 -- , 4.0, 10.0, 0.0 ]
55 --
56 -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
57
58
59 -- * Other products
60 outer, kronecker, cross,
61 scale,
62 sumElements, prodElements, absSum,
63
64 -- * Linear Systems
65 (<\>),
66 linearSolve,
67 linearSolveLS,
68 linearSolveSVD,
69 luSolve,
70 cholSolve,
71 cgSolve,
72 cgSolve',
73
74 -- * Inverse and pseudoinverse
75 inv, pinv, pinvTol,
76
77 -- * Determinant and rank
78 rcond, rank, ranksv,
79 det, invlndet,
80
81 -- * Singular value decomposition
82 svd,
83 fullSVD,
84 thinSVD,
85 compactSVD,
86 singularValues,
87 leftSV, rightSV,
88
89 -- * Eigensystems
90 eig, eigSH, eigSH',
91 eigenvalues, eigenvaluesSH, eigenvaluesSH',
92 geigSH',
93
94 -- * QR
95 qr, rq, qrRaw, qrgr,
96
97 -- * Cholesky
98 chol, cholSH, mbCholSH,
99
100 -- * Hessenberg
101 hess,
102
103 -- * Schur
104 schur,
105
106 -- * LU
107 lu, luPacked,
108
109 -- * Matrix functions
110 expm,
111 sqrtm,
112 matFunc,
113 12
114 -- * Nullspace 13module Numeric.LinearAlgebra (
115 nullspacePrec, 14 module Numeric.Container,
116 nullVector, 15 module Numeric.LinearAlgebra.Algorithms
117 nullspaceSVD,
118 null1, null1sym,
119
120 orth,
121
122 -- * Norms
123 norm1, norm2, normInf, pnorm, NormType(..),
124
125 -- * Correlation and convolution
126 corr, conv, corrMin, corr2, conv2,
127
128 -- * Random arrays
129
130 RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
131
132 -- * Misc
133 meanCov, peps, relativeError, haussholder, optimiseMult, dot, udot, mXm, mXv, smXv, (<>), (◇), Seed, checkT,
134 -- * Auxiliary classes
135 Element, Container, Product, Contraction, LSDiv,
136 Complexable(), RealElement(),
137 RealOf, ComplexOf, SingleOf, DoubleOf,
138 IndexOf,
139 Field, Normed,
140 CGMat, Transposable,
141 R,V
142) where 16) where
143 17
144import Numeric.LinearAlgebra.Data
145
146import Numeric.Matrix()
147import Numeric.Vector()
148import Numeric.Container 18import Numeric.Container
149import Numeric.LinearAlgebra.Algorithms 19import Numeric.LinearAlgebra.Algorithms
150import Numeric.LinearAlgebra.Util 20import Numeric.Matrix()
151import Numeric.LinearAlgebra.Random 21import Numeric.Vector()
152import Numeric.Sparse(smXv)
153import Numeric.LinearAlgebra.Util.CG
154
155 22
diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
index c7e7043..bbcc513 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs
@@ -76,12 +76,12 @@ module Numeric.LinearAlgebra.Algorithms (
76) where 76) where
77 77
78 78
79import Data.Packed.Development hiding ((//))
80import Data.Packed 79import Data.Packed
81import Numeric.LinearAlgebra.LAPACK as LAPACK 80import Numeric.LinearAlgebra.LAPACK as LAPACK
82import Data.List(foldl1') 81import Data.List(foldl1')
83import Data.Array 82import Data.Array
84import Data.Packed.Internal.Numeric 83import Data.Packed.Internal.Numeric
84import Data.Packed.Internal(shSize)
85 85
86 86
87{- | Generic linear algebra functions for double precision real and complex matrices. 87{- | Generic linear algebra functions for double precision real and complex matrices.
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index e3cbe31..89bebbe 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -16,13 +16,19 @@ module Numeric.LinearAlgebra.Data(
16 -- * Vector 16 -- * Vector
17 -- | 1D arrays are storable vectors from the vector package. 17 -- | 1D arrays are storable vectors from the vector package.
18 18
19 Vector, (|>), dim, (@>), 19 vect, (|>),
20 20
21 -- * Matrix 21 -- * Matrix
22 Matrix, (><), size, (@@>), trans, ctrans, 22
23 mat, (><), tr,
24
25 -- * Indexing
26
27 size,
28 Indexable(..),
23 29
24 -- * Construction 30 -- * Construction
25 scalar, konst, build, assoc, accum, linspace, -- ones, zeros, 31 scalar, Konst(..), Build(..), assoc, accum, linspace, -- ones, zeros,
26 32
27 -- * Diagonal 33 -- * Diagonal
28 ident, diag, diagl, diagRect, takeDiag, 34 ident, diag, diagl, diagRect, takeDiag,
@@ -62,11 +68,13 @@ module Numeric.LinearAlgebra.Data(
62 68
63 module Data.Complex, 69 module Data.Complex,
64 70
71 Vector, Matrix
72
65) where 73) where
66 74
67import Data.Packed.Vector 75import Data.Packed.Vector
68import Data.Packed.Matrix 76import Data.Packed.Matrix
69import Numeric.Container 77import Data.Packed.Numeric
70import Numeric.LinearAlgebra.Util 78import Numeric.LinearAlgebra.Util
71import Data.Complex 79import Data.Complex
72import Numeric.Sparse 80import Numeric.Sparse
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs
index 2f91e18..a7d6946 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs
@@ -1,4 +1,9 @@
1{-# LANGUAGE FlexibleContexts #-} 1{-# LANGUAGE FlexibleContexts #-}
2{-# LANGUAGE FlexibleInstances #-}
3{-# LANGUAGE TypeFamilies #-}
4{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE FunctionalDependencies #-}
6
2----------------------------------------------------------------------------- 7-----------------------------------------------------------------------------
3{- | 8{- |
4Module : Numeric.LinearAlgebra.Util 9Module : Numeric.LinearAlgebra.Util
@@ -14,19 +19,24 @@ Stability : provisional
14module Numeric.LinearAlgebra.Util( 19module Numeric.LinearAlgebra.Util(
15 20
16 -- * Convenience functions 21 -- * Convenience functions
17 size, disp, 22 vect, mat,
23 disp,
18 zeros, ones, 24 zeros, ones,
19 diagl, 25 diagl,
20 row, 26 row,
21 col, 27 col,
22 (&), (¦), (——), (#), 28 (&), (¦), (——), (#),
23 (?), (¿), 29 (?), (¿),
30 Indexable(..), size,
31 rand, randn,
24 cross, 32 cross,
25 norm, 33 norm,
34 ℕ,ℤ,ℝ,ℂ,ℝn,ℂn,𝑖,i_C, --ℍ
35 norm_1, norm_2, norm_0, norm_Inf, norm_Frob, norm_nuclear,
36 mnorm_1, mnorm_2, mnorm_0, mnorm_Inf,
26 unitary, 37 unitary,
27 mt, 38 mt,
28 pairwiseD2, 39 pairwiseD2,
29 meanCov,
30 rowOuters, 40 rowOuters,
31 null1, 41 null1,
32 null1sym, 42 null1sym,
@@ -48,13 +58,49 @@ module Numeric.LinearAlgebra.Util(
48 vtrans 58 vtrans
49) where 59) where
50 60
51import Numeric.Container 61import Data.Packed.Numeric
52import Numeric.LinearAlgebra.Algorithms hiding (i) 62import Numeric.LinearAlgebra.Algorithms hiding (i)
53import Numeric.Matrix() 63import Numeric.Matrix()
54import Numeric.Vector() 64import Numeric.Vector()
55 65import Numeric.LinearAlgebra.Random
56import Numeric.LinearAlgebra.Util.Convolution 66import Numeric.LinearAlgebra.Util.Convolution
57 67
68type ℝ = Double
69type ℕ = Int
70type ℤ = Int
71type ℂ = Complex Double
72type ℝn = Vector ℝ
73type ℂn = Vector ℂ
74--newtype ℍ m = H m
75
76i_C, 𝑖 :: ℂ
77𝑖 = 0:+1
78i_C = 𝑖
79
80{- | create a real vector
81
82>>> vect [1..5]
83fromList [1.0,2.0,3.0,4.0,5.0]
84
85-}
86vect :: [ℝ] -> ℝn
87vect = fromList
88
89{- | create a real matrix
90
91>>> mat 5 [1..15]
92(3><5)
93 [ 1.0, 2.0, 3.0, 4.0, 5.0
94 , 6.0, 7.0, 8.0, 9.0, 10.0
95 , 11.0, 12.0, 13.0, 14.0, 15.0 ]
96
97-}
98mat
99 :: Int -- ^ columns
100 -> [ℝ] -- ^ elements
101 -> Matrix ℝ
102mat c = reshape c . fromList
103
58{- | print a real matrix with given number of digits after the decimal point 104{- | print a real matrix with given number of digits after the decimal point
59 105
60>>> disp 5 $ ident 2 / 3 106>>> disp 5 $ ident 2 / 3
@@ -175,38 +221,97 @@ norm :: Vector Double -> Double
175-- ^ 2-norm of real vector 221-- ^ 2-norm of real vector
176norm = pnorm PNorm2 222norm = pnorm PNorm2
177 223
224norm_2 :: Normed Vector t => Vector t -> RealOf t
225norm_2 = pnorm PNorm2
226
227norm_1 :: Normed Vector t => Vector t -> RealOf t
228norm_1 = pnorm PNorm1
229
230norm_Inf :: Normed Vector t => Vector t -> RealOf t
231norm_Inf = pnorm Infinity
232
233norm_0 :: Vector ℝ -> ℝ
234norm_0 v = sumElements (step (abs v - scalar (eps*mx)))
235 where
236 mx = norm_Inf v
237
238norm_Frob :: Normed Matrix t => Matrix t -> RealOf t
239norm_Frob = pnorm Frobenius
240
241norm_nuclear :: Field t => Matrix t -> ℝ
242norm_nuclear = sumElements . singularValues
243
244mnorm_2 :: Normed Matrix t => Matrix t -> RealOf t
245mnorm_2 = pnorm PNorm2
246
247mnorm_1 :: Normed Matrix t => Matrix t -> RealOf t
248mnorm_1 = pnorm PNorm1
249
250mnorm_Inf :: Normed Matrix t => Matrix t -> RealOf t
251mnorm_Inf = pnorm Infinity
252
253mnorm_0 :: Matrix ℝ -> ℝ
254mnorm_0 = norm_0 . flatten
178 255
179-- | Obtains a vector in the same direction with 2-norm=1 256-- | Obtains a vector in the same direction with 2-norm=1
180unitary :: Vector Double -> Vector Double 257unitary :: Vector Double -> Vector Double
181unitary v = v / scalar (norm v) 258unitary v = v / scalar (norm v)
182 259
183-- | ('rows' &&& 'cols')
184size :: Matrix t -> (Int, Int)
185size m = (rows m, cols m)
186 260
187-- | trans . inv 261-- | trans . inv
188mt :: Matrix Double -> Matrix Double 262mt :: Matrix Double -> Matrix Double
189mt = trans . inv 263mt = trans . inv
190 264
191-------------------------------------------------------------------------------- 265--------------------------------------------------------------------------------
266{- |
267
268>>> size $ fromList[1..10::Double]
26910
270>>> size $ (2><5)[1..10::Double]
271(2,5)
272
273-}
274size :: Container c t => c t -> IndexOf c
275size = size'
192 276
193{- | Compute mean vector and covariance matrix of the rows of a matrix. 277{- |
278
279>>> vect [1..10] ! 3
2804.0
281
282>>> mat 5 [1..15] ! 1
283fromList [6.0,7.0,8.0,9.0,10.0]
194 284
195>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3]) 285>>> mat 5 [1..15] ! 1 ! 3
196(fromList [4.010341078059521,5.0197204699640405], 2869.0
197(2><2)
198 [ 1.9862461923890056, -1.0127225830525157e-2
199 , -1.0127225830525157e-2, 3.0373954915729318 ])
200 287
201-} 288-}
202meanCov :: Matrix Double -> (Vector Double, Matrix Double) 289class Indexable c t | c -> t , t -> c
203meanCov x = (med,cov) where 290 where
204 r = rows x 291 infixl 9 !
205 k = 1 / fromIntegral r 292 (!) :: c -> Int -> t
206 med = konst k r `vXm` x 293
207 meds = konst 1 r `outer` med 294instance Indexable (Vector Double) Double
208 xc = x `sub` meds 295 where
209 cov = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc) 296 (!) = (@>)
297
298instance Indexable (Vector Float) Float
299 where
300 (!) = (@>)
301
302instance Indexable (Vector (Complex Double)) (Complex Double)
303 where
304 (!) = (@>)
305
306instance Indexable (Vector (Complex Float)) (Complex Float)
307 where
308 (!) = (@>)
309
310instance Element t => Indexable (Matrix t) (Vector t)
311 where
312 m!j = subVector (j*c) c (flatten m)
313 where
314 c = cols m
210 315
211-------------------------------------------------------------------------------- 316--------------------------------------------------------------------------------
212 317
@@ -220,7 +325,7 @@ pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y
220 ox = one (rows x) 325 ox = one (rows x)
221 oy = one (rows y) 326 oy = one (rows y)
222 oc = one (cols x) 327 oc = one (cols x)
223 one k = constant 1 k 328 one k = konst 1 k
224 x2 = x * x <> oc 329 x2 = x * x <> oc
225 y2 = y * y <> oc 330 y2 = y * y <> oc
226 ok = cols x == cols y 331 ok = cols x == cols y
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
index d21602d..5e2ea84 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util/CG.hs
@@ -6,7 +6,7 @@ module Numeric.LinearAlgebra.Util.CG(
6 CGMat, CGState(..), R, V 6 CGMat, CGState(..), R, V
7) where 7) where
8 8
9import Numeric.Container 9import Data.Packed.Numeric
10import Numeric.Vector() 10import Numeric.Vector()
11 11
12{- 12{-
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs
index e4cba8f..c8c7536 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs
@@ -16,16 +16,18 @@ module Numeric.LinearAlgebra.Util.Convolution(
16 corr2, conv2, separable 16 corr2, conv2, separable
17) where 17) where
18 18
19import Numeric.Container 19import Data.Packed.Numeric
20 20
21 21
22vectSS :: Element t => Int -> Vector t -> Matrix t 22vectSS :: Element t => Int -> Vector t -> Matrix t
23vectSS n v = fromRows [ subVector k n v | k <- [0 .. dim v - n] ] 23vectSS n v = fromRows [ subVector k n v | k <- [0 .. dim v - n] ]
24 24
25 25
26corr :: Product t => Vector t -- ^ kernel 26corr
27 -> Vector t -- ^ source 27 :: (Container Vector t, Product t)
28 -> Vector t 28 => Vector t -- ^ kernel
29 -> Vector t -- ^ source
30 -> Vector t
29{- ^ correlation 31{- ^ correlation
30 32
31>>> corr (fromList[1,2,3]) (fromList [1..10]) 33>>> corr (fromList[1,2,3]) (fromList [1..10])
@@ -33,12 +35,12 @@ fromList [14.0,20.0,26.0,32.0,38.0,44.0,50.0,56.0]
33 35
34-} 36-}
35corr ker v 37corr ker v
36 | dim ker == 0 = constant 0 (dim v) 38 | dim ker == 0 = konst 0 (dim v)
37 | dim ker <= dim v = vectSS (dim ker) v <> ker 39 | dim ker <= dim v = vectSS (dim ker) v <> ker
38 | otherwise = error $ "corr: dim kernel ("++show (dim ker)++") > dim vector ("++show (dim v)++")" 40 | otherwise = error $ "corr: dim kernel ("++show (dim ker)++") > dim vector ("++show (dim v)++")"
39 41
40 42
41conv :: (Product t, Num t) => Vector t -> Vector t -> Vector t 43conv :: (Container Vector t, Product t, Num t) => Vector t -> Vector t -> Vector t
42{- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product) 44{- ^ convolution ('corr' with reversed kernel and padded input, equivalent to polynomial product)
43 45
44>>> conv (fromList[1,1]) (fromList [-1,1]) 46>>> conv (fromList[1,1]) (fromList [-1,1])
@@ -46,12 +48,12 @@ fromList [-1.0,0.0,1.0]
46 48
47-} 49-}
48conv ker v 50conv ker v
49 | dim ker == 0 = constant 0 (dim v) 51 | dim ker == 0 = konst 0 (dim v)
50 | otherwise = corr ker' v' 52 | otherwise = corr ker' v'
51 where 53 where
52 ker' = (flatten.fliprl.asRow) ker 54 ker' = (flatten.fliprl.asRow) ker
53 v' = vjoin [z,v,z] 55 v' = vjoin [z,v,z]
54 z = constant 0 (dim ker -1) 56 z = konst 0 (dim ker -1)
55 57
56corrMin :: (Container Vector t, RealElement t, Product t) 58corrMin :: (Container Vector t, RealElement t, Product t)
57 => Vector t 59 => Vector t
diff --git a/packages/base/src/Numeric/Sparse.hs b/packages/base/src/Numeric/Sparse.hs
index 1957d3a..2df4578 100644
--- a/packages/base/src/Numeric/Sparse.hs
+++ b/packages/base/src/Numeric/Sparse.hs
@@ -10,7 +10,7 @@ module Numeric.Sparse(
10 smXv 10 smXv
11)where 11)where
12 12
13import Numeric.Container 13import Data.Packed.Numeric
14import qualified Data.Vector.Storable as V 14import qualified Data.Vector.Storable as V
15import Data.Function(on) 15import Data.Function(on)
16import Control.Arrow((***)) 16import Control.Arrow((***))
@@ -133,10 +133,6 @@ instance Transposable SMatrix
133 tr (CSC vs rs cs n m) = CSR vs rs cs m n 133 tr (CSC vs rs cs n m) = CSR vs rs cs m n
134 tr (Diag v n m) = Diag v m n 134 tr (Diag v n m) = Diag v m n
135 135
136instance Transposable (Matrix Double)
137 where
138 tr = trans
139
140 136
141instance CGMat SMatrix 137instance CGMat SMatrix
142instance CGMat (Matrix Double) 138instance CGMat (Matrix Double)
diff --git a/packages/gsl/src/Numeric/GSL/Random.hs b/packages/gsl/src/Numeric/GSL/Random.hs
index 07fa64d..f1f49e5 100644
--- a/packages/gsl/src/Numeric/GSL/Random.hs
+++ b/packages/gsl/src/Numeric/GSL/Random.hs
@@ -21,7 +21,7 @@ module Numeric.GSL.Random (
21) where 21) where
22 22
23import Numeric.GSL.Vector 23import Numeric.GSL.Vector
24import Numeric.LinearAlgebra(Seed,cholSH) 24import Numeric.LinearAlgebra(cholSH)
25import Numeric.Container hiding ( 25import Numeric.Container hiding (
26 randomVector, 26 randomVector,
27 gaussianSample, 27 gaussianSample,
@@ -30,7 +30,7 @@ import Numeric.Container hiding (
30import System.Random(randomIO) 30import System.Random(randomIO)
31 31
32 32
33 33type Seed = Int
34 34
35-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate 35-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
36-- Gaussian distribution. 36-- Gaussian distribution.
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 5fdc5af..3803f3b 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -24,8 +24,8 @@ module Numeric.LinearAlgebra.Tests(
24--, runBigTests 24--, runBigTests
25) where 25) where
26 26
27--import Data.Packed.Random
28import Numeric.LinearAlgebra 27import Numeric.LinearAlgebra
28import Numeric.HMatrix
29import Numeric.LinearAlgebra.Util(col,row) 29import Numeric.LinearAlgebra.Util(col,row)
30import Data.Packed 30import Data.Packed
31import Numeric.LinearAlgebra.LAPACK 31import Numeric.LinearAlgebra.LAPACK
@@ -56,8 +56,8 @@ import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
56 56
57import Test.QuickCheck.Test(isSuccess) 57import Test.QuickCheck.Test(isSuccess)
58 58
59eps = peps :: Double 59--eps = peps :: Double
60i = 0:+1 :: Complex Double 60--i = 0:+1 :: Complex Double
61 61
62qCheck n x = do 62qCheck n x = do
63 r <- quickCheckWithResult stdArgs {maxSize = n} x 63 r <- quickCheckWithResult stdArgs {maxSize = n} x
@@ -760,7 +760,7 @@ mkVecBench = do
760subBench = do 760subBench = do
761 putStrLn "" 761 putStrLn ""
762 let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v)) 762 let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v))
763 time "0.1M subVector " (g (constant 1 (1+10^5) :: Vector Double) @> 0) 763 time "0.1M subVector " (g (konst 1 (1+10^5) :: Vector Double) @> 0)
764 let f = foldl1' (.) (replicate (10^5) (fromRows.toRows)) 764 let f = foldl1' (.) (replicate (10^5) (fromRows.toRows))
765 time "subVector-join 3" (f (ident 3 :: Matrix Double) @@>(0,0)) 765 time "subVector-join 3" (f (ident 3 :: Matrix Double) @@>(0,0))
766 time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0)) 766 time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0))
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
index 4f5a5d2..b306f94 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -30,7 +30,7 @@ import System.Random
30 30
31import Numeric.LinearAlgebra 31import Numeric.LinearAlgebra
32import Numeric.LinearAlgebra.Devel 32import Numeric.LinearAlgebra.Devel
33import Numeric.Container hiding (Testable) 33import Numeric.Container
34import Control.Monad(replicateM) 34import Control.Monad(replicateM)
35import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector 35import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
36 ,sized,classify,Testable,Property 36 ,sized,classify,Testable,Property
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
index 423edaa..f4f1b89 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -43,7 +43,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
43 linearSolveProp, linearSolveProp2 43 linearSolveProp, linearSolveProp2
44) where 44) where
45 45
46import Numeric.Container hiding (Testable) 46import Numeric.Container
47import Numeric.LinearAlgebra --hiding (real,complex) 47import Numeric.LinearAlgebra --hiding (real,complex)
48import Numeric.LinearAlgebra.LAPACK 48import Numeric.LinearAlgebra.LAPACK
49import Debug.Trace 49import Debug.Trace