From 60dc8245539263899c083b2760310a86b14d367b Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Mon, 1 Jun 2015 18:51:17 +0200 Subject: generic gaussElim, modularTest --- packages/base/src/Data/Packed/Internal/Vector.hs | 2 +- packages/base/src/Numeric/LinearAlgebra.hs | 2 +- packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | 2 + packages/base/src/Numeric/LinearAlgebra/Util.hs | 74 ++++++++++++++++++++-- .../base/src/Numeric/LinearAlgebra/Util/Modular.hs | 26 ++++++-- packages/tests/src/Numeric/LinearAlgebra/Tests.hs | 5 ++ 6 files changed, 99 insertions(+), 12 deletions(-) diff --git a/packages/base/src/Data/Packed/Internal/Vector.hs b/packages/base/src/Data/Packed/Internal/Vector.hs index ac019a8..8cb77b0 100644 --- a/packages/base/src/Data/Packed/Internal/Vector.hs +++ b/packages/base/src/Data/Packed/Internal/Vector.hs @@ -106,7 +106,7 @@ inlinePerformIO :: IO a -> a inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r {-# INLINE inlinePerformIO #-} -{- | extracts the Vector elements to a list +{- extracts the Vector elements to a list >>> toList (linspace 5 (1,10)) [1.0,3.25,5.5,7.75,10.0] diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs index c7eb491..9dd7177 100644 --- a/packages/base/src/Numeric/LinearAlgebra.hs +++ b/packages/base/src/Numeric/LinearAlgebra.hs @@ -134,7 +134,7 @@ module Numeric.LinearAlgebra ( Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, -- * Misc - meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, + meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, gaussElim, ℝ,ℂ,iC, -- * Auxiliary classes Element, Container, Product, Numeric, LSDiv, 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 ] diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs index b138c2b..2be75de 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs @@ -383,6 +383,10 @@ intTest = utest "int ops" (fst $ checkT (undefined :: Matrix I)) -------------------------------------------------------------------------------- +modularTest = utest "modular ops" (fst $ checkT (undefined :: Matrix (F 13))) + +-------------------------------------------------------------------------------- + indexProp g f x = a1 == g a2 && a2 == a3 && b1 == g b2 && b2 == b3 where l = map g (toList (f x)) @@ -572,6 +576,7 @@ runTests n = do , sparseTest , staticTest , intTest + , modularTest ] when (errors c + failures c > 0) exitFailure return () -- cgit v1.2.3