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