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