generic gaussElim, modularTest

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
 Maintainer  :  Alberto Ruiz
 Stability   :  provisional
 
+compatibility with previous version, to be removed
+
 -}
 --------------------------------------------------------------------------------
 
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(
     -- ** 1D
     corr, conv, corrMin,
     -- ** 2D
-    corr2, conv2, separable
+    corr2, conv2, separable,
+    gaussElim
 ) where
 
 import Data.Packed.Numeric
@@ -65,7 +66,9 @@ import Numeric.LinearAlgebra.Util.Convolution
 import Control.Monad(when)
 import Text.Printf
 import Data.List.Split(splitOn)
-import Data.List(intercalate)
+import Data.List(intercalate,sortBy)
+import Data.Function(on)
+import Control.Arrow((&&&))
 
 type ℝ = Double
 type ℕ = Int
@@ -76,7 +79,7 @@ type ℂ = Complex Double
 iC :: ℂ
 iC = 0:+1
 
-{- | create a real vector
+{- | Create a real vector.
 
 >>> vector [1..5]
 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]
 vector :: [ℝ] -> Vector ℝ
 vector = fromList
 
-{- | create a real matrix
+{- | Create a real matrix.
 
 >>> matrix 5 [1..15]
 (3><5)
@@ -95,8 +98,8 @@ vector = fromList
 
 -}
 matrix
-  :: Int -- ^ columns
-  -> [ℝ] -- ^ elements
+  :: Int -- ^ number of columns
+  -> [ℝ] -- ^ elements in row order
   -> Matrix ℝ
 matrix c = reshape c . fromList
 
@@ -179,10 +182,22 @@ infixl 2 #
 a # b = fromBlocks [[a],[b]]
 
 -- | create a single row real matrix from a list
+--
+-- >>> row [2,3,1,8]
+-- (1><4)
+--  [ 2.0, 3.0, 1.0, 8.0 ]
+--
 row :: [Double] -> Matrix Double
 row = asRow . fromList
 
 -- | create a single column real matrix from a list
+--
+-- >>> col [7,-2,4]
+-- (3><1)
+--  [  7.0
+--  , -2.0
+--  ,  4.0 ]
+--
 col :: [Double] -> Matrix Double
 col = asColumn . fromList
 
@@ -477,6 +492,53 @@ dispShort maxr maxc dec m =
 
 --------------------------------------------------------------------------------
 
+-- | generic reference implementation of gaussian elimination
+--
+-- @a <> gauss a b = b@
+--
+gaussElim
+  :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
+  => Matrix t -> Matrix t -> Matrix t
+
+gaussElim x y = dropColumns (rows x) (flipud $ fromRows s2)
+  where
+    rs = toRows $ fromBlocks [[x , y]]
+    s1 = pivotDown (rows x) 0 rs
+    s2 = pivotUp (rows x-1) (reverse s1)
+
+pivotDown t n xs
+    | t == n    = []
+    | otherwise = y : pivotDown t (n+1) ys
+  where
+    y:ys = redu (pivot n xs)
+
+    pivot k = (const k &&& id)
+            . reverse . sortBy (compare `on` (abs. (!k)))  -- FIXME
+
+    redu (k,x:zs)
+        | p == 0 = error "gauss: singular!"  -- FIXME
+        | otherwise = u : map f zs
+      where
+        p = x!k
+        u = scale (recip (x!k)) x
+        f z = z - scale (z!k) u
+    redu (_,[]) = []
+
+
+pivotUp n xs
+    | n == -1 = []
+    | otherwise = y : pivotUp (n-1) ys
+  where
+    y:ys = redu' (n,xs)
+
+    redu' (k,x:zs) = u : map f zs
+      where
+        u = x
+        f z = z - scale (z!k) u
+    redu' (_,[]) = []
+
+--------------------------------------------------------------------------------
+
 instance Testable (Matrix I) where
    checkT _ = test
 
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(
 ) where
 
 import Data.Packed.Numeric
-import Numeric.LinearAlgebra.Util(Indexable(..),size)
+import Numeric.LinearAlgebra.Util(Indexable(..),gaussElim)
 import GHC.TypeLits
 import Data.Proxy(Proxy)
 import Foreign.ForeignPtr(castForeignPtr)
@@ -65,11 +65,14 @@ instance KnownNat m => Integral (F m)
       where
          (q,r) = quotRem (unMod a) (unMod b)
 
--- | this instance is only valid for prime m (not checked)
+-- | this instance is only valid for prime m
 instance KnownNat m => Fractional (F m)
   where
-    recip x = x^(m'-2)
+    recip x
+        | x*r == 1  = r
+        | otherwise = error $ show x ++" does not have a multiplicative inverse mod "++show m'
       where
+        r = x^(m'-2)
         m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
     fromRational x = fromInteger (numerator x) / fromInteger (denominator x)
 
@@ -124,7 +127,7 @@ instance Element (F n)
 instance forall m . KnownNat m => Container Vector (F m)
   where
     conj' = id
-    size' = size
+    size' = dim
     scale' s x = fromInt (scale (unMod s) (toInt x))
     addConstant c x = fromInt (addConstant (unMod c) (toInt x))
     add a b = fromInt (add (toInt a) (toInt b))
@@ -216,6 +219,18 @@ test = (ok, info)
   where
     v = fromList [3,-5,75] :: V 11
     m = (3><3) [1..]   :: M 11
+    
+    a = (3><3) [1,2 , 3
+               ,4,5 , 6
+               ,0,10,-3] :: Matrix I
+
+    b = (3><2) [0..] :: Matrix I
+    
+    am = fromInt a :: Matrix (F 13)
+    bm = fromInt b :: Matrix (F 13)
+    ad = fromInt a :: Matrix Double
+    bd = fromInt b :: Matrix Double
+    
     info = do
         print v
         print m
@@ -225,9 +240,12 @@ test = (ok, info)
         print $ m <> m
         print $ m #> v
 
+        print $ am <> gaussElim am bm - bm
+        print $ ad <> gaussElim ad bd - bd
 
     ok = and
       [ toInt (m #> v) == cmod 11 (toInt m #> toInt v )
+      , am <> gaussElim am bm == bm
       ]