module Agda.Termination.SparseMatrix
(
Matrix
, matrixInvariant
, Size(..)
, sizeInvariant
, MIx (..)
, mIxInvariant
, fromLists
, fromIndexList
, toLists
, matrix
, matrixUsingRowGen
, size
, square
, isEmpty
, isSingleton
, add, intersectWith
, mul
, transpose
, diagonal
, addRow
, addColumn
, Agda.Termination.SparseMatrix.tests
) where
import Data.Array
import qualified Data.List as List
import Data.Maybe
import Data.Monoid
import Agda.Utils.Pretty hiding (isEmpty)
import Agda.Utils.QuickCheck
import Agda.Utils.TestHelpers
import Agda.Termination.Semiring (HasZero(..), SemiRing, Semiring)
import qualified Agda.Termination.Semiring as Semiring
#include "../undefined.h"
import Agda.Utils.Impossible
type TM = Matrix Integer Integer
data Size i = Size { rows :: i, cols :: i }
deriving (Eq, Ord, Show)
sizeInvariant :: (Ord i, Num i) => Size i -> Bool
sizeInvariant sz = rows sz >= 0 && cols sz >= 0
instance (Arbitrary i, Integral i) => Arbitrary (Size i) where
arbitrary = do
r <- natural
c <- natural
return $ Size { rows = fromInteger r, cols = fromInteger c }
instance CoArbitrary i => CoArbitrary (Size i) where
coarbitrary (Size rs cs) = coarbitrary rs . coarbitrary cs
prop_Arbitrary_Size :: Size Integer -> Bool
prop_Arbitrary_Size = sizeInvariant
toBounds :: Num i => Size i -> (MIx i, MIx i)
toBounds sz = (MIx { row = 1, col = 1 }, MIx { row = rows sz, col = cols sz })
data MIx i = MIx { row, col :: i }
deriving (Eq, Show, Ix, Ord)
instance (Arbitrary i, Integral i) => Arbitrary (MIx i) where
arbitrary = do
r <- positive
c <- positive
return $ MIx { row = r, col = c }
instance CoArbitrary i => CoArbitrary (MIx i) where
coarbitrary (MIx r c) = coarbitrary r . coarbitrary c
mIxInvariant :: (Ord i, Num i) => MIx i -> Bool
mIxInvariant i = row i >= 1 && col i >= 1
prop_Arbitrary_MIx :: MIx Integer -> Bool
prop_Arbitrary_MIx = mIxInvariant
data Matrix i b = M { size :: Size i, unM :: [(MIx i, b)] }
deriving (Eq, Ord, Functor)
matrixInvariant :: (Num i, Ix i) => Matrix i b -> Bool
matrixInvariant m = all (\ (MIx i j, b) -> 1 <= i && i <= rows sz
&& 1 <= j && j <= cols sz) (unM m)
&& strictlySorted (MIx 0 0) (unM m)
&& sizeInvariant sz
where sz = size m
strictlySorted :: (Ord i) => i -> [(i, b)] -> Bool
strictlySorted i [] = True
strictlySorted i ((i', b) : l) = i < i' && strictlySorted i' l
instance (Ord i, Integral i, Enum i, Ix i, Show i, Show b, HasZero b) => Show (Matrix i b) where
showsPrec _ m =
showString "Agda.Termination.Matrix.fromLists " . shows (size m) .
showString " " . shows (toLists m)
instance (Show i, Integral i, Ix i, HasZero b, Pretty b) =>
Pretty (Matrix i b) where
pretty = vcat . map (hsep . map pretty) . toLists
instance (Arbitrary i, Num i, Integral i, Arbitrary b, HasZero b)
=> Arbitrary (Matrix i b) where
arbitrary = matrix =<< arbitrary
instance (Show i, Ord i, Integral i, Enum i, Ix i, CoArbitrary b, HasZero b) => CoArbitrary (Matrix i b) where
coarbitrary m = coarbitrary (toLists m)
prop_Arbitrary_Matrix :: TM -> Bool
prop_Arbitrary_Matrix = matrixInvariant
matrixUsingRowGen :: (Arbitrary i, Integral i, Arbitrary b, HasZero b)
=> Size i
-> (i -> Gen [b])
-> Gen (Matrix i b)
matrixUsingRowGen sz rowGen = do
rows <- vectorOf (fromIntegral $ rows sz) (rowGen $ cols sz)
return $ fromLists sz rows
matrix :: (Arbitrary i, Integral i, Arbitrary b, HasZero b)
=> Size i -> Gen (Matrix i b)
matrix sz = matrixUsingRowGen sz (\n -> vectorOf (fromIntegral n) arbitrary)
prop_matrix sz = forAll (matrix sz :: Gen TM) $ \m ->
size m == sz
fromIndexList :: (Ord i, HasZero b) => Size i -> [(MIx i, b)] -> Matrix i b
fromIndexList sz = M sz . List.sortBy (\ (i,_) (j,_) -> compare i j) . filter (\ (i,b) -> b /= zeroElement)
prop_fromIndexList :: TM -> Bool
prop_fromIndexList m = matrixInvariant m' && m' == m
where vs = unM m
m' = fromIndexList (size m) vs
fromLists :: (Ord i, Num i, Enum i, HasZero b) => Size i -> [[b]] -> Matrix i b
fromLists sz bs = fromIndexList sz $ zip ([ MIx i j | i <- [1..rows sz]
, j <- [1..cols sz]]) (concat bs)
toSparseRows :: (Eq i, Num i, Enum i) => Matrix i b -> [(i,[(i,b)])]
toSparseRows m = aux 1 [] (unM m)
where aux i' [] [] = []
aux i' row [] = [(i', reverse row)]
aux i' row ((MIx i j, b) : m)
| i' == i = aux i' ((j,b):row) m
| otherwise = (i', reverse row) : aux i [(j,b)] m
blowUpSparseVec :: (Show i, Ord i, Num i, Enum i) => b -> i -> [(i,b)] -> [b]
blowUpSparseVec zero n l = aux 1 l
where aux i [] | i > n = []
| otherwise = zero : aux (i+1) []
aux i ((j,b):l) | i <= n && j == i = b : aux (succ i) l
aux i ((j,b):l) | i <= n && j >= i = zero : aux (succ i) ((j,b):l)
aux i l = __IMPOSSIBLE__
toLists :: (Show i, Ord i, Integral i, Enum i, Ix i, HasZero b) => Matrix i b -> [[b]]
toLists m =
blowUpSparseVec emptyRow nr $
map (\ (i,r) -> (i, blowUpSparseVec zeroElement nc r)) $ toSparseRows m
where
sz = size m
nr = rows sz
nc = cols sz
emptyRow = take (fromIntegral nc) $ repeat zeroElement
prop_fromLists_toLists :: TM -> Bool
prop_fromLists_toLists m = fromLists (size m) (toLists m) == m
prop_size :: TM -> Bool
prop_size m = sizeInvariant (size m)
prop_size_fromIndexList :: Size Int -> Bool
prop_size_fromIndexList sz =
size (fromIndexList sz ([] :: [(MIx Int, Integer)])) == sz
square :: Ix i => Matrix i b -> Bool
square m = rows (size m) == cols (size m)
isEmpty :: (Num i, Ix i) => Matrix i b -> Bool
isEmpty m = rows sz <= 0 || cols sz <= 0
where sz = size m
isSingleton :: (Num i, Ix i) => Matrix i b -> Maybe b
isSingleton m = if (rows sz == 1 || cols sz == 1) then
case unM m of
[(_,b)] -> Just b
_ -> __IMPOSSIBLE__
else Nothing
where sz = size m
transposeSize (Size { rows = n, cols = m }) = Size { rows = m, cols = n }
transpose m = M { size = transposeSize (size m)
, unM = List.sortBy (\ (i,a) (j,b) -> compare i j) $
map (\(MIx i j, b) -> (MIx j i, b)) $ unM m }
prop_transpose :: TM -> Bool
prop_transpose m = matrixInvariant m' && m == transpose m'
where m' = transpose m
add :: (Ord i) => (a -> a -> a) -> Matrix i a -> Matrix i a -> Matrix i a
add plus m1 m2 = M (supSize m1 m2) $ mergeAssocWith plus (unM m1) (unM m2)
supSize m1 m2 = Size { rows = r, cols = c }
where
sz1 = size m1
sz2 = size m2
r = max (rows sz1) (rows sz2)
c = max (cols sz1) (cols sz2)
mergeAssocWith :: (Ord i) => (a -> a -> a) -> [(i,a)] -> [(i,a)] -> [(i,a)]
mergeAssocWith f [] m = m
mergeAssocWith f l [] = l
mergeAssocWith f l@((i,a):l') m@((j,b):m')
| i < j = (i,a) : mergeAssocWith f l' m
| i > j = (j,b) : mergeAssocWith f l m'
| otherwise = (i, f a b) : mergeAssocWith f l' m'
intersectWith :: (Ord i) => (a -> a -> a) -> Matrix i a -> Matrix i a -> Matrix i a
intersectWith f m1 m2 = M (infSize m1 m2) $ interAssocWith f (unM m1) (unM m2)
infSize m1 m2 = Size { rows = r, cols = c }
where
sz1 = size m1
sz2 = size m2
r = min (rows sz1) (rows sz2)
c = min (cols sz1) (cols sz2)
interAssocWith :: (Ord i) => (a -> a -> a) -> [(i,a)] -> [(i,a)] -> [(i,a)]
interAssocWith f [] m = []
interAssocWith f l [] = []
interAssocWith f l@((i,a):l') m@((j,b):m')
| i < j = interAssocWith f l' m
| i > j = interAssocWith f l m'
| otherwise = (i, f a b) : interAssocWith f l' m'
prop_add sz =
forAll (three (matrix sz :: Gen TM)) $ \(m1, m2, m3) ->
let m' = add (+) m1 m2 in
associative (add (+)) m1 m2 m3 &&
commutative (add (+)) m1 m2 &&
matrixInvariant m' &&
size m' == size m1
mul :: (Enum i, Num i, Ix i, Eq a)
=> Semiring a -> Matrix i a -> Matrix i a -> Matrix i a
mul semiring m1 m2 = M (Size { rows = rows (size m1), cols = cols (size m2) }) $
filter (\ (i,b) -> b /= Semiring.zero semiring) $
[ (MIx i j, List.foldl' (Semiring.add semiring) (Semiring.zero semiring) $
map snd $ interAssocWith (Semiring.mul semiring) v w)
| (i,v) <- toSparseRows m1
, (j,w) <- toSparseRows $ transpose m2 ]
prop_mul sz =
sized $ \n -> resize (n `div` 2) $
forAll (two natural) $ \(c2, c3) ->
forAll (matrix sz :: Gen TM) $ \m1 ->
forAll (matrix (Size { rows = cols sz, cols = c2 })) $ \m2 ->
forAll (matrix (Size { rows = c2, cols = c3 })) $ \m3 ->
let m' = mult m1 m2 in
associative mult m1 m2 m3 &&
matrixInvariant m' &&
size m' == Size { rows = rows sz, cols = c2 }
where mult = mul Semiring.integerSemiring
diagonal :: (Show i, Enum i, Num i, Ix i, HasZero b) => Matrix i b -> Array i b
diagonal m =
listArray (1, n) $ blowUpSparseVec zeroElement n $
mapMaybe (\ ((MIx i j),b) -> if i==j then Just (i,b) else Nothing) $ unM m
where
sz = size m
r = rows sz
c = cols sz
n = max r c
prop_diagonal =
forAll natural $ \n ->
forAll (matrix (Size n n) :: Gen TM) $ \m ->
bounds (diagonal m) == (1, n)
addColumn :: (Num i, HasZero b) => b -> Matrix i b -> Matrix i b
addColumn x m | x == zeroElement = m { size = (size m) { cols = cols (size m) + 1 }}
| otherwise = __IMPOSSIBLE__
prop_addColumn :: TM -> Bool
prop_addColumn m =
matrixInvariant m'
&&
map init (toLists m') == toLists m
where
m' = addColumn zeroElement m
addRow :: (Num i, HasZero b) => b -> Matrix i b -> Matrix i b
addRow x m | x == zeroElement = m { size = (size m) { rows = rows (size m) + 1 }}
| otherwise = __IMPOSSIBLE__
prop_addRow :: TM -> Bool
prop_addRow m =
matrixInvariant m'
&&
init (toLists m') == toLists m
where
m' = addRow zeroElement m
tests :: IO Bool
tests = runTests "Agda.Termination.Matrix"
[ quickCheck' prop_transpose
, quickCheck' prop_Arbitrary_Size
, quickCheck' prop_Arbitrary_Matrix
, quickCheck' prop_Arbitrary_MIx
, quickCheck' prop_fromIndexList
, quickCheck' prop_matrix
, quickCheck' prop_size
, quickCheck' prop_size_fromIndexList
, quickCheck' prop_fromLists_toLists
, quickCheck' prop_add
, quickCheck' prop_mul
, quickCheck' prop_diagonal
, quickCheck' prop_addColumn
, quickCheck' prop_addRow
]