diff options
Diffstat (limited to 'packages/base/src/Numeric/LinearAlgebra')
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | 2 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Util.hs | 74 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs | 26 |
3 files changed, 92 insertions, 10 deletions
diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs index 54f066b..a6383c1 100644 --- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs +++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | |||
@@ -6,6 +6,8 @@ License : BSD3 | |||
6 | Maintainer : Alberto Ruiz | 6 | Maintainer : Alberto Ruiz |
7 | Stability : provisional | 7 | Stability : provisional |
8 | 8 | ||
9 | compatibility with previous version, to be removed | ||
10 | |||
9 | -} | 11 | -} |
10 | -------------------------------------------------------------------------------- | 12 | -------------------------------------------------------------------------------- |
11 | 13 | ||
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs index 8d9f842..37cdc0f 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Util.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs | |||
@@ -52,7 +52,8 @@ module Numeric.LinearAlgebra.Util( | |||
52 | -- ** 1D | 52 | -- ** 1D |
53 | corr, conv, corrMin, | 53 | corr, conv, corrMin, |
54 | -- ** 2D | 54 | -- ** 2D |
55 | corr2, conv2, separable | 55 | corr2, conv2, separable, |
56 | gaussElim | ||
56 | ) where | 57 | ) where |
57 | 58 | ||
58 | import Data.Packed.Numeric | 59 | import Data.Packed.Numeric |
@@ -65,7 +66,9 @@ import Numeric.LinearAlgebra.Util.Convolution | |||
65 | import Control.Monad(when) | 66 | import Control.Monad(when) |
66 | import Text.Printf | 67 | import Text.Printf |
67 | import Data.List.Split(splitOn) | 68 | import Data.List.Split(splitOn) |
68 | import Data.List(intercalate) | 69 | import Data.List(intercalate,sortBy) |
70 | import Data.Function(on) | ||
71 | import Control.Arrow((&&&)) | ||
69 | 72 | ||
70 | type ℝ = Double | 73 | type ℝ = Double |
71 | type ℕ = Int | 74 | type ℕ = Int |
@@ -76,7 +79,7 @@ type ℂ = Complex Double | |||
76 | iC :: ℂ | 79 | iC :: ℂ |
77 | iC = 0:+1 | 80 | iC = 0:+1 |
78 | 81 | ||
79 | {- | create a real vector | 82 | {- | Create a real vector. |
80 | 83 | ||
81 | >>> vector [1..5] | 84 | >>> vector [1..5] |
82 | fromList [1.0,2.0,3.0,4.0,5.0] | 85 | fromList [1.0,2.0,3.0,4.0,5.0] |
@@ -85,7 +88,7 @@ fromList [1.0,2.0,3.0,4.0,5.0] | |||
85 | vector :: [ℝ] -> Vector ℝ | 88 | vector :: [ℝ] -> Vector ℝ |
86 | vector = fromList | 89 | vector = fromList |
87 | 90 | ||
88 | {- | create a real matrix | 91 | {- | Create a real matrix. |
89 | 92 | ||
90 | >>> matrix 5 [1..15] | 93 | >>> matrix 5 [1..15] |
91 | (3><5) | 94 | (3><5) |
@@ -95,8 +98,8 @@ vector = fromList | |||
95 | 98 | ||
96 | -} | 99 | -} |
97 | matrix | 100 | matrix |
98 | :: Int -- ^ columns | 101 | :: Int -- ^ number of columns |
99 | -> [ℝ] -- ^ elements | 102 | -> [ℝ] -- ^ elements in row order |
100 | -> Matrix ℝ | 103 | -> Matrix ℝ |
101 | matrix c = reshape c . fromList | 104 | matrix c = reshape c . fromList |
102 | 105 | ||
@@ -179,10 +182,22 @@ infixl 2 # | |||
179 | a # b = fromBlocks [[a],[b]] | 182 | a # b = fromBlocks [[a],[b]] |
180 | 183 | ||
181 | -- | create a single row real matrix from a list | 184 | -- | create a single row real matrix from a list |
185 | -- | ||
186 | -- >>> row [2,3,1,8] | ||
187 | -- (1><4) | ||
188 | -- [ 2.0, 3.0, 1.0, 8.0 ] | ||
189 | -- | ||
182 | row :: [Double] -> Matrix Double | 190 | row :: [Double] -> Matrix Double |
183 | row = asRow . fromList | 191 | row = asRow . fromList |
184 | 192 | ||
185 | -- | create a single column real matrix from a list | 193 | -- | create a single column real matrix from a list |
194 | -- | ||
195 | -- >>> col [7,-2,4] | ||
196 | -- (3><1) | ||
197 | -- [ 7.0 | ||
198 | -- , -2.0 | ||
199 | -- , 4.0 ] | ||
200 | -- | ||
186 | col :: [Double] -> Matrix Double | 201 | col :: [Double] -> Matrix Double |
187 | col = asColumn . fromList | 202 | col = asColumn . fromList |
188 | 203 | ||
@@ -477,6 +492,53 @@ dispShort maxr maxc dec m = | |||
477 | 492 | ||
478 | -------------------------------------------------------------------------------- | 493 | -------------------------------------------------------------------------------- |
479 | 494 | ||
495 | -- | generic reference implementation of gaussian elimination | ||
496 | -- | ||
497 | -- @a <> gauss a b = b@ | ||
498 | -- | ||
499 | gaussElim | ||
500 | :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t) | ||
501 | => Matrix t -> Matrix t -> Matrix t | ||
502 | |||
503 | gaussElim x y = dropColumns (rows x) (flipud $ fromRows s2) | ||
504 | where | ||
505 | rs = toRows $ fromBlocks [[x , y]] | ||
506 | s1 = pivotDown (rows x) 0 rs | ||
507 | s2 = pivotUp (rows x-1) (reverse s1) | ||
508 | |||
509 | pivotDown t n xs | ||
510 | | t == n = [] | ||
511 | | otherwise = y : pivotDown t (n+1) ys | ||
512 | where | ||
513 | y:ys = redu (pivot n xs) | ||
514 | |||
515 | pivot k = (const k &&& id) | ||
516 | . reverse . sortBy (compare `on` (abs. (!k))) -- FIXME | ||
517 | |||
518 | redu (k,x:zs) | ||
519 | | p == 0 = error "gauss: singular!" -- FIXME | ||
520 | | otherwise = u : map f zs | ||
521 | where | ||
522 | p = x!k | ||
523 | u = scale (recip (x!k)) x | ||
524 | f z = z - scale (z!k) u | ||
525 | redu (_,[]) = [] | ||
526 | |||
527 | |||
528 | pivotUp n xs | ||
529 | | n == -1 = [] | ||
530 | | otherwise = y : pivotUp (n-1) ys | ||
531 | where | ||
532 | y:ys = redu' (n,xs) | ||
533 | |||
534 | redu' (k,x:zs) = u : map f zs | ||
535 | where | ||
536 | u = x | ||
537 | f z = z - scale (z!k) u | ||
538 | redu' (_,[]) = [] | ||
539 | |||
540 | -------------------------------------------------------------------------------- | ||
541 | |||
480 | instance Testable (Matrix I) where | 542 | instance Testable (Matrix I) where |
481 | checkT _ = test | 543 | checkT _ = test |
482 | 544 | ||
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs index 18b635d..ea4a668 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs | |||
@@ -26,7 +26,7 @@ module Numeric.LinearAlgebra.Util.Modular( | |||
26 | ) where | 26 | ) where |
27 | 27 | ||
28 | import Data.Packed.Numeric | 28 | import Data.Packed.Numeric |
29 | import Numeric.LinearAlgebra.Util(Indexable(..),size) | 29 | import Numeric.LinearAlgebra.Util(Indexable(..),gaussElim) |
30 | import GHC.TypeLits | 30 | import GHC.TypeLits |
31 | import Data.Proxy(Proxy) | 31 | import Data.Proxy(Proxy) |
32 | import Foreign.ForeignPtr(castForeignPtr) | 32 | import Foreign.ForeignPtr(castForeignPtr) |
@@ -65,11 +65,14 @@ instance KnownNat m => Integral (F m) | |||
65 | where | 65 | where |
66 | (q,r) = quotRem (unMod a) (unMod b) | 66 | (q,r) = quotRem (unMod a) (unMod b) |
67 | 67 | ||
68 | -- | this instance is only valid for prime m (not checked) | 68 | -- | this instance is only valid for prime m |
69 | instance KnownNat m => Fractional (F m) | 69 | instance KnownNat m => Fractional (F m) |
70 | where | 70 | where |
71 | recip x = x^(m'-2) | 71 | recip x |
72 | | x*r == 1 = r | ||
73 | | otherwise = error $ show x ++" does not have a multiplicative inverse mod "++show m' | ||
72 | where | 74 | where |
75 | r = x^(m'-2) | ||
73 | m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int | 76 | m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int |
74 | fromRational x = fromInteger (numerator x) / fromInteger (denominator x) | 77 | fromRational x = fromInteger (numerator x) / fromInteger (denominator x) |
75 | 78 | ||
@@ -124,7 +127,7 @@ instance Element (F n) | |||
124 | instance forall m . KnownNat m => Container Vector (F m) | 127 | instance forall m . KnownNat m => Container Vector (F m) |
125 | where | 128 | where |
126 | conj' = id | 129 | conj' = id |
127 | size' = size | 130 | size' = dim |
128 | scale' s x = fromInt (scale (unMod s) (toInt x)) | 131 | scale' s x = fromInt (scale (unMod s) (toInt x)) |
129 | addConstant c x = fromInt (addConstant (unMod c) (toInt x)) | 132 | addConstant c x = fromInt (addConstant (unMod c) (toInt x)) |
130 | add a b = fromInt (add (toInt a) (toInt b)) | 133 | add a b = fromInt (add (toInt a) (toInt b)) |
@@ -216,6 +219,18 @@ test = (ok, info) | |||
216 | where | 219 | where |
217 | v = fromList [3,-5,75] :: V 11 | 220 | v = fromList [3,-5,75] :: V 11 |
218 | m = (3><3) [1..] :: M 11 | 221 | m = (3><3) [1..] :: M 11 |
222 | |||
223 | a = (3><3) [1,2 , 3 | ||
224 | ,4,5 , 6 | ||
225 | ,0,10,-3] :: Matrix I | ||
226 | |||
227 | b = (3><2) [0..] :: Matrix I | ||
228 | |||
229 | am = fromInt a :: Matrix (F 13) | ||
230 | bm = fromInt b :: Matrix (F 13) | ||
231 | ad = fromInt a :: Matrix Double | ||
232 | bd = fromInt b :: Matrix Double | ||
233 | |||
219 | info = do | 234 | info = do |
220 | print v | 235 | print v |
221 | print m | 236 | print m |
@@ -225,9 +240,12 @@ test = (ok, info) | |||
225 | print $ m <> m | 240 | print $ m <> m |
226 | print $ m #> v | 241 | print $ m #> v |
227 | 242 | ||
243 | print $ am <> gaussElim am bm - bm | ||
244 | print $ ad <> gaussElim ad bd - bd | ||
228 | 245 | ||
229 | ok = and | 246 | ok = and |
230 | [ toInt (m #> v) == cmod 11 (toInt m #> toInt v ) | 247 | [ toInt (m #> v) == cmod 11 (toInt m #> toInt v ) |
248 | , am <> gaussElim am bm == bm | ||
231 | ] | 249 | ] |
232 | 250 | ||
233 | 251 | ||