summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--packages/base/src/Data/Packed/Internal/Vector.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/HMatrix.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Util.hs74
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs26
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests.hs5
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
106inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r 106inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
107{-# INLINE inlinePerformIO #-} 107{-# INLINE inlinePerformIO #-}
108 108
109{- | extracts the Vector elements to a list 109{- extracts the Vector elements to a list
110 110
111>>> toList (linspace 5 (1,10)) 111>>> toList (linspace 5 (1,10))
112[1.0,3.25,5.5,7.75,10.0] 112[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 (
134 Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, 134 Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
135 135
136 -- * Misc 136 -- * Misc
137 meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, 137 meanCov, rowOuters, pairwiseD2, unitary, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, gaussElim,
138 ℝ,ℂ,iC, 138 ℝ,ℂ,iC,
139 -- * Auxiliary classes 139 -- * Auxiliary classes
140 Element, Container, Product, Numeric, LSDiv, 140 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
6Maintainer : Alberto Ruiz 6Maintainer : Alberto Ruiz
7Stability : provisional 7Stability : provisional
8 8
9compatibility 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
58import Data.Packed.Numeric 59import Data.Packed.Numeric
@@ -65,7 +66,9 @@ import Numeric.LinearAlgebra.Util.Convolution
65import Control.Monad(when) 66import Control.Monad(when)
66import Text.Printf 67import Text.Printf
67import Data.List.Split(splitOn) 68import Data.List.Split(splitOn)
68import Data.List(intercalate) 69import Data.List(intercalate,sortBy)
70import Data.Function(on)
71import Control.Arrow((&&&))
69 72
70type ℝ = Double 73type ℝ = Double
71type ℕ = Int 74type ℕ = Int
@@ -76,7 +79,7 @@ type ℂ = Complex Double
76iC :: ℂ 79iC :: ℂ
77iC = 0:+1 80iC = 0:+1
78 81
79{- | create a real vector 82{- | Create a real vector.
80 83
81>>> vector [1..5] 84>>> vector [1..5]
82fromList [1.0,2.0,3.0,4.0,5.0] 85fromList [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]
85vector :: [ℝ] -> Vector ℝ 88vector :: [ℝ] -> Vector ℝ
86vector = fromList 89vector = 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-}
97matrix 100matrix
98 :: Int -- ^ columns 101 :: Int -- ^ number of columns
99 -> [ℝ] -- ^ elements 102 -> [ℝ] -- ^ elements in row order
100 -> Matrix ℝ 103 -> Matrix ℝ
101matrix c = reshape c . fromList 104matrix c = reshape c . fromList
102 105
@@ -179,10 +182,22 @@ infixl 2 #
179a # b = fromBlocks [[a],[b]] 182a # 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--
182row :: [Double] -> Matrix Double 190row :: [Double] -> Matrix Double
183row = asRow . fromList 191row = 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--
186col :: [Double] -> Matrix Double 201col :: [Double] -> Matrix Double
187col = asColumn . fromList 202col = 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--
499gaussElim
500 :: (Fractional t, Num (Vector t), Ord t, Indexable (Vector t) t, Numeric t)
501 => Matrix t -> Matrix t -> Matrix t
502
503gaussElim 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
509pivotDown 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
528pivotUp 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
480instance Testable (Matrix I) where 542instance 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
28import Data.Packed.Numeric 28import Data.Packed.Numeric
29import Numeric.LinearAlgebra.Util(Indexable(..),size) 29import Numeric.LinearAlgebra.Util(Indexable(..),gaussElim)
30import GHC.TypeLits 30import GHC.TypeLits
31import Data.Proxy(Proxy) 31import Data.Proxy(Proxy)
32import Foreign.ForeignPtr(castForeignPtr) 32import 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
69instance KnownNat m => Fractional (F m) 69instance 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)
124instance forall m . KnownNat m => Container Vector (F m) 127instance 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
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))
383 383
384-------------------------------------------------------------------------------- 384--------------------------------------------------------------------------------
385 385
386modularTest = utest "modular ops" (fst $ checkT (undefined :: Matrix (F 13)))
387
388--------------------------------------------------------------------------------
389
386indexProp g f x = a1 == g a2 && a2 == a3 && b1 == g b2 && b2 == b3 390indexProp g f x = a1 == g a2 && a2 == a3 && b1 == g b2 && b2 == b3
387 where 391 where
388 l = map g (toList (f x)) 392 l = map g (toList (f x))
@@ -572,6 +576,7 @@ runTests n = do
572 , sparseTest 576 , sparseTest
573 , staticTest 577 , staticTest
574 , intTest 578 , intTest
579 , modularTest
575 ] 580 ]
576 when (errors c + failures c > 0) exitFailure 581 when (errors c + failures c > 0) exitFailure
577 return () 582 return ()