summaryrefslogtreecommitdiff
path: root/packages/base/src/Numeric
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Numeric')
-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
9 files changed, 351 insertions, 424 deletions
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)