summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2015-06-15 12:52:46 +0200
committerAlberto Ruiz <aruiz@um.es>2015-06-15 12:52:46 +0200
commit7376d022b12a27db5a396f89806a709555c1c522 (patch)
treef5bab9687cf775185f20dbc22713eddb360b495a
parent57487d828065ea219cdb33c9dc177b67c60b34c7 (diff)
documentation, more general cond, remove some unicode, minor changes
-rw-r--r--packages/base/CHANGELOG2
-rw-r--r--packages/base/src/Internal/CG.hs8
-rw-r--r--packages/base/src/Internal/Container.hs85
-rw-r--r--packages/base/src/Internal/Element.hs38
-rw-r--r--packages/base/src/Internal/Matrix.hs1
-rw-r--r--packages/base/src/Internal/Modular.hs9
-rw-r--r--packages/base/src/Internal/Numeric.hs12
-rw-r--r--packages/base/src/Internal/Util.hs38
-rw-r--r--packages/base/src/Numeric/LinearAlgebra.hs48
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Data.hs19
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/HMatrix.hs5
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Static.hs16
12 files changed, 217 insertions, 64 deletions
diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG
index 93b2594..c3e118c 100644
--- a/packages/base/CHANGELOG
+++ b/packages/base/CHANGELOG
@@ -5,7 +5,7 @@
5 5
6 * basic support of Int32 and Int64 elements 6 * basic support of Int32 and Int64 elements
7 7
8 * remap, ccompare, sortIndex 8 * remap, more general cond
9 9
10 * experimental support of type safe modular arithmetic, including linear 10 * experimental support of type safe modular arithmetic, including linear
11 systems and lu factorization 11 systems and lu factorization
diff --git a/packages/base/src/Internal/CG.hs b/packages/base/src/Internal/CG.hs
index fd14212..758d130 100644
--- a/packages/base/src/Internal/CG.hs
+++ b/packages/base/src/Internal/CG.hs
@@ -45,13 +45,13 @@ cg sym at a (CGState p r r2 x _) = CGState p' r' r'2 x' rdx
45 ap1 = a p 45 ap1 = a p
46 ap | sym = ap1 46 ap | sym = ap1
47 | otherwise = at ap1 47 | otherwise = at ap1
48 pap | sym = p <·> ap1 48 pap | sym = p <.> ap1
49 | otherwise = norm2 ap1 ** 2 49 | otherwise = norm2 ap1 ** 2
50 alpha = r2 / pap 50 alpha = r2 / pap
51 dx = scale alpha p 51 dx = scale alpha p
52 x' = x + dx 52 x' = x + dx
53 r' = r - scale alpha ap 53 r' = r - scale alpha ap
54 r'2 = r' <·> r' 54 r'2 = r' <.> r'
55 beta = r'2 / r2 55 beta = r'2 / r2
56 p' = r' + scale beta p 56 p' = r' + scale beta p
57 57
@@ -75,9 +75,9 @@ solveG mat ma meth rawb x0' ϵb ϵx
75 b = mat rawb 75 b = mat rawb
76 x0 = if x0' == 0 then konst 0 (dim b) else x0' 76 x0 = if x0' == 0 then konst 0 (dim b) else x0'
77 r0 = b - a x0 77 r0 = b - a x0
78 r20 = r0 <·> r0 78 r20 = r0 <.> r0
79 p0 = r0 79 p0 = r0
80 nb2 = b <·> b 80 nb2 = b <.> b
81 ok CGState {..} 81 ok CGState {..}
82 = cgr2 <nb2*ϵb**2 82 = cgr2 <nb2*ϵb**2
83 || cgdx < ϵx 83 || cgdx < ϵx
diff --git a/packages/base/src/Internal/Container.hs b/packages/base/src/Internal/Container.hs
index 7fe5758..1c158ff 100644
--- a/packages/base/src/Internal/Container.hs
+++ b/packages/base/src/Internal/Container.hs
@@ -53,28 +53,23 @@ linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n
53 53
54-------------------------------------------------------------------------------- 54--------------------------------------------------------------------------------
55 55
56infixl 7 <.> 56infixr 8 <.>
57-- | An infix synonym for 'dot' 57{- | An infix synonym for 'dot'
58(<.>) :: Numeric t => Vector t -> Vector t -> t
59(<.>) = dot
60 58
59>>> vector [1,2,3,4] <.> vector [-2,0,1,1]
605.0
61 61
62infixr 8 <·>, #> 62>>> let 𝑖 = 0:+1 :: C
63>>> fromList [1+𝑖,1] <.> fromList [1,1+𝑖]
642.0 :+ 0.0
63 65
64{- | infix synonym for 'dot' 66-}
65 67
66>>> vector [1,2,3,4] > vector [-2,0,1,1] 68(<.>) :: Numeric t => Vector t -> Vector t -> t
675.0 69(<.>) = dot
68 70
69>>> let 𝑖 = 0:+1 :: ℂ
70>>> fromList [1+𝑖,1] <·> fromList [1,1+𝑖]
712.0 :+ 0.0
72 71
73(the dot symbol "·" is obtained by Alt-Gr .)
74 72
75-}
76(<·>) :: Numeric t => Vector t -> Vector t -> t
77(<·>) = dot
78 73
79 74
80{- | infix synonym for 'app' 75{- | infix synonym for 'app'
@@ -91,6 +86,7 @@ infixr 8 <·>, #>
91fromList [140.0,320.0] 86fromList [140.0,320.0]
92 87
93-} 88-}
89infixr 8 #>
94(#>) :: Numeric t => Matrix t -> Vector t -> Vector t 90(#>) :: Numeric t => Matrix t -> Vector t -> Vector t
95(#>) = mXv 91(#>) = mXv
96 92
@@ -264,6 +260,24 @@ instance Numeric Z
264sortVector :: (Ord t, Element t) => Vector t -> Vector t 260sortVector :: (Ord t, Element t) => Vector t -> Vector t
265sortVector = sortV 261sortVector = sortV
266 262
263{- |
264
265>>> m <- randn 4 10
266>>> disp 2 m
2674x10
268-0.31 0.41 0.43 -0.19 -0.17 -0.23 -0.17 -1.04 -0.07 -1.24
269 0.26 0.19 0.14 0.83 -1.54 -0.09 0.37 -0.63 0.71 -0.50
270-0.11 -0.10 -1.29 -1.40 -1.04 -0.89 -0.68 0.35 -1.46 1.86
271 1.04 -0.29 0.19 -0.75 -2.20 -0.01 1.06 0.11 -2.09 -1.58
272
273>>> disp 2 $ m ?? (All, Pos $ sortIndex (m!1))
2744x10
275-0.17 -1.04 -1.24 -0.23 0.43 0.41 -0.31 -0.17 -0.07 -0.19
276-1.54 -0.63 -0.50 -0.09 0.14 0.19 0.26 0.37 0.71 0.83
277-1.04 0.35 1.86 -0.89 -1.29 -0.10 -0.11 -0.68 -1.46 -1.40
278-2.20 0.11 -1.58 -0.01 0.19 -0.29 1.04 1.06 -2.09 -0.75
279
280-}
267sortIndex :: (Ord t, Element t) => Vector t -> Vector I 281sortIndex :: (Ord t, Element t) => Vector t -> Vector I
268sortIndex = sortI 282sortIndex = sortI
269 283
@@ -273,11 +287,50 @@ ccompare = ccompare'
273cselect :: (Container c t) => c I -> c t -> c t -> c t -> c t 287cselect :: (Container c t) => c I -> c t -> c t -> c t -> c t
274cselect = cselect' 288cselect = cselect'
275 289
290{- | Extract elements from positions given in matrices of rows and columns.
291
292>>> r
293(3><3)
294 [ 1, 1, 1
295 , 1, 2, 2
296 , 1, 2, 3 ]
297>>> c
298(3><3)
299 [ 0, 1, 5
300 , 2, 2, 1
301 , 4, 4, 1 ]
302>>> m
303(4><6)
304 [ 0, 1, 2, 3, 4, 5
305 , 6, 7, 8, 9, 10, 11
306 , 12, 13, 14, 15, 16, 17
307 , 18, 19, 20, 21, 22, 23 ]
308
309>>> remap r c m
310(3><3)
311 [ 6, 7, 11
312 , 8, 14, 13
313 , 10, 16, 19 ]
314
315The indexes are autoconformable.
316
317>>> c'
318(3><1)
319 [ 1
320 , 2
321 , 4 ]
322>>> remap r c' m
323(3><3)
324 [ 7, 7, 7
325 , 8, 14, 14
326 , 10, 16, 22 ]
327
328-}
276remap :: Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t 329remap :: Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t
277remap i j m 330remap i j m
278 | minElement i >= 0 && maxElement i < fromIntegral (rows m) && 331 | minElement i >= 0 && maxElement i < fromIntegral (rows m) &&
279 minElement j >= 0 && maxElement j < fromIntegral (cols m) = remapM i' j' m 332 minElement j >= 0 && maxElement j < fromIntegral (cols m) = remapM i' j' m
280 | otherwise = error $ "out of range index in rmap" 333 | otherwise = error $ "out of range index in remap"
281 where 334 where
282 [i',j'] = conformMs [i,j] 335 [i',j'] = conformMs [i,j]
283 336
diff --git a/packages/base/src/Internal/Element.hs b/packages/base/src/Internal/Element.hs
index 4007491..51d5686 100644
--- a/packages/base/src/Internal/Element.hs
+++ b/packages/base/src/Internal/Element.hs
@@ -80,7 +80,7 @@ breakAt c l = (a++[c],tail b) where
80 (a,b) = break (==c) l 80 (a,b) = break (==c) l
81 81
82-------------------------------------------------------------------------------- 82--------------------------------------------------------------------------------
83 83-- | Specification of indexes for the operator '??'.
84data Extractor 84data Extractor
85 = All 85 = All
86 | Range Int Int Int 86 | Range Int Int Int
@@ -102,7 +102,32 @@ ppext (Drop n) = printf "Drop %d" n
102ppext (TakeLast n) = printf "TakeLast %d" n 102ppext (TakeLast n) = printf "TakeLast %d" n
103ppext (DropLast n) = printf "DropLast %d" n 103ppext (DropLast n) = printf "DropLast %d" n
104 104
105{- | General matrix slicing.
106
107>>> m
108(4><5)
109 [ 0, 1, 2, 3, 4
110 , 5, 6, 7, 8, 9
111 , 10, 11, 12, 13, 14
112 , 15, 16, 17, 18, 19 ]
113
114>>> m ?? (Take 3, DropLast 2)
115(3><3)
116 [ 0, 1, 2
117 , 5, 6, 7
118 , 10, 11, 12 ]
119
120>>> m ?? (Pos (idxs[2,1]), All)
121(2><5)
122 [ 10, 11, 12, 13, 14
123 , 5, 6, 7, 8, 9 ]
124
125>>> m ?? (PosCyc (idxs[-7,80]), Range 4 (-2) 0)
126(2><3)
127 [ 9, 7, 5
128 , 4, 2, 0 ]
105 129
130-}
106infixl 9 ?? 131infixl 9 ??
107(??) :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t 132(??) :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t
108 133
@@ -328,27 +353,30 @@ r >< c = f where
328 353
329---------------------------------------------------------------- 354----------------------------------------------------------------
330 355
331-- | Creates a matrix with the first n rows of another matrix
332takeRows :: Element t => Int -> Matrix t -> Matrix t 356takeRows :: Element t => Int -> Matrix t -> Matrix t
333takeRows n mt = subMatrix (0,0) (n, cols mt) mt 357takeRows n mt = subMatrix (0,0) (n, cols mt) mt
358
334-- | Creates a matrix with the last n rows of another matrix 359-- | Creates a matrix with the last n rows of another matrix
335takeLastRows :: Element t => Int -> Matrix t -> Matrix t 360takeLastRows :: Element t => Int -> Matrix t -> Matrix t
336takeLastRows n mt = subMatrix (rows mt - n, 0) (n, cols mt) mt 361takeLastRows n mt = subMatrix (rows mt - n, 0) (n, cols mt) mt
337-- | Creates a copy of a matrix without the first n rows 362
338dropRows :: Element t => Int -> Matrix t -> Matrix t 363dropRows :: Element t => Int -> Matrix t -> Matrix t
339dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt 364dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt
365
340-- | Creates a copy of a matrix without the last n rows 366-- | Creates a copy of a matrix without the last n rows
341dropLastRows :: Element t => Int -> Matrix t -> Matrix t 367dropLastRows :: Element t => Int -> Matrix t -> Matrix t
342dropLastRows n mt = subMatrix (0,0) (rows mt - n, cols mt) mt 368dropLastRows n mt = subMatrix (0,0) (rows mt - n, cols mt) mt
343-- |Creates a matrix with the first n columns of another matrix 369
344takeColumns :: Element t => Int -> Matrix t -> Matrix t 370takeColumns :: Element t => Int -> Matrix t -> Matrix t
345takeColumns n mt = subMatrix (0,0) (rows mt, n) mt 371takeColumns n mt = subMatrix (0,0) (rows mt, n) mt
372
346-- |Creates a matrix with the last n columns of another matrix 373-- |Creates a matrix with the last n columns of another matrix
347takeLastColumns :: Element t => Int -> Matrix t -> Matrix t 374takeLastColumns :: Element t => Int -> Matrix t -> Matrix t
348takeLastColumns n mt = subMatrix (0, cols mt - n) (rows mt, n) mt 375takeLastColumns n mt = subMatrix (0, cols mt - n) (rows mt, n) mt
349-- | Creates a copy of a matrix without the first n columns 376
350dropColumns :: Element t => Int -> Matrix t -> Matrix t 377dropColumns :: Element t => Int -> Matrix t -> Matrix t
351dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt 378dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt
379
352-- | Creates a copy of a matrix without the last n columns 380-- | Creates a copy of a matrix without the last n columns
353dropLastColumns :: Element t => Int -> Matrix t -> Matrix t 381dropLastColumns :: Element t => Int -> Matrix t -> Matrix t
354dropLastColumns n mt = subMatrix (0,0) (rows mt, cols mt - n) mt 382dropLastColumns n mt = subMatrix (0,0) (rows mt, cols mt - n) mt
diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs
index e4b1226..75e92a5 100644
--- a/packages/base/src/Internal/Matrix.hs
+++ b/packages/base/src/Internal/Matrix.hs
@@ -378,7 +378,6 @@ foreign import ccall unsafe "transL" ctransL :: TMM Z
378 378
379---------------------------------------------------------------------- 379----------------------------------------------------------------------
380 380
381-- | Extracts a submatrix from a matrix.
382subMatrix :: Element a 381subMatrix :: Element a
383 => (Int,Int) -- ^ (r0,c0) starting position 382 => (Int,Int) -- ^ (r0,c0) starting position
384 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix 383 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs
index 3b27310..0d906bb 100644
--- a/packages/base/src/Internal/Modular.hs
+++ b/packages/base/src/Internal/Modular.hs
@@ -9,8 +9,8 @@
9{-# LANGUAGE FlexibleInstances #-} 9{-# LANGUAGE FlexibleInstances #-}
10{-# LANGUAGE UndecidableInstances #-} 10{-# LANGUAGE UndecidableInstances #-}
11{-# LANGUAGE GADTs #-} 11{-# LANGUAGE GADTs #-}
12{-# LANGUAGE TypeFamilies #-} 12{-# LANGUAGE TypeFamilies #-}
13 13{-# LANGUAGE TypeOperators #-}
14 14
15{- | 15{- |
16Module : Internal.Modular 16Module : Internal.Modular
@@ -23,7 +23,7 @@ Proof of concept of statically checked modular arithmetic.
23-} 23-}
24 24
25module Internal.Modular( 25module Internal.Modular(
26 Mod 26 Mod, type (./.)
27) where 27) where
28 28
29import Internal.Vector 29import Internal.Vector
@@ -49,6 +49,9 @@ import Data.Complex
49newtype Mod (n :: Nat) t = Mod {unMod:: t} 49newtype Mod (n :: Nat) t = Mod {unMod:: t}
50 deriving (Storable) 50 deriving (Storable)
51 51
52infixr 5 ./.
53type (./.) x n = Mod n x
54
52instance (Integral t, Enum t, KnownNat m) => Enum (Mod m t) 55instance (Integral t, Enum t, KnownNat m) => Enum (Mod m t)
53 where 56 where
54 toEnum = l0 (\m x -> fromIntegral $ x `mod` (fromIntegral m)) 57 toEnum = l0 (\m x -> fromIntegral $ x `mod` (fromIntegral m))
diff --git a/packages/base/src/Internal/Numeric.hs b/packages/base/src/Internal/Numeric.hs
index ca17c23..efcde2c 100644
--- a/packages/base/src/Internal/Numeric.hs
+++ b/packages/base/src/Internal/Numeric.hs
@@ -481,14 +481,16 @@ step = step'
481-- , 0.0, 100.0, 7.0, 8.0 481-- , 0.0, 100.0, 7.0, 8.0
482-- , 0.0, 0.0, 100.0, 12.0 ] 482-- , 0.0, 0.0, 100.0, 12.0 ]
483-- 483--
484-- >>> let chop x = cond (abs x) 1E-6 0 0 x
485--
484cond 486cond
485 :: (Ord e, Container c e) 487 :: (Ord e, Container c e, Container c x)
486 => c e -- ^ a 488 => c e -- ^ a
487 -> c e -- ^ b 489 -> c e -- ^ b
488 -> c e -- ^ l 490 -> c x -- ^ l
489 -> c e -- ^ e 491 -> c x -- ^ e
490 -> c e -- ^ g 492 -> c x -- ^ g
491 -> c e -- ^ result 493 -> c x -- ^ result
492cond a b l e g = cselect' (ccompare' a b) l e g 494cond a b l e g = cselect' (ccompare' a b) l e g
493 495
494 496
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs
index 09ba21c..f08f710 100644
--- a/packages/base/src/Internal/Util.hs
+++ b/packages/base/src/Internal/Util.hs
@@ -85,7 +85,7 @@ type ℤ = Int
85type ℂ = Complex Double 85type ℂ = Complex Double
86 86
87-- | imaginary unit 87-- | imaginary unit
88iC :: 88iC :: C
89iC = 0:+1 89iC = 0:+1
90 90
91{- | Create a real vector. 91{- | Create a real vector.
@@ -94,7 +94,7 @@ iC = 0:+1
94fromList [1.0,2.0,3.0,4.0,5.0] 94fromList [1.0,2.0,3.0,4.0,5.0]
95 95
96-} 96-}
97vector :: [] -> Vector 97vector :: [R] -> Vector R
98vector = fromList 98vector = fromList
99 99
100{- | Create a real matrix. 100{- | Create a real matrix.
@@ -108,8 +108,8 @@ vector = fromList
108-} 108-}
109matrix 109matrix
110 :: Int -- ^ number of columns 110 :: Int -- ^ number of columns
111 -> [] -- ^ elements in row order 111 -> [R] -- ^ elements in row order
112 -> Matrix 112 -> Matrix R
113matrix c = reshape c . fromList 113matrix c = reshape c . fromList
114 114
115 115
@@ -260,34 +260,34 @@ norm = pnorm PNorm2
260 260
261class Normed a 261class Normed a
262 where 262 where
263 norm_0 :: a -> 263 norm_0 :: a -> R
264 norm_1 :: a -> 264 norm_1 :: a -> R
265 norm_2 :: a -> 265 norm_2 :: a -> R
266 norm_Inf :: a -> 266 norm_Inf :: a -> R
267 267
268 268
269instance Normed (Vector ) 269instance Normed (Vector R)
270 where 270 where
271 norm_0 v = sumElements (step (abs v - scalar (eps*normInf v))) 271 norm_0 v = sumElements (step (abs v - scalar (eps*normInf v)))
272 norm_1 = pnorm PNorm1 272 norm_1 = pnorm PNorm1
273 norm_2 = pnorm PNorm2 273 norm_2 = pnorm PNorm2
274 norm_Inf = pnorm Infinity 274 norm_Inf = pnorm Infinity
275 275
276instance Normed (Vector ) 276instance Normed (Vector C)
277 where 277 where
278 norm_0 v = sumElements (step (fst (fromComplex (abs v)) - scalar (eps*normInf v))) 278 norm_0 v = sumElements (step (fst (fromComplex (abs v)) - scalar (eps*normInf v)))
279 norm_1 = pnorm PNorm1 279 norm_1 = pnorm PNorm1
280 norm_2 = pnorm PNorm2 280 norm_2 = pnorm PNorm2
281 norm_Inf = pnorm Infinity 281 norm_Inf = pnorm Infinity
282 282
283instance Normed (Matrix ) 283instance Normed (Matrix R)
284 where 284 where
285 norm_0 = norm_0 . flatten 285 norm_0 = norm_0 . flatten
286 norm_1 = pnorm PNorm1 286 norm_1 = pnorm PNorm1
287 norm_2 = pnorm PNorm2 287 norm_2 = pnorm PNorm2
288 norm_Inf = pnorm Infinity 288 norm_Inf = pnorm Infinity
289 289
290instance Normed (Matrix ) 290instance Normed (Matrix C)
291 where 291 where
292 norm_0 = norm_0 . flatten 292 norm_0 = norm_0 . flatten
293 norm_1 = pnorm PNorm1 293 norm_1 = pnorm PNorm1
@@ -323,12 +323,22 @@ instance Normed (Vector (Complex Float))
323 norm_Inf = norm_Inf . double 323 norm_Inf = norm_Inf . double
324 324
325 325
326norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> 326norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R
327norm_Frob = norm_2 . flatten 327norm_Frob = norm_2 . flatten
328 328
329norm_nuclear :: Field t => Matrix t -> 329norm_nuclear :: Field t => Matrix t -> R
330norm_nuclear = sumElements . singularValues 330norm_nuclear = sumElements . singularValues
331 331
332{- | Check if the absolute value or complex magnitude is greater than a given threshold
333
334>>> magnit 1E-6 (1E-12 :: R)
335False
336>>> magnit 1E-6 (3+iC :: C)
337True
338>>> magnit 0 (3 :: I ./. 5)
339True
340
341-}
332magnit :: (Element t, Normed (Vector t)) => R -> t -> Bool 342magnit :: (Element t, Normed (Vector t)) => R -> t -> Bool
333magnit e x = norm_1 (fromList [x]) > e 343magnit e x = norm_1 (fromList [x]) > e
334 344
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index fe524cc..2e6b8ca 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -47,7 +47,7 @@ module Numeric.LinearAlgebra (
47 47
48 -- * Products 48 -- * Products
49 -- ** dot 49 -- ** dot
50 dot, (<·>), 50 dot, (<.>),
51 -- ** matrix-vector 51 -- ** matrix-vector
52 app, (#>), (<#), (!#>), 52 app, (#>), (<#), (!#>),
53 -- ** matrix-matrix 53 -- ** matrix-matrix
@@ -240,7 +240,53 @@ nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m)
240-- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'. 240-- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'.
241orth m = orthSVD (Left (1*eps)) m (leftSV m) 241orth m = orthSVD (Left (1*eps)) m (leftSV m)
242 242
243{- | Experimental implementation of LU factorization
244 working on any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.
245
246>>> let m = ident 5 + (5><5) [0..] :: Matrix (Z ./. 17)
247(5><5)
248 [ 1, 1, 2, 3, 4
249 , 5, 7, 7, 8, 9
250 , 10, 11, 13, 13, 14
251 , 15, 16, 0, 2, 2
252 , 3, 4, 5, 6, 8 ]
253
254>>> let (l,u,p,s) = luFact $ luPacked' m
255>>> l
256(5><5)
257 [ 1, 0, 0, 0, 0
258 , 6, 1, 0, 0, 0
259 , 12, 7, 1, 0, 0
260 , 7, 10, 7, 1, 0
261 , 8, 2, 6, 11, 1 ]
262>>> u
263(5><5)
264 [ 15, 16, 0, 2, 2
265 , 0, 13, 7, 13, 14
266 , 0, 0, 15, 0, 11
267 , 0, 0, 0, 15, 15
268 , 0, 0, 0, 0, 1 ]
269
270-}
243luPacked' x = mutable (luST (magnit 0)) x 271luPacked' x = mutable (luST (magnit 0)) x
244 272
273{- | Experimental implementation of gaussian elimination
274 working on any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'.
275
276>>> let a = (2><2) [1,2,3,5] :: Matrix (Z ./. 13)
277(2><2)
278 [ 1, 2
279 , 3, 5 ]
280>>> b
281(2><3)
282 [ 5, 1, 3
283 , 8, 6, 3 ]
284
285>>> let x = linearSolve' a b
286(2><3)
287 [ 4, 7, 4
288 , 7, 10, 6 ]
289
290-}
245linearSolve' x y = gaussElim x y 291linearSolve' x y = gaussElim x y
246 292
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index fffc2bd..2a45771 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -1,3 +1,5 @@
1{-# LANGUAGE TypeOperators #-}
2
1-------------------------------------------------------------------------------- 3--------------------------------------------------------------------------------
2{- | 4{- |
3Module : Numeric.LinearAlgebra.Data 5Module : Numeric.LinearAlgebra.Data
@@ -17,11 +19,11 @@ module Numeric.LinearAlgebra.Data(
17 -- | 1D arrays are storable vectors from the vector package. There is no distinction 19 -- | 1D arrays are storable vectors from the vector package. There is no distinction
18 -- between row and column vectors. 20 -- between row and column vectors.
19 21
20 fromList, toList, vector, (|>), 22 fromList, toList, (|>), vector, range, idxs,
21 23
22 -- * Matrix 24 -- * Matrix
23 25
24 matrix, (><), tr, tr', 26 (><), matrix, tr, tr',
25 27
26 -- * Dimensions 28 -- * Dimensions
27 29
@@ -43,7 +45,7 @@ module Numeric.LinearAlgebra.Data(
43 Indexable(..), 45 Indexable(..),
44 46
45 -- * Construction 47 -- * Construction
46 scalar, Konst(..), Build(..), assoc, accum, linspace, range, idxs, -- ones, zeros, 48 scalar, Konst(..), Build(..), assoc, accum, linspace, -- ones, zeros,
47 49
48 -- * Diagonal 50 -- * Diagonal
49 ident, diag, diagl, diagRect, takeDiag, 51 ident, diag, diagl, diagRect, takeDiag,
@@ -53,16 +55,19 @@ module Numeric.LinearAlgebra.Data(
53 55
54 -- * Matrix extraction 56 -- * Matrix extraction
55 Extractor(..), (??), 57 Extractor(..), (??),
58
56 takeRows, dropRows, takeColumns, dropColumns, 59 takeRows, dropRows, takeColumns, dropColumns,
57 subMatrix, (?), (¿), fliprl, flipud, remap, 60 subMatrix, (?), (¿), fliprl, flipud,
61
62 remap,
58 63
59 -- * Block matrix 64 -- * Block matrix
60 fromBlocks, (|||), (===), diagBlock, repmat, toBlocks, toBlocksEvery, 65 fromBlocks, (|||), (===), diagBlock, repmat, toBlocks, toBlocksEvery,
61 66
62 -- * Mapping functions 67 -- * Mapping functions
63 conj, cmap, cmod, 68 conj, cmap, cmod,
64 69
65 step, cond, ccompare, cselect, 70 step, cond,
66 71
67 -- * Find elements 72 -- * Find elements
68 find, maxIndex, minIndex, maxElement, minElement, 73 find, maxIndex, minIndex, maxElement, minElement,
@@ -87,7 +92,7 @@ module Numeric.LinearAlgebra.Data(
87 separable, 92 separable,
88 fromArray2D, 93 fromArray2D,
89 module Data.Complex, 94 module Data.Complex,
90 R,C,I,Z,Mod, 95 R,C,I,Z,Mod, type(./.),
91 Vector, Matrix, GMatrix, nRows, nCols 96 Vector, Matrix, GMatrix, nRows, nCols
92 97
93) where 98) where
diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
index 11c2487..8adaaaf 100644
--- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs
@@ -13,10 +13,13 @@ compatibility with previous version, to be removed
13 13
14module Numeric.LinearAlgebra.HMatrix ( 14module Numeric.LinearAlgebra.HMatrix (
15 module Numeric.LinearAlgebra, 15 module Numeric.LinearAlgebra,
16 (¦),(——),ℝ,ℂ, 16 (¦),(——),ℝ,ℂ,(<·>)
17) where 17) where
18 18
19import Numeric.LinearAlgebra 19import Numeric.LinearAlgebra
20import Internal.Util 20import Internal.Util
21 21
22infixr 8 <·>
23(<·>) :: Numeric t => Vector t -> Vector t -> t
24(<·>) = dot
22 25
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
index a657bd0..d0a790d 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs
@@ -44,7 +44,7 @@ module Numeric.LinearAlgebra.Static(
44 -- * Complex 44 -- * Complex
45 C, M, Her, her, 𝑖, 45 C, M, Her, her, 𝑖,
46 -- * Products 46 -- * Products
47 (<>),(#>),(<·>), 47 (<>),(#>),(<.>),
48 -- * Linear Systems 48 -- * Linear Systems
49 linSolve, (<\>), 49 linSolve, (<\>),
50 -- * Factorizations 50 -- * Factorizations
@@ -55,13 +55,13 @@ module Numeric.LinearAlgebra.Static(
55 Disp(..), Domain(..), 55 Disp(..), Domain(..),
56 withVector, withMatrix, 56 withVector, withMatrix,
57 toRows, toColumns, 57 toRows, toColumns,
58 Sized(..), Diag(..), Sym, sym, mTm, unSym 58 Sized(..), Diag(..), Sym, sym, mTm, unSym, (<·>)
59) where 59) where
60 60
61 61
62import GHC.TypeLits 62import GHC.TypeLits
63import Numeric.LinearAlgebra hiding ( 63import Numeric.LinearAlgebra hiding (
64 (<>),(#>),(<·>),Konst(..),diag, disp,(===),(|||), 64 (<>),(#>),(<.>),Konst(..),diag, disp,(===),(|||),
65 row,col,vector,matrix,linspace,toRows,toColumns, 65 row,col,vector,matrix,linspace,toRows,toColumns,
66 (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH', 66 (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',
67 eigenvalues,eigenvaluesSH,eigenvaluesSH',build, 67 eigenvalues,eigenvaluesSH,eigenvaluesSH',build,
@@ -207,6 +207,10 @@ infixr 8 <·>
207(<·>) :: R n -> R n -> ℝ 207(<·>) :: R n -> R n -> ℝ
208(<·>) = dotR 208(<·>) = dotR
209 209
210infixr 8 <.>
211(<.>) :: R n -> R n -> ℝ
212(<.>) = dotR
213
210-------------------------------------------------------------------------------- 214--------------------------------------------------------------------------------
211 215
212class Diag m d | m -> d 216class Diag m d | m -> d
@@ -496,7 +500,7 @@ appC m v = mkC (extract m LA.#> extract v)
496dotC :: KnownNat n => C n -> C n -> ℂ 500dotC :: KnownNat n => C n -> C n -> ℂ
497dotC (unwrap -> u) (unwrap -> v) 501dotC (unwrap -> u) (unwrap -> v)
498 | singleV u || singleV v = sumElements (conj u * v) 502 | singleV u || singleV v = sumElements (conj u * v)
499 | otherwise = u LA.<·> v 503 | otherwise = u LA.<.> v
500 504
501 505
502crossC :: C 3 -> C 3 -> C 3 506crossC :: C 3 -> C 3 -> C 3
@@ -584,12 +588,12 @@ test = (ok,info)
584 where 588 where
585 q = tm :: L 10 3 589 q = tm :: L 10 3
586 590
587 thingD = vjoin [ud1 u, 1] LA.<·> tr m LA.#> m LA.#> ud1 v 591 thingD = vjoin [ud1 u, 1] LA.<.> tr m LA.#> m LA.#> ud1 v
588 where 592 where
589 m = LA.matrix 3 [1..30] 593 m = LA.matrix 3 [1..30]
590 594
591 precS = (1::Double) + (2::Double) * ((1 :: R 3) * (u & 6)) <·> konst 2 #> v 595 precS = (1::Double) + (2::Double) * ((1 :: R 3) * (u & 6)) <·> konst 2 #> v
592 precD = 1 + 2 * vjoin[ud1 u, 6] LA.<·> LA.konst 2 (LA.size (ud1 u) +1, LA.size (ud1 v)) LA.#> ud1 v 596 precD = 1 + 2 * vjoin[ud1 u, 6] LA.<.> LA.konst 2 (LA.size (ud1 u) +1, LA.size (ud1 v)) LA.#> ud1 v
593 597
594 598
595splittest 599splittest