summaryrefslogtreecommitdiff
path: root/packages/base/src
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src')
-rw-r--r--packages/base/src/Data/Packed/Numeric.hs4
-rw-r--r--packages/base/src/Numeric/Complex.hs (renamed from packages/base/src/Numeric/LinearAlgebra/Complex.hs)30
-rw-r--r--packages/base/src/Numeric/HMatrix.hs634
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Data.hs4
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/HMatrix.hs201
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Real.hs480
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Static.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Util.hs14
-rw-r--r--packages/base/src/Numeric/Sparse.hs2
9 files changed, 690 insertions, 681 deletions
diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs
index d2a20be..e59c1cd 100644
--- a/packages/base/src/Data/Packed/Numeric.hs
+++ b/packages/base/src/Data/Packed/Numeric.hs
@@ -106,7 +106,7 @@ infixr 8 <ยท>, #>
106 106
107{- | dot product 107{- | dot product
108 108
109>>> vect [1,2,3,4] <ยท> vect [-2,0,1,1] 109>>> vector [1,2,3,4] <ยท> vector [-2,0,1,1]
1105.0 1105.0
111 111
112>>> let ๐‘– = 0:+1 :: โ„‚ 112>>> let ๐‘– = 0:+1 :: โ„‚
@@ -128,7 +128,7 @@ infixr 8 <ยท>, #>
128 [ 1.0, 2.0, 3.0 128 [ 1.0, 2.0, 3.0
129 , 4.0, 5.0, 6.0 ] 129 , 4.0, 5.0, 6.0 ]
130 130
131>>> let v = vect [10,20,30] 131>>> let v = vector [10,20,30]
132 132
133>>> m #> v 133>>> m #> v
134fromList [140.0,320.0] 134fromList [140.0,320.0]
diff --git a/packages/base/src/Numeric/LinearAlgebra/Complex.hs b/packages/base/src/Numeric/Complex.hs
index 17bc397..d194af3 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Complex.hs
+++ b/packages/base/src/Numeric/Complex.hs
@@ -14,43 +14,35 @@
14 14
15 15
16{- | 16{- |
17Module : Numeric.LinearAlgebra.Complex 17Module : Numeric.HMatrix.Static.Complex
18Copyright : (c) Alberto Ruiz 2006-14 18Copyright : (c) Alberto Ruiz 2006-14
19License : BSD3 19License : BSD3
20Stability : experimental 20Stability : experimental
21 21
22-} 22-}
23 23
24module Numeric.LinearAlgebra.Complex( 24module Numeric.Complex(
25 C, 25 C, M,
26 vec2, vec3, vec4, (&), (#), 26 vec2, vec3, vec4, (&), (#),
27 vect, 27 vect,
28 R 28 Her, her, ๐‘–,
29) where 29) where
30 30
31import GHC.TypeLits 31import GHC.TypeLits
32import Numeric.HMatrix hiding ( 32import Numeric.LinearAlgebra.Util(โ„‚,iC)
33 (<>),(#>),(<ยท>),Konst(..),diag, disp,(ยฆ),(โ€”โ€”),row,col,vect,mat,linspace) 33import qualified Numeric.LinearAlgebra.HMatrix as LA
34import qualified Numeric.HMatrix as LA
35import Data.Proxy(Proxy)
36import Numeric.LinearAlgebra.Static 34import Numeric.LinearAlgebra.Static
37 35
38 36
37๐‘– :: Sized โ„‚ s c => s
38๐‘– = konst iC
39 39
40instance forall n . KnownNat n => Show (C n) 40newtype Her n = Her (M n n)
41 where
42 show (ud1 -> v)
43 | size v == 1 = "("++show (v!0)++" :: C "++show d++")"
44 | otherwise = "(vect"++ drop 8 (show v)++" :: C "++show d++")"
45 where
46 d = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
47 41
42her :: KnownNat n => M n n -> Her n
43her m = Her $ (m + LA.tr m)/2
48 44
49ud1 :: C n -> Vector โ„‚
50ud1 (C (Dim v)) = v
51 45
52mkC :: Vector โ„‚ -> C n
53mkC = C . Dim
54 46
55 47
56infixl 4 & 48infixl 4 &
diff --git a/packages/base/src/Numeric/HMatrix.hs b/packages/base/src/Numeric/HMatrix.hs
index ec96bfc..421333a 100644
--- a/packages/base/src/Numeric/HMatrix.hs
+++ b/packages/base/src/Numeric/HMatrix.hs
@@ -1,201 +1,497 @@
1----------------------------------------------------------------------------- 1{-# LANGUAGE DataKinds #-}
2{-# LANGUAGE KindSignatures #-}
3{-# LANGUAGE GeneralizedNewtypeDeriving #-}
4{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE FlexibleContexts #-}
7{-# LANGUAGE ScopedTypeVariables #-}
8{-# LANGUAGE EmptyDataDecls #-}
9{-# LANGUAGE Rank2Types #-}
10{-# LANGUAGE FlexibleInstances #-}
11{-# LANGUAGE TypeOperators #-}
12{-# LANGUAGE ViewPatterns #-}
13{-# LANGUAGE GADTs #-}
14{-# LANGUAGE OverlappingInstances #-}
15{-# LANGUAGE TypeFamilies #-}
16
17
2{- | 18{- |
3Module : Numeric.HMatrix 19Module : Numeric.HMatrix
4Copyright : (c) Alberto Ruiz 2006-14 20Copyright : (c) Alberto Ruiz 2006-14
5License : BSD3 21License : BSD3
6Maintainer : Alberto Ruiz 22Stability : experimental
7Stability : provisional 23
24Experimental interface for real arrays with statically checked dimensions.
8 25
9-} 26-}
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 -- >>> vect [1,2,3] * vect [3,0,-2]
21 -- fromList [3.0,0.0,-6.0]
22 --
23 -- >>> mat 3 [1..9] * ident 3
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 -- >>> mat 3 [1..9] + mat 1 [10,20,30]
40 -- (3><3)
41 -- [ 11.0, 12.0, 13.0
42 -- , 24.0, 25.0, 26.0
43 -- , 37.0, 38.0, 39.0 ]
44 --
45 27
28module Numeric.HMatrix(
29 -- * Vector
30 R,
31 vec2, vec3, vec4, (&), (#), split, headTail,
32 vector,
33 linspace, range, dim,
34 -- * Matrix
35 L, Sq, build,
36 row, col, (ยฆ),(โ€”โ€”), splitRows, splitCols,
37 unrow, uncol,
38
39 eye,
40 diagR, diag,
41 blockAt,
42 matrix,
46 -- * Products 43 -- * Products
47 -- ** dot 44 (<>),(#>),(<ยท>),
48 (<ยท>),
49 -- ** matrix-vector
50 (#>), (!#>),
51 -- ** matrix-matrix
52 (<>),
53 -- | The matrix x matrix product is also implemented in the "Data.Monoid" instance, where
54 -- single-element matrices (created from numeric literals or using 'scalar')
55 -- are used for scaling.
56 --
57 -- >>> import Data.Monoid as M
58 -- >>> let m = mat 3 [1..6]
59 -- >>> m M.<> 2 M.<> diagl[0.5,1,0]
60 -- (2><3)
61 -- [ 1.0, 4.0, 0.0
62 -- , 4.0, 10.0, 0.0 ]
63 --
64 -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
65
66
67 -- ** other
68 outer, kronecker, cross,
69 scale,
70 sumElements, prodElements,
71
72 -- * Linear Systems 45 -- * Linear Systems
73 (<\>), 46 linSolve, (<\>),
74 linearSolve, 47 -- * Factorizations
75 linearSolveLS, 48 svd, svdTall, svdFlat, Eigen(..),
76 linearSolveSVD, 49 withNullspace,
77 luSolve, 50 -- * Misc
78 cholSolve, 51 Disp(..),
79 cgSolve, 52 withVector, withMatrix,
80 cgSolve', 53 toRows, toColumns,
54 Sized(..), Diag(..), Sym, sym,
55 module Numeric.LinearAlgebra.HMatrix
56) where
81 57
82 -- * Inverse and pseudoinverse
83 inv, pinv, pinvTol,
84 58
85 -- * Determinant and rank 59import GHC.TypeLits
86 rcond, rank, 60import Numeric.LinearAlgebra.HMatrix hiding (
87 det, invlndet, 61 (<>),(#>),(<ยท>),Konst(..),diag, disp,(ยฆ),(โ€”โ€”),row,col,vector,matrix,linspace,toRows,toColumns,
62 (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',eigenvalues,eigenvaluesSH,eigenvaluesSH',build)
63import qualified Numeric.LinearAlgebra.HMatrix as LA
64import Data.Proxy(Proxy)
65import Numeric.LinearAlgebra.Static
66import Control.Arrow((***))
88 67
89 -- * Norms
90 Normed(..),
91 norm_Frob, norm_nuclear,
92 68
93 -- * Nullspace and range
94 orth,
95 nullspace, null1, null1sym,
96 69
97 -- * SVD
98 svd,
99 fullSVD,
100 thinSVD,
101 compactSVD,
102 singularValues,
103 leftSV, rightSV,
104 70
105 -- * Eigensystems
106 eig, eigSH, eigSH',
107 eigenvalues, eigenvaluesSH, eigenvaluesSH',
108 geigSH',
109 71
110 -- * QR 72ud1 :: R n -> Vector โ„
111 qr, rq, qrRaw, qrgr, 73ud1 (R (Dim v)) = v
112 74
113 -- * Cholesky
114 chol, cholSH, mbCholSH,
115 75
116 -- * Hessenberg 76infixl 4 &
117 hess, 77(&) :: forall n . KnownNat n
78 => R n -> โ„ -> R (n+1)
79u & x = u # (konst x :: R 1)
118 80
119 -- * Schur 81infixl 4 #
120 schur, 82(#) :: forall n m . (KnownNat n, KnownNat m)
83 => R n -> R m -> R (n+m)
84(R u) # (R v) = R (vconcat u v)
121 85
122 -- * LU
123 lu, luPacked,
124 86
125 -- * Matrix functions
126 expm,
127 sqrtm,
128 matFunc,
129 87
130 -- * Correlation and convolution 88vec2 :: โ„ -> โ„ -> R 2
131 corr, conv, corrMin, corr2, conv2, 89vec2 a b = R (gvec2 a b)
132 90
133 -- * Random arrays 91vec3 :: โ„ -> โ„ -> โ„ -> R 3
92vec3 a b c = R (gvec3 a b c)
93
94
95vec4 :: โ„ -> โ„ -> โ„ -> โ„ -> R 4
96vec4 a b c d = R (gvec4 a b c d)
97
98vector :: KnownNat n => [โ„] -> R n
99vector = fromList
100
101matrix :: (KnownNat m, KnownNat n) => [โ„] -> L m n
102matrix = fromList
103
104linspace :: forall n . KnownNat n => (โ„,โ„) -> R n
105linspace (a,b) = mkR (LA.linspace d (a,b))
106 where
107 d = fromIntegral . natVal $ (undefined :: Proxy n)
108
109range :: forall n . KnownNat n => R n
110range = mkR (LA.linspace d (1,fromIntegral d))
111 where
112 d = fromIntegral . natVal $ (undefined :: Proxy n)
113
114dim :: forall n . KnownNat n => R n
115dim = mkR (scalar d)
116 where
117 d = fromIntegral . natVal $ (undefined :: Proxy n)
118
119
120--------------------------------------------------------------------------------
121
122
123ud2 :: L m n -> Matrix โ„
124ud2 (L (Dim (Dim x))) = x
125
126
127--------------------------------------------------------------------------------
128--------------------------------------------------------------------------------
129
130diagR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => โ„ -> R k -> L m n
131diagR x v = mkL (asRow (vjoin [scalar x, ev, zeros]))
132 where
133 ev = extract v
134 zeros = LA.konst x (max 0 ((min m' n') - size ev))
135 m' = fromIntegral . natVal $ (undefined :: Proxy m)
136 n' = fromIntegral . natVal $ (undefined :: Proxy n)
137
138diag :: KnownNat n => R n -> Sq n
139diag = diagR 0
140
141eye :: KnownNat n => Sq n
142eye = diag 1
143
144--------------------------------------------------------------------------------
145
146blockAt :: forall m n . (KnownNat m, KnownNat n) => โ„ -> Int -> Int -> Matrix Double -> L m n
147blockAt x r c a = mkL res
148 where
149 z = scalar x
150 z1 = LA.konst x (r,c)
151 z2 = LA.konst x (max 0 (m'-(ra+r)), max 0 (n'-(ca+c)))
152 ra = min (rows a) . max 0 $ m'-r
153 ca = min (cols a) . max 0 $ n'-c
154 sa = subMatrix (0,0) (ra, ca) a
155 m' = fromIntegral . natVal $ (undefined :: Proxy m)
156 n' = fromIntegral . natVal $ (undefined :: Proxy n)
157 res = fromBlocks [[z1,z,z],[z,sa,z],[z,z,z2]]
158
159
160
161
162
163--------------------------------------------------------------------------------
164
165
166row :: R n -> L 1 n
167row = mkL . asRow . ud1
168
169--col :: R n -> L n 1
170col v = tr . row $ v
171
172unrow :: L 1 n -> R n
173unrow = mkR . head . LA.toRows . ud2
174
175--uncol :: L n 1 -> R n
176uncol v = unrow . tr $ v
177
178
179infixl 2 โ€”โ€”
180(โ€”โ€”) :: (KnownNat r1, KnownNat r2, KnownNat c) => L r1 c -> L r2 c -> L (r1+r2) c
181a โ€”โ€” b = mkL (extract a LA.โ€”โ€” extract b)
182
183
184infixl 3 ยฆ
185-- (ยฆ) :: (KnownNat r, KnownNat c1, KnownNat c2) => L r c1 -> L r c2 -> L r (c1+c2)
186a ยฆ b = tr (tr a โ€”โ€” tr b)
187
188
189type Sq n = L n n
190--type CSq n = CL n n
191
192type GL = (KnownNat n, KnownNat m) => L m n
193type GSq = KnownNat n => Sq n
194
195isKonst :: forall m n . (KnownNat m, KnownNat n) => L m n -> Maybe (โ„,(Int,Int))
196isKonst (unwrap -> x)
197 | singleM x = Just (x `atIndex` (0,0), (m',n'))
198 | otherwise = Nothing
199 where
200 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
201 n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
134 202
135 Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
136 203
137 -- * Misc
138 meanCov, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv,
139 โ„,โ„‚,iC,
140 -- * Auxiliary classes
141 Element, Container, Product, Numeric, LSDiv,
142 Complexable, RealElement,
143 RealOf, ComplexOf, SingleOf, DoubleOf,
144 IndexOf,
145 Field,
146-- Normed,
147 Transposable,
148 CGState(..),
149 Testable(..)
150) where
151 204
152import Numeric.LinearAlgebra.Data
153
154import Numeric.Matrix()
155import Numeric.Vector()
156import Data.Packed.Numeric hiding ((<>))
157import Numeric.LinearAlgebra.Algorithms hiding (linearSolve,Normed,orth)
158import qualified Numeric.LinearAlgebra.Algorithms as A
159import Numeric.LinearAlgebra.Util
160import Numeric.LinearAlgebra.Random
161import Numeric.Sparse((!#>))
162import Numeric.LinearAlgebra.Util.CG
163
164{- | dense matrix product
165
166>>> let a = (3><5) [1..]
167>>> a
168(3><5)
169 [ 1.0, 2.0, 3.0, 4.0, 5.0
170 , 6.0, 7.0, 8.0, 9.0, 10.0
171 , 11.0, 12.0, 13.0, 14.0, 15.0 ]
172
173>>> let b = (5><2) [1,3, 0,2, -1,5, 7,7, 6,0]
174>>> b
175(5><2)
176 [ 1.0, 3.0
177 , 0.0, 2.0
178 , -1.0, 5.0
179 , 7.0, 7.0
180 , 6.0, 0.0 ]
181
182>>> a <> b
183(3><2)
184 [ 56.0, 50.0
185 , 121.0, 135.0
186 , 186.0, 220.0 ]
187 205
188-}
189(<>) :: Numeric t => Matrix t -> Matrix t -> Matrix t
190(<>) = mXm
191infixr 8 <> 206infixr 8 <>
207(<>) :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => L m k -> L k n -> L m n
208
209(isKonst -> Just (a,(_,k))) <> (isKonst -> Just (b,_)) = konst (a * b * fromIntegral k)
210
211(isDiag -> Just (0,a,_)) <> (isDiag -> Just (0,b,_)) = diagR 0 (mkR v :: R k)
212 where
213 v = a' * b'
214 n = min (size a) (size b)
215 a' = subVector 0 n a
216 b' = subVector 0 n b
217
218(isDiag -> Just (0,a,_)) <> (extract -> b) = mkL (asColumn a * takeRows (size a) b)
219
220(extract -> a) <> (isDiag -> Just (0,b,_)) = mkL (takeColumns (size b) a * asRow b)
221
222a <> b = mkL (extract a LA.<> extract b)
223
224infixr 8 #>
225(#>) :: (KnownNat m, KnownNat n) => L m n -> R n -> R m
226(isDiag -> Just (0, w, _)) #> v = mkR (w * subVector 0 (size w) (extract v))
227m #> v = mkR (extract m LA.#> extract v)
228
229
230infixr 8 <ยท>
231(<ยท>) :: R n -> R n -> โ„
232(ud1 -> u) <ยท> (ud1 -> v)
233 | singleV u || singleV v = sumElements (u * v)
234 | otherwise = udot u v
235
236--------------------------------------------------------------------------------
237
238{-
239class Minim (n :: Nat) (m :: Nat)
240 where
241 type Mini n m :: Nat
242
243instance forall (n :: Nat) . Minim n n
244 where
245 type Mini n n = n
246
247
248instance forall (n :: Nat) (m :: Nat) . (n <= m+1) => Minim n m
249 where
250 type Mini n m = n
251
252instance forall (n :: Nat) (m :: Nat) . (m <= n+1) => Minim n m
253 where
254 type Mini n m = m
255-}
256
257class Diag m d | m -> d
258 where
259 takeDiag :: m -> d
260
261
262
263instance forall n . (KnownNat n) => Diag (L n n) (R n)
264 where
265 takeDiag m = mkR (LA.takeDiag (extract m))
266
267
268instance forall m n . (KnownNat m, KnownNat n, m <= n+1) => Diag (L m n) (R m)
269 where
270 takeDiag m = mkR (LA.takeDiag (extract m))
271
272
273instance forall m n . (KnownNat m, KnownNat n, n <= m+1) => Diag (L m n) (R n)
274 where
275 takeDiag m = mkR (LA.takeDiag (extract m))
276
277
278--------------------------------------------------------------------------------
279
280linSolve :: (KnownNat m, KnownNat n) => L m m -> L m n -> Maybe (L m n)
281linSolve (extract -> a) (extract -> b) = fmap mkL (LA.linearSolve a b)
282
283(<\>) :: (KnownNat m, KnownNat n, KnownNat r) => L m n -> L m r -> L n r
284(extract -> a) <\> (extract -> b) = mkL (a LA.<\> b)
285
286svd :: (KnownNat m, KnownNat n) => L m n -> (L m m, R n, L n n)
287svd (extract -> m) = (mkL u, mkR s', mkL v)
288 where
289 (u,s,v) = LA.svd m
290 s' = vjoin [s, z]
291 z = LA.konst 0 (max 0 (cols m - size s))
292
293
294svdTall :: (KnownNat m, KnownNat n, n <= m) => L m n -> (L m n, R n, L n n)
295svdTall (extract -> m) = (mkL u, mkR s, mkL v)
296 where
297 (u,s,v) = LA.thinSVD m
298
299
300svdFlat :: (KnownNat m, KnownNat n, m <= n) => L m n -> (L m m, R m, L n m)
301svdFlat (extract -> m) = (mkL u, mkR s, mkL v)
302 where
303 (u,s,v) = LA.thinSVD m
304
305--------------------------------------------------------------------------------
306
307class Eigen m l v | m -> l, m -> v
308 where
309 eigensystem :: m -> (l,v)
310 eigenvalues :: m -> l
311
312newtype Sym n = Sym (Sq n) deriving Show
313
314
315sym :: KnownNat n => Sq n -> Sym n
316sym m = Sym $ (m + tr m)/2
317
318
319
320instance KnownNat n => Eigen (Sym n) (R n) (L n n)
321 where
322 eigenvalues (Sym (extract -> m)) = mkR . LA.eigenvaluesSH' $ m
323 eigensystem (Sym (extract -> m)) = (mkR l, mkL v)
324 where
325 (l,v) = LA.eigSH' m
326
327instance KnownNat n => Eigen (Sq n) (C n) (M n n)
328 where
329 eigenvalues (extract -> m) = mkC . LA.eigenvalues $ m
330 eigensystem (extract -> m) = (mkC l, mkM v)
331 where
332 (l,v) = LA.eig m
333
334--------------------------------------------------------------------------------
335
336
337withNullspace
338 :: forall m n z . (KnownNat m, KnownNat n)
339 => L m n
340 -> (forall k . (KnownNat k) => L n k -> z)
341 -> z
342withNullspace (LA.nullspace . extract -> a) f =
343 case someNatVal $ fromIntegral $ cols a of
344 Nothing -> error "static/dynamic mismatch"
345 Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k)
346
347--------------------------------------------------------------------------------
348
349split :: forall p n . (KnownNat p, KnownNat n, p<=n) => R n -> (R p, R (n-p))
350split (extract -> v) = ( mkR (subVector 0 p' v) ,
351 mkR (subVector p' (size v - p') v) )
352 where
353 p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int
354
355
356headTail :: (KnownNat n, 1<=n) => R n -> (โ„, R (n-1))
357headTail = ((!0) . extract *** id) . split
358
359
360splitRows :: forall p m n . (KnownNat p, KnownNat m, KnownNat n, p<=m) => L m n -> (L p n, L (m-p) n)
361splitRows (extract -> x) = ( mkL (takeRows p' x) ,
362 mkL (dropRows p' x) )
363 where
364 p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int
365
366splitCols :: forall p m n. (KnownNat p, KnownNat m, KnownNat n, KnownNat (n-p), p<=n) => L m n -> (L m p, L m (n-p))
367splitCols = (tr *** tr) . splitRows . tr
368
369
370toRows :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R n]
371toRows (LA.toRows . extract -> vs) = map mkR vs
372
373
374toColumns :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R m]
375toColumns (LA.toColumns . extract -> vs) = map mkR vs
376
377
378splittest
379 = do
380 let v = range :: R 7
381 a = snd (split v) :: R 4
382 print $ a
383 print $ snd . headTail . snd . headTail $ v
384 print $ first (vec3 1 2 3)
385 print $ second (vec3 1 2 3)
386 print $ third (vec3 1 2 3)
387 print $ (snd $ splitRows eye :: L 4 6)
388 where
389 first v = fst . headTail $ v
390 second v = first . snd . headTail $ v
391 third v = first . snd . headTail . snd . headTail $ v
392
393--------------------------------------------------------------------------------
394
395build
396 :: forall m n . (KnownNat n, KnownNat m)
397 => (โ„ -> โ„ -> โ„)
398 -> L m n
399build f = mkL $ LA.build (m',n') f
400 where
401 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
402 n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
192 403
193-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. 404--------------------------------------------------------------------------------
194linearSolve m b = A.mbLinearSolve m b
195 405
196-- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'. 406withVector
197nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m) 407 :: forall z
408 . Vector โ„
409 -> (forall n . (KnownNat n) => R n -> z)
410 -> z
411withVector v f =
412 case someNatVal $ fromIntegral $ size v of
413 Nothing -> error "static/dynamic mismatch"
414 Just (SomeNat (_ :: Proxy m)) -> f (mkR v :: R m)
415
416
417withMatrix
418 :: forall z
419 . Matrix โ„
420 -> (forall m n . (KnownNat m, KnownNat n) => L m n -> z)
421 -> z
422withMatrix a f =
423 case someNatVal $ fromIntegral $ rows a of
424 Nothing -> error "static/dynamic mismatch"
425 Just (SomeNat (_ :: Proxy m)) ->
426 case someNatVal $ fromIntegral $ cols a of
427 Nothing -> error "static/dynamic mismatch"
428 Just (SomeNat (_ :: Proxy n)) ->
429 f (mkL a :: L m n)
430
431
432--------------------------------------------------------------------------------
433
434test :: (Bool, IO ())
435test = (ok,info)
436 where
437 ok = extract (eye :: Sq 5) == ident 5
438 && unwrap (mTm sm :: Sq 3) == tr ((3><3)[1..]) LA.<> (3><3)[1..]
439 && unwrap (tm :: L 3 5) == LA.matrix 5 [1..15]
440 && thingS == thingD
441 && precS == precD
442 && withVector (LA.vector [1..15]) sumV == sumElements (LA.fromList [1..15])
443
444 info = do
445 print $ u
446 print $ v
447 print (eye :: Sq 3)
448 print $ ((u & 5) + 1) <ยท> v
449 print (tm :: L 2 5)
450 print (tm <> sm :: L 2 3)
451 print thingS
452 print thingD
453 print precS
454 print precD
455 print $ withVector (LA.vector [1..15]) sumV
456 splittest
457
458 sumV w = w <ยท> konst 1
459
460 u = vec2 3 5
461
462 ๐•ง x = vector [x] :: R 1
463
464 v = ๐•ง 2 & 4 & 7
465
466-- mTm :: L n m -> Sq m
467 mTm a = tr a <> a
468
469 tm :: GL
470 tm = lmat 0 [1..]
471
472 lmat :: forall m n . (KnownNat m, KnownNat n) => โ„ -> [โ„] -> L m n
473 lmat z xs = mkL . reshape n' . LA.fromList . take (m'*n') $ xs ++ repeat z
474 where
475 m' = fromIntegral . natVal $ (undefined :: Proxy m)
476 n' = fromIntegral . natVal $ (undefined :: Proxy n)
477
478 sm :: GSq
479 sm = lmat 0 [1..]
480
481 thingS = (u & 1) <ยท> tr q #> q #> v
482 where
483 q = tm :: L 10 3
484
485 thingD = vjoin [ud1 u, 1] LA.<ยท> tr m LA.#> m LA.#> ud1 v
486 where
487 m = LA.matrix 3 [1..30]
488
489 precS = (1::Double) + (2::Double) * ((1 :: R 3) * (u & 6)) <ยท> konst 2 #> v
490 precD = 1 + 2 * vjoin[ud1 u, 6] LA.<ยท> LA.konst 2 (size (ud1 u) +1, size (ud1 v)) LA.#> ud1 v
491
492
493instance (KnownNat n', KnownNat m') => Testable (L n' m')
494 where
495 checkT _ = test
198 496
199-- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'.
200orth m = orthSVD (Left (1*eps)) m (leftSV m)
201 497
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index 20f3b81..4e4868a 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -16,11 +16,11 @@ 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 vect, (|>), 19 vector, (|>),
20 20
21 -- * Matrix 21 -- * Matrix
22 22
23 mat, (><), tr, 23 matrix, (><), tr,
24 24
25 -- * Indexing 25 -- * Indexing
26 26
diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
new file mode 100644
index 0000000..54ddd68
--- /dev/null
+++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
@@ -0,0 +1,201 @@
1-----------------------------------------------------------------------------
2{- |
3Module : Numeric.LinearAlgebra.HMatrix
4Copyright : (c) Alberto Ruiz 2006-14
5License : BSD3
6Maintainer : Alberto Ruiz
7Stability : provisional
8
9-}
10-----------------------------------------------------------------------------
11module Numeric.LinearAlgebra.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 -- >>> vector [1,2,3] * vector [3,0,-2]
21 -- fromList [3.0,0.0,-6.0]
22 --
23 -- >>> matrix 3 [1..9] * ident 3
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 -- >>> matrix 3 [1..9] + matrix 1 [10,20,30]
40 -- (3><3)
41 -- [ 11.0, 12.0, 13.0
42 -- , 24.0, 25.0, 26.0
43 -- , 37.0, 38.0, 39.0 ]
44 --
45
46 -- * Products
47 -- ** dot
48 (<ยท>),
49 -- ** matrix-vector
50 (#>), (!#>),
51 -- ** matrix-matrix
52 (<>),
53 -- | The matrix x matrix product is also implemented in the "Data.Monoid" instance, where
54 -- single-element matrices (created from numeric literals or using 'scalar')
55 -- are used for scaling.
56 --
57 -- >>> import Data.Monoid as M
58 -- >>> let m = matrix 3 [1..6]
59 -- >>> m M.<> 2 M.<> diagl[0.5,1,0]
60 -- (2><3)
61 -- [ 1.0, 4.0, 0.0
62 -- , 4.0, 10.0, 0.0 ]
63 --
64 -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
65
66
67 -- ** other
68 outer, kronecker, cross,
69 scale,
70 sumElements, prodElements,
71
72 -- * Linear Systems
73 (<\>),
74 linearSolve,
75 linearSolveLS,
76 linearSolveSVD,
77 luSolve,
78 cholSolve,
79 cgSolve,
80 cgSolve',
81
82 -- * Inverse and pseudoinverse
83 inv, pinv, pinvTol,
84
85 -- * Determinant and rank
86 rcond, rank,
87 det, invlndet,
88
89 -- * Norms
90 Normed(..),
91 norm_Frob, norm_nuclear,
92
93 -- * Nullspace and range
94 orth,
95 nullspace, null1, null1sym,
96
97 -- * SVD
98 svd,
99 fullSVD,
100 thinSVD,
101 compactSVD,
102 singularValues,
103 leftSV, rightSV,
104
105 -- * Eigensystems
106 eig, eigSH, eigSH',
107 eigenvalues, eigenvaluesSH, eigenvaluesSH',
108 geigSH',
109
110 -- * QR
111 qr, rq, qrRaw, qrgr,
112
113 -- * Cholesky
114 chol, cholSH, mbCholSH,
115
116 -- * Hessenberg
117 hess,
118
119 -- * Schur
120 schur,
121
122 -- * LU
123 lu, luPacked,
124
125 -- * Matrix functions
126 expm,
127 sqrtm,
128 matFunc,
129
130 -- * Correlation and convolution
131 corr, conv, corrMin, corr2, conv2,
132
133 -- * Random arrays
134
135 Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
136
137 -- * Misc
138 meanCov, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv,
139 โ„,โ„‚,iC,
140 -- * Auxiliary classes
141 Element, Container, Product, Numeric, LSDiv,
142 Complexable, RealElement,
143 RealOf, ComplexOf, SingleOf, DoubleOf,
144 IndexOf,
145 Field,
146-- Normed,
147 Transposable,
148 CGState(..),
149 Testable(..)
150) where
151
152import Numeric.LinearAlgebra.Data
153
154import Numeric.Matrix()
155import Numeric.Vector()
156import Data.Packed.Numeric hiding ((<>))
157import Numeric.LinearAlgebra.Algorithms hiding (linearSolve,Normed,orth)
158import qualified Numeric.LinearAlgebra.Algorithms as A
159import Numeric.LinearAlgebra.Util
160import Numeric.LinearAlgebra.Random
161import Numeric.Sparse((!#>))
162import Numeric.LinearAlgebra.Util.CG
163
164{- | dense matrix product
165
166>>> let a = (3><5) [1..]
167>>> a
168(3><5)
169 [ 1.0, 2.0, 3.0, 4.0, 5.0
170 , 6.0, 7.0, 8.0, 9.0, 10.0
171 , 11.0, 12.0, 13.0, 14.0, 15.0 ]
172
173>>> let b = (5><2) [1,3, 0,2, -1,5, 7,7, 6,0]
174>>> b
175(5><2)
176 [ 1.0, 3.0
177 , 0.0, 2.0
178 , -1.0, 5.0
179 , 7.0, 7.0
180 , 6.0, 0.0 ]
181
182>>> a <> b
183(3><2)
184 [ 56.0, 50.0
185 , 121.0, 135.0
186 , 186.0, 220.0 ]
187
188-}
189(<>) :: Numeric t => Matrix t -> Matrix t -> Matrix t
190(<>) = mXm
191infixr 8 <>
192
193-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'.
194linearSolve m b = A.mbLinearSolve m b
195
196-- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'.
197nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m)
198
199-- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'.
200orth m = orthSVD (Left (1*eps)) m (leftSV m)
201
diff --git a/packages/base/src/Numeric/LinearAlgebra/Real.hs b/packages/base/src/Numeric/LinearAlgebra/Real.hs
deleted file mode 100644
index 97c462e..0000000
--- a/packages/base/src/Numeric/LinearAlgebra/Real.hs
+++ /dev/null
@@ -1,480 +0,0 @@
1{-# LANGUAGE DataKinds #-}
2{-# LANGUAGE KindSignatures #-}
3{-# LANGUAGE GeneralizedNewtypeDeriving #-}
4{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE FlexibleContexts #-}
7{-# LANGUAGE ScopedTypeVariables #-}
8{-# LANGUAGE EmptyDataDecls #-}
9{-# LANGUAGE Rank2Types #-}
10{-# LANGUAGE FlexibleInstances #-}
11{-# LANGUAGE TypeOperators #-}
12{-# LANGUAGE ViewPatterns #-}
13{-# LANGUAGE GADTs #-}
14{-# LANGUAGE OverlappingInstances #-}
15{-# LANGUAGE TypeFamilies #-}
16
17
18{- |
19Module : Numeric.LinearAlgebra.Real
20Copyright : (c) Alberto Ruiz 2006-14
21License : BSD3
22Stability : experimental
23
24Experimental interface for real arrays with statically checked dimensions.
25
26-}
27
28module Numeric.LinearAlgebra.Real(
29 -- * Vector
30 R, C,
31 vec2, vec3, vec4, (&), (#), split, headTail,
32 vector,
33 linspace, range, dim,
34 -- * Matrix
35 L, Sq, M, def,
36 row, col, (ยฆ),(โ€”โ€”), splitRows, splitCols,
37 unrow, uncol,
38
39 eye,
40 diagR, diag,
41 blockAt,
42 matrix,
43 -- * Products
44 (<>),(#>),(<ยท>),
45 -- * Linear Systems
46 linSolve, (<\>),
47 -- * Factorizations
48 svd, svdTall, svdFlat, Eigen(..),
49 -- * Pretty printing
50 Disp(..),
51 -- * Misc
52 withVector, withMatrix,
53 Sized(..), Diag(..), Sym, sym, Her, her, ๐‘–,
54 module Numeric.HMatrix
55) where
56
57
58import GHC.TypeLits
59import Numeric.HMatrix hiding (
60 (<>),(#>),(<ยท>),Konst(..),diag, disp,(ยฆ),(โ€”โ€”),row,col,vect,mat,linspace,
61 (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',eigenvalues,eigenvaluesSH,eigenvaluesSH',build)
62import qualified Numeric.HMatrix as LA
63import Data.Proxy(Proxy)
64import Numeric.LinearAlgebra.Static
65import Control.Arrow((***))
66
67
68๐‘– :: Sized โ„‚ s c => s
69๐‘– = konst iC
70
71
72
73
74
75ud1 :: R n -> Vector โ„
76ud1 (R (Dim v)) = v
77
78
79infixl 4 &
80(&) :: forall n . KnownNat n
81 => R n -> โ„ -> R (n+1)
82u & x = u # (konst x :: R 1)
83
84infixl 4 #
85(#) :: forall n m . (KnownNat n, KnownNat m)
86 => R n -> R m -> R (n+m)
87(R u) # (R v) = R (vconcat u v)
88
89
90
91vec2 :: โ„ -> โ„ -> R 2
92vec2 a b = R (gvec2 a b)
93
94vec3 :: โ„ -> โ„ -> โ„ -> R 3
95vec3 a b c = R (gvec3 a b c)
96
97
98vec4 :: โ„ -> โ„ -> โ„ -> โ„ -> R 4
99vec4 a b c d = R (gvec4 a b c d)
100
101vector :: KnownNat n => [โ„] -> R n
102vector = fromList
103
104matrix :: (KnownNat m, KnownNat n) => [โ„] -> L m n
105matrix = fromList
106
107linspace :: forall n . KnownNat n => (โ„,โ„) -> R n
108linspace (a,b) = mkR (LA.linspace d (a,b))
109 where
110 d = fromIntegral . natVal $ (undefined :: Proxy n)
111
112range :: forall n . KnownNat n => R n
113range = mkR (LA.linspace d (1,fromIntegral d))
114 where
115 d = fromIntegral . natVal $ (undefined :: Proxy n)
116
117dim :: forall n . KnownNat n => R n
118dim = mkR (scalar d)
119 where
120 d = fromIntegral . natVal $ (undefined :: Proxy n)
121
122
123--------------------------------------------------------------------------------
124
125
126ud2 :: L m n -> Matrix โ„
127ud2 (L (Dim (Dim x))) = x
128
129
130--------------------------------------------------------------------------------
131--------------------------------------------------------------------------------
132
133diagR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => โ„ -> R k -> L m n
134diagR x v = mkL (asRow (vjoin [scalar x, ev, zeros]))
135 where
136 ev = extract v
137 zeros = LA.konst x (max 0 ((min m' n') - size ev))
138 m' = fromIntegral . natVal $ (undefined :: Proxy m)
139 n' = fromIntegral . natVal $ (undefined :: Proxy n)
140
141diag :: KnownNat n => R n -> Sq n
142diag = diagR 0
143
144eye :: KnownNat n => Sq n
145eye = diag 1
146
147--------------------------------------------------------------------------------
148
149blockAt :: forall m n . (KnownNat m, KnownNat n) => โ„ -> Int -> Int -> Matrix Double -> L m n
150blockAt x r c a = mkL res
151 where
152 z = scalar x
153 z1 = LA.konst x (r,c)
154 z2 = LA.konst x (max 0 (m'-(ra+r)), max 0 (n'-(ca+c)))
155 ra = min (rows a) . max 0 $ m'-r
156 ca = min (cols a) . max 0 $ n'-c
157 sa = subMatrix (0,0) (ra, ca) a
158 m' = fromIntegral . natVal $ (undefined :: Proxy m)
159 n' = fromIntegral . natVal $ (undefined :: Proxy n)
160 res = fromBlocks [[z1,z,z],[z,sa,z],[z,z,z2]]
161
162
163
164
165
166--------------------------------------------------------------------------------
167
168
169row :: R n -> L 1 n
170row = mkL . asRow . ud1
171
172--col :: R n -> L n 1
173col v = tr . row $ v
174
175unrow :: L 1 n -> R n
176unrow = mkR . head . toRows . ud2
177
178--uncol :: L n 1 -> R n
179uncol v = unrow . tr $ v
180
181
182infixl 2 โ€”โ€”
183(โ€”โ€”) :: (KnownNat r1, KnownNat r2, KnownNat c) => L r1 c -> L r2 c -> L (r1+r2) c
184a โ€”โ€” b = mkL (extract a LA.โ€”โ€” extract b)
185
186
187infixl 3 ยฆ
188-- (ยฆ) :: (KnownNat r, KnownNat c1, KnownNat c2) => L r c1 -> L r c2 -> L r (c1+c2)
189a ยฆ b = tr (tr a โ€”โ€” tr b)
190
191
192type Sq n = L n n
193--type CSq n = CL n n
194
195type GL = (KnownNat n, KnownNat m) => L m n
196type GSq = KnownNat n => Sq n
197
198isKonst :: forall m n . (KnownNat m, KnownNat n) => L m n -> Maybe (โ„,(Int,Int))
199isKonst (unwrap -> x)
200 | singleM x = Just (x `atIndex` (0,0), (m',n'))
201 | otherwise = Nothing
202 where
203 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
204 n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
205
206
207
208
209infixr 8 <>
210(<>) :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => L m k -> L k n -> L m n
211
212(isKonst -> Just (a,(_,k))) <> (isKonst -> Just (b,_)) = konst (a * b * fromIntegral k)
213
214(isDiag -> Just (0,a,_)) <> (isDiag -> Just (0,b,_)) = diagR 0 (mkR v :: R k)
215 where
216 v = a' * b'
217 n = min (size a) (size b)
218 a' = subVector 0 n a
219 b' = subVector 0 n b
220
221(isDiag -> Just (0,a,_)) <> (extract -> b) = mkL (asColumn a * takeRows (size a) b)
222
223(extract -> a) <> (isDiag -> Just (0,b,_)) = mkL (takeColumns (size b) a * asRow b)
224
225a <> b = mkL (extract a LA.<> extract b)
226
227infixr 8 #>
228(#>) :: (KnownNat m, KnownNat n) => L m n -> R n -> R m
229(isDiag -> Just (0, w, _)) #> v = mkR (w * subVector 0 (size w) (extract v))
230m #> v = mkR (extract m LA.#> extract v)
231
232
233infixr 8 <ยท>
234(<ยท>) :: R n -> R n -> โ„
235(ud1 -> u) <ยท> (ud1 -> v)
236 | singleV u || singleV v = sumElements (u * v)
237 | otherwise = udot u v
238
239--------------------------------------------------------------------------------
240
241{-
242class Minim (n :: Nat) (m :: Nat)
243 where
244 type Mini n m :: Nat
245
246instance forall (n :: Nat) . Minim n n
247 where
248 type Mini n n = n
249
250
251instance forall (n :: Nat) (m :: Nat) . (n <= m+1) => Minim n m
252 where
253 type Mini n m = n
254
255instance forall (n :: Nat) (m :: Nat) . (m <= n+1) => Minim n m
256 where
257 type Mini n m = m
258-}
259
260class Diag m d | m -> d
261 where
262 takeDiag :: m -> d
263
264
265
266instance forall n . (KnownNat n) => Diag (L n n) (R n)
267 where
268 takeDiag m = mkR (LA.takeDiag (extract m))
269
270
271instance forall m n . (KnownNat m, KnownNat n, m <= n+1) => Diag (L m n) (R m)
272 where
273 takeDiag m = mkR (LA.takeDiag (extract m))
274
275
276instance forall m n . (KnownNat m, KnownNat n, n <= m+1) => Diag (L m n) (R n)
277 where
278 takeDiag m = mkR (LA.takeDiag (extract m))
279
280
281--------------------------------------------------------------------------------
282
283linSolve :: (KnownNat m, KnownNat n) => L m m -> L m n -> Maybe (L m n)
284linSolve (extract -> a) (extract -> b) = fmap mkL (LA.linearSolve a b)
285
286(<\>) :: (KnownNat m, KnownNat n, KnownNat r) => L m n -> L m r -> L n r
287(extract -> a) <\> (extract -> b) = mkL (a LA.<\> b)
288
289svd :: (KnownNat m, KnownNat n) => L m n -> (L m m, R n, L n n)
290svd (extract -> m) = (mkL u, mkR s', mkL v)
291 where
292 (u,s,v) = LA.svd m
293 s' = vjoin [s, z]
294 z = LA.konst 0 (max 0 (cols m - size s))
295
296
297svdTall :: (KnownNat m, KnownNat n, n <= m) => L m n -> (L m n, R n, L n n)
298svdTall (extract -> m) = (mkL u, mkR s, mkL v)
299 where
300 (u,s,v) = LA.thinSVD m
301
302
303svdFlat :: (KnownNat m, KnownNat n, m <= n) => L m n -> (L m m, R m, L m n)
304svdFlat (extract -> m) = (mkL u, mkR s, mkL v)
305 where
306 (u,s,v) = LA.thinSVD m
307
308--------------------------------------------------------------------------------
309
310class Eigen m l v | m -> l, m -> v
311 where
312 eigensystem :: m -> (l,v)
313 eigenvalues :: m -> l
314
315newtype Sym n = Sym (Sq n) deriving Show
316
317newtype Her n = Her (M n n)
318
319sym :: KnownNat n => Sq n -> Sym n
320sym m = Sym $ (m + tr m)/2
321
322her :: KnownNat n => M n n -> Her n
323her m = Her $ (m + tr m)/2
324
325instance KnownNat n => Eigen (Sym n) (R n) (L n n)
326 where
327 eigenvalues (Sym (extract -> m)) = mkR . LA.eigenvaluesSH' $ m
328 eigensystem (Sym (extract -> m)) = (mkR l, mkL v)
329 where
330 (l,v) = LA.eigSH' m
331
332instance KnownNat n => Eigen (Sq n) (C n) (M n n)
333 where
334 eigenvalues (extract -> m) = mkC . LA.eigenvalues $ m
335 eigensystem (extract -> m) = (mkC l, mkM v)
336 where
337 (l,v) = LA.eig m
338
339--------------------------------------------------------------------------------
340
341split :: forall p n . (KnownNat p, KnownNat n, p<=n) => R n -> (R p, R (n-p))
342split (extract -> v) = ( mkR (subVector 0 p' v) ,
343 mkR (subVector p' (size v - p') v) )
344 where
345 p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int
346
347
348headTail :: (KnownNat n, 1<=n) => R n -> (โ„, R (n-1))
349headTail = ((!0) . extract *** id) . split
350
351
352splitRows :: forall p m n. (KnownNat p, KnownNat m, KnownNat n, p<=m) => L m n -> (L p n, L (m-p) n)
353splitRows (extract -> x) = ( mkL (takeRows p' x) ,
354 mkL (dropRows p' x) )
355 where
356 p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int
357
358splitCols :: forall p m n. (KnownNat p, KnownNat m, KnownNat n, KnownNat (n-p), p<=n) => L m n -> (L m p, L m (n-p))
359splitCols = (tr *** tr) . splitRows . tr
360
361
362splittest
363 = do
364 let v = range :: R 7
365 a = snd (split v) :: R 4
366 print $ a
367 print $ snd . headTail . snd . headTail $ v
368 print $ first (vec3 1 2 3)
369 print $ second (vec3 1 2 3)
370 print $ third (vec3 1 2 3)
371 print $ (snd $ splitRows eye :: L 4 6)
372 where
373 first v = fst . headTail $ v
374 second v = first . snd . headTail $ v
375 third v = first . snd . headTail . snd . headTail $ v
376
377--------------------------------------------------------------------------------
378
379def
380 :: forall m n . (KnownNat n, KnownNat m)
381 => (โ„ -> โ„ -> โ„)
382 -> L m n
383def f = mkL $ LA.build (m',n') f
384 where
385 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
386 n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int
387
388--------------------------------------------------------------------------------
389
390withVector
391 :: forall z
392 . Vector โ„
393 -> (forall n . (KnownNat n) => R n -> z)
394 -> z
395withVector v f =
396 case someNatVal $ fromIntegral $ size v of
397 Nothing -> error "static/dynamic mismatch"
398 Just (SomeNat (_ :: Proxy m)) -> f (mkR v :: R m)
399
400
401withMatrix
402 :: forall z
403 . Matrix โ„
404 -> (forall m n . (KnownNat m, KnownNat n) => L m n -> z)
405 -> z
406withMatrix a f =
407 case someNatVal $ fromIntegral $ rows a of
408 Nothing -> error "static/dynamic mismatch"
409 Just (SomeNat (_ :: Proxy m)) ->
410 case someNatVal $ fromIntegral $ cols a of
411 Nothing -> error "static/dynamic mismatch"
412 Just (SomeNat (_ :: Proxy n)) ->
413 f (mkL a :: L n m)
414
415--------------------------------------------------------------------------------
416
417test :: (Bool, IO ())
418test = (ok,info)
419 where
420 ok = extract (eye :: Sq 5) == ident 5
421 && unwrap (mTm sm :: Sq 3) == tr ((3><3)[1..]) LA.<> (3><3)[1..]
422 && unwrap (tm :: L 3 5) == LA.mat 5 [1..15]
423 && thingS == thingD
424 && precS == precD
425 && withVector (LA.vect [1..15]) sumV == sumElements (LA.fromList [1..15])
426
427 info = do
428 print $ u
429 print $ v
430 print (eye :: Sq 3)
431 print $ ((u & 5) + 1) <ยท> v
432 print (tm :: L 2 5)
433 print (tm <> sm :: L 2 3)
434 print thingS
435 print thingD
436 print precS
437 print precD
438 print $ withVector (LA.vect [1..15]) sumV
439 splittest
440
441 sumV w = w <ยท> konst 1
442
443 u = vec2 3 5
444
445 ๐•ง x = vector [x] :: R 1
446
447 v = ๐•ง 2 & 4 & 7
448
449-- mTm :: L n m -> Sq m
450 mTm a = tr a <> a
451
452 tm :: GL
453 tm = lmat 0 [1..]
454
455 lmat :: forall m n . (KnownNat m, KnownNat n) => โ„ -> [โ„] -> L m n
456 lmat z xs = mkL . reshape n' . LA.fromList . take (m'*n') $ xs ++ repeat z
457 where
458 m' = fromIntegral . natVal $ (undefined :: Proxy m)
459 n' = fromIntegral . natVal $ (undefined :: Proxy n)
460
461 sm :: GSq
462 sm = lmat 0 [1..]
463
464 thingS = (u & 1) <ยท> tr q #> q #> v
465 where
466 q = tm :: L 10 3
467
468 thingD = vjoin [ud1 u, 1] LA.<ยท> tr m LA.#> m LA.#> ud1 v
469 where
470 m = LA.mat 3 [1..30]
471
472 precS = (1::Double) + (2::Double) * ((1 :: R 3) * (u & 6)) <ยท> konst 2 #> v
473 precD = 1 + 2 * vjoin[ud1 u, 6] LA.<ยท> LA.konst 2 (size (ud1 u) +1, size (ud1 v)) LA.#> ud1 v
474
475
476instance (KnownNat n', KnownNat m') => Testable (L n' m')
477 where
478 checkT _ = test
479
480
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
index 2647f20..daf8d80 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs
@@ -25,7 +25,7 @@ module Numeric.LinearAlgebra.Static where
25 25
26 26
27import GHC.TypeLits 27import GHC.TypeLits
28import Numeric.HMatrix as LA 28import Numeric.LinearAlgebra.HMatrix as LA
29import Data.Packed as D 29import Data.Packed as D
30import Data.Packed.ST 30import Data.Packed.ST
31import Data.Proxy(Proxy) 31import Data.Proxy(Proxy)
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs
index 4824af4..0ac4634 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Util.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs
@@ -19,7 +19,7 @@ Stability : provisional
19module Numeric.LinearAlgebra.Util( 19module Numeric.LinearAlgebra.Util(
20 20
21 -- * Convenience functions 21 -- * Convenience functions
22 vect, mat, 22 vector, matrix,
23 disp, 23 disp,
24 zeros, ones, 24 zeros, ones,
25 diagl, 25 diagl,
@@ -79,27 +79,27 @@ iC = 0:+1
79 79
80{- | create a real vector 80{- | create a real vector
81 81
82>>> vect [1..5] 82>>> vector [1..5]
83fromList [1.0,2.0,3.0,4.0,5.0] 83fromList [1.0,2.0,3.0,4.0,5.0]
84 84
85-} 85-}
86vect :: [โ„] -> Vector โ„ 86vector :: [โ„] -> Vector โ„
87vect = fromList 87vector = fromList
88 88
89{- | create a real matrix 89{- | create a real matrix
90 90
91>>> mat 5 [1..15] 91>>> matrix 5 [1..15]
92(3><5) 92(3><5)
93 [ 1.0, 2.0, 3.0, 4.0, 5.0 93 [ 1.0, 2.0, 3.0, 4.0, 5.0
94 , 6.0, 7.0, 8.0, 9.0, 10.0 94 , 6.0, 7.0, 8.0, 9.0, 10.0
95 , 11.0, 12.0, 13.0, 14.0, 15.0 ] 95 , 11.0, 12.0, 13.0, 14.0, 15.0 ]
96 96
97-} 97-}
98mat 98matrix
99 :: Int -- ^ columns 99 :: Int -- ^ columns
100 -> [โ„] -- ^ elements 100 -> [โ„] -- ^ elements
101 -> Matrix โ„ 101 -> Matrix โ„
102mat c = reshape c . fromList 102matrix c = reshape c . fromList
103 103
104 104
105{- | print a real matrix with given number of digits after the decimal point 105{- | print a real matrix with given number of digits after the decimal point
diff --git a/packages/base/src/Numeric/Sparse.hs b/packages/base/src/Numeric/Sparse.hs
index f495e3a..f1516ec 100644
--- a/packages/base/src/Numeric/Sparse.hs
+++ b/packages/base/src/Numeric/Sparse.hs
@@ -168,7 +168,7 @@ gmXv Dense{..} v
168{- | general matrix - vector product 168{- | general matrix - vector product
169 169
170>>> let m = mkSparse [((0,999),1.0),((1,1999),2.0)] 170>>> let m = mkSparse [((0,999),1.0),((1,1999),2.0)]
171>>> m !#> vect[1..2000] 171>>> m !#> vector [1..2000]
172fromList [1000.0,4000.0] 172fromList [1000.0,4000.0]
173 173
174-} 174-}
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-11 09:59:51 +0000