diff options
-rw-r--r-- | packages/base/CHANGELOG | 6 | ||||
-rw-r--r-- | packages/base/src/Internal/Modular.hs | 6 | ||||
-rw-r--r-- | packages/base/src/Internal/Util.hs | 31 |
3 files changed, 24 insertions, 19 deletions
diff --git a/packages/base/CHANGELOG b/packages/base/CHANGELOG index b1d85fb..27c4d31 100644 --- a/packages/base/CHANGELOG +++ b/packages/base/CHANGELOG | |||
@@ -1,11 +1,11 @@ | |||
1 | 0.17.0.0 | 1 | 0.17.0.0 |
2 | -------- | 2 | -------- |
3 | 3 | ||
4 | * initial support of (C)Int elements | ||
5 | |||
6 | * improved matrix extraction using vectors of indexes (??) | 4 | * improved matrix extraction using vectors of indexes (??) |
5 | |||
6 | * basic support of Int32 and Int64 elements | ||
7 | 7 | ||
8 | * remap, ccompare | 8 | * remap, ccompare, sortIndex |
9 | 9 | ||
10 | * experimental support of type safe modular arithmetic | 10 | * experimental support of type safe modular arithmetic |
11 | 11 | ||
diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs index 6b34010..1289a21 100644 --- a/packages/base/src/Internal/Modular.hs +++ b/packages/base/src/Internal/Modular.hs | |||
@@ -137,7 +137,7 @@ instance forall m . KnownNat m => Container Vector (Mod m I) | |||
137 | scalar' x = fromList [x] | 137 | scalar' x = fromList [x] |
138 | konst' x = i2f . konst (unMod x) | 138 | konst' x = i2f . konst (unMod x) |
139 | build' n f = build n (fromIntegral . f) | 139 | build' n f = build n (fromIntegral . f) |
140 | cmap' = cmap | 140 | cmap' = mapVector |
141 | atIndex' x k = fromIntegral (atIndex (f2i x) k) | 141 | atIndex' x k = fromIntegral (atIndex (f2i x) k) |
142 | minIndex' = minIndex . f2i | 142 | minIndex' = minIndex . f2i |
143 | maxIndex' = maxIndex . f2i | 143 | maxIndex' = maxIndex . f2i |
@@ -177,7 +177,7 @@ instance forall m . KnownNat m => Container Vector (Mod m Z) | |||
177 | scalar' x = fromList [x] | 177 | scalar' x = fromList [x] |
178 | konst' x = i2f . konst (unMod x) | 178 | konst' x = i2f . konst (unMod x) |
179 | build' n f = build n (fromIntegral . f) | 179 | build' n f = build n (fromIntegral . f) |
180 | cmap' = cmap | 180 | cmap' = mapVector |
181 | atIndex' x k = fromIntegral (atIndex (f2i x) k) | 181 | atIndex' x k = fromIntegral (atIndex (f2i x) k) |
182 | minIndex' = minIndex . f2i | 182 | minIndex' = minIndex . f2i |
183 | maxIndex' = maxIndex . f2i | 183 | maxIndex' = maxIndex . f2i |
@@ -311,7 +311,7 @@ test = (ok, info) | |||
311 | 311 | ||
312 | print $ am <> gaussElim am bm - bm | 312 | print $ am <> gaussElim am bm - bm |
313 | print $ ad <> gaussElim ad bd - bd | 313 | print $ ad <> gaussElim ad bd - bd |
314 | 314 | ||
315 | print g | 315 | print g |
316 | print $ g <> g | 316 | print $ g <> g |
317 | print gm | 317 | print gm |
diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs index 1b639e1..b1fb800 100644 --- a/packages/base/src/Internal/Util.hs +++ b/packages/base/src/Internal/Util.hs | |||
@@ -512,7 +512,7 @@ block3x3 r nr c nc m = [[m ?? (er !! i, ec !! j) | j <- [0..2] ] | i <- [0..2] ] | |||
512 | er = [ Range 0 1 (r-1), Range r 1 (r+nr-1), Drop (nr+r) ] | 512 | er = [ Range 0 1 (r-1), Range r 1 (r+nr-1), Drop (nr+r) ] |
513 | ec = [ Range 0 1 (c-1), Range c 1 (c+nc-1), Drop (nc+c) ] | 513 | ec = [ Range 0 1 (c-1), Range c 1 (c+nc-1), Drop (nc+c) ] |
514 | 514 | ||
515 | view1 :: Numeric t => Matrix t -> Maybe (t, Vector t, Vector t, Matrix t) | 515 | view1 :: Numeric t => Matrix t -> Maybe (View1 t) |
516 | view1 m | 516 | view1 m |
517 | | rows m > 0 && cols m > 0 = Just (e, flatten m12, flatten m21 , m22) | 517 | | rows m > 0 && cols m > 0 = Just (e, flatten m12, flatten m21 , m22) |
518 | | otherwise = Nothing | 518 | | otherwise = Nothing |
@@ -520,21 +520,25 @@ view1 m | |||
520 | [[m11,m12],[m21,m22]] = block2x2 1 1 m | 520 | [[m11,m12],[m21,m22]] = block2x2 1 1 m |
521 | e = m11 `atIndex` (0, 0) | 521 | e = m11 `atIndex` (0, 0) |
522 | 522 | ||
523 | unView1 :: Numeric t => (t, Vector t, Vector t, Matrix t) -> Matrix t | 523 | unView1 :: Numeric t => View1 t -> Matrix t |
524 | unView1 (e,r,c,m) = fromBlocks [[scalar e, asRow r],[asColumn c, m]] | 524 | unView1 (e,r,c,m) = fromBlocks [[scalar e, asRow r],[asColumn c, m]] |
525 | 525 | ||
526 | type View1 t = (t, Vector t, Vector t, Matrix t) | ||
526 | 527 | ||
528 | foldMatrix :: Numeric t => (Matrix t -> Matrix t) -> (View1 t -> View1 t) -> (Matrix t -> Matrix t) | ||
527 | foldMatrix g f ( (f <$>) . view1 . g -> Just (e,r,c,m)) = unView1 (e, r, c, foldMatrix g f m) | 529 | foldMatrix g f ( (f <$>) . view1 . g -> Just (e,r,c,m)) = unView1 (e, r, c, foldMatrix g f m) |
528 | foldMatrix _ _ m = m | 530 | foldMatrix _ _ m = m |
529 | 531 | ||
530 | sortRowsBy h j m = m ?? (Pos (sortIndex (h (tr m ! j))), All) | ||
531 | |||
532 | splitColsAt n = (takeColumns n &&& dropColumns n) | ||
533 | 532 | ||
533 | swapMax k m | ||
534 | | rows m > 0 && j>0 = (j, m ?? (Pos (idxs swapped), All)) | ||
535 | | otherwise = (0,m) | ||
536 | where | ||
537 | j = maxIndex $ abs (tr m ! k) | ||
538 | swapped = j:[1..j-1] ++ 0:[j+1..rows m-1] | ||
534 | 539 | ||
535 | down a = foldMatrix g f a | 540 | down g a = foldMatrix g f a |
536 | where | 541 | where |
537 | g = sortRowsBy (negate.abs) 0 | ||
538 | f (e,r,c,m) | 542 | f (e,r,c,m) |
539 | | e /= 0 = (1, r', 0, m - outer c r') | 543 | | e /= 0 = (1, r', 0, m - outer c r') |
540 | | otherwise = error "singular!" | 544 | | otherwise = error "singular!" |
@@ -547,15 +551,16 @@ down a = foldMatrix g f a | |||
547 | -- @a <> gaussElim a b = b@ | 551 | -- @a <> gaussElim a b = b@ |
548 | -- | 552 | -- |
549 | gaussElim | 553 | gaussElim |
550 | :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t) | 554 | :: (Eq t, Fractional t, Num (Vector t), Numeric t) |
551 | => Matrix t -> Matrix t -> Matrix t | 555 | => Matrix t -> Matrix t -> Matrix t |
552 | 556 | ||
553 | gaussElim a b = r | 557 | gaussElim a b = flipudrl r |
554 | where | 558 | where |
555 | go x y = splitColsAt (cols a) (down $ fromBlocks [[x,y]]) | 559 | flipudrl = flipud . fliprl |
556 | (a1,b1) = go a b | 560 | splitColsAt n = (takeColumns n &&& dropColumns n) |
557 | ( _, r) = go (flipud . fliprl $ a1) (flipud . fliprl $ b1) | 561 | go f x y = splitColsAt (cols a) (down f $ fromBlocks [[x,y]]) |
558 | 562 | (a1,b1) = go (snd . swapMax 0) a b | |
563 | ( _, r) = go id (flipudrl $ a1) (flipudrl $ b1) | ||
559 | 564 | ||
560 | -------------------------------------------------------------------------------- | 565 | -------------------------------------------------------------------------------- |
561 | 566 | ||