module Agda.Termination.SparseMatrix
(
Matrix
, matrixInvariant
, Size(..)
, sizeInvariant
, MIx (..)
, mIxInvariant
, fromLists
, fromIndexList
, toLists
, matrix
, matrixUsingRowGen
, size
, square
, isEmpty
, isSingleton
, zipMatrices
, add, intersectWith
, mul
, transpose
, Diagonal(..)
, addRow
, addColumn
, Agda.Termination.SparseMatrix.tests
) where
import Data.Array
import Data.Function
import qualified Data.List as List
import Data.Maybe
import Data.Monoid
import Data.Foldable (Foldable)
import qualified Data.Foldable as Fold
import Data.Traversable (Traversable)
import qualified Text.PrettyPrint.Boxes as Boxes
import Agda.Termination.Semiring (HasZero(..), Semiring)
import qualified Agda.Termination.Semiring as Semiring
import Agda.Utils.Functor
import Agda.Utils.List
import Agda.Utils.Maybe
import Agda.Utils.Monad
import Agda.Utils.PartialOrd
import Agda.Utils.Pretty hiding (isEmpty)
import Agda.Utils.QuickCheck
import Agda.Utils.TestHelpers
import Agda.Utils.Tuple
#include "undefined.h"
import Agda.Utils.Impossible
data Size i = Size
{ rows :: i
, cols :: i
}
deriving (Eq, Ord, Show)
data MIx i = MIx
{ row :: i
, col :: i
}
deriving (Eq, Ord, Show, Ix)
toBounds :: Num i => Size i -> (MIx i, MIx i)
toBounds sz = (MIx { row = 1, col = 1 }, MIx { row = rows sz, col = cols sz })
data Matrix i b = Matrix
{ size :: Size i
, unM :: [(MIx i, b)]
}
deriving (Eq, Ord, Functor, Foldable, Traversable)
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
supSize :: Ord i => Matrix i a -> Matrix i b -> Size i
supSize (Matrix (Size r1 c1) _) (Matrix (Size r2 c2) _) =
Size (max r1 r2) (max c1 c2)
infSize :: Ord i => Matrix i a -> Matrix i b -> Size i
infSize (Matrix (Size r1 c1) _) (Matrix (Size r2 c2) _) =
Size (min r1 r2) (min c1 c2)
fromIndexList :: (Ord i, HasZero b) => Size i -> [(MIx i, b)] -> Matrix i b
fromIndexList sz = Matrix sz
. List.sortBy (compare `on` fst)
. filter ((zeroElement /=) . snd)
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) => Matrix i b -> [(i,[(i,b)])]
toSparseRows (Matrix _ []) = []
toSparseRows (Matrix _ ((MIx i j, b) : m)) = aux i [(j,b)] 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
toSparseRows' :: (Eq i) => Matrix i b -> [(i,[(i,b)])]
toSparseRows' (Matrix _ m) =
for (List.groupBy ((==) `on` (row . fst)) m) $ \ ((MIx i j, b) : vs) ->
(i, (j,b) : map (mapFst col) vs)
prop_toSparseRows :: TM -> Bool
prop_toSparseRows m = toSparseRows m == toSparseRows' m
blowUpSparseVec :: (Integral i) => b -> i -> [(i,b)] -> [b]
blowUpSparseVec zero n l = aux 1 l
where aux i [] = List.genericReplicate (n + 1 i) zero
aux i l@((j,b):l')
| i > j || i > n = __IMPOSSIBLE__
| i == j = b : aux (i + 1) l'
| otherwise = zero : aux (i + 1) l
blowUpSparseVec' :: (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 :: (Integral i, HasZero b) => Matrix i b -> [[b]]
toLists m@(Matrix size@(Size nrows ncols) _) =
blowUpSparseVec emptyRow nrows $
map (mapSnd (blowUpSparseVec zeroElement ncols)) $ toSparseRows m
where
emptyRow = List.genericReplicate ncols zeroElement
isSingleton :: (Eq i, Num i, HasZero b) => Matrix i b -> Maybe b
isSingleton (Matrix (Size 1 1) [(_,b)]) = Just b
isSingleton (Matrix (Size 1 1) [] ) = Just zeroElement
isSingleton (Matrix (Size 1 1) _ ) = __IMPOSSIBLE__
isSingleton _ = Nothing
class Diagonal m e | m -> e where
diagonal :: m -> [e]
instance (Integral i, HasZero b) => Diagonal (Matrix i b) b where
diagonal (Matrix (Size r c) m) =
blowUpSparseVec zeroElement (min r c) $
mapMaybe (\ ((MIx i j), b) -> if i==j then Just (i,b) else Nothing) m
class Transpose a where
transpose :: a -> a
instance Transpose (Size i) where
transpose (Size n m) = Size m n
instance Transpose (MIx i) where
transpose (MIx i j) = MIx j i
instance Ord i => Transpose (Matrix i b) where
transpose (Matrix size m) =
Matrix (transpose size) $
List.sortBy (compare `on` fst) $
map (mapFst transpose) m
zipAssocWith :: (Ord i)
=> ([(i,a)] -> [(i,c)])
-> ([(i,b)] -> [(i,c)])
-> (a -> Maybe c)
-> (b -> Maybe c)
-> (a -> b -> Maybe c)
-> [(i,a)] -> [(i,b)] -> [(i,c)]
zipAssocWith fs gs f g h = merge
where
merge m1 [] = mapMaybe (\ (i, a) -> (i,) <$> f a) m1
merge [] m2 = mapMaybe (\ (i, b) -> (i,) <$> g b) m2
merge m1@((i,a):m1') m2@((j,b):m2') =
case compare i j of
LT -> mcons ((i,) <$> f a) $ merge m1' m2
GT -> mcons ((j,) <$> g b) $ merge m1 m2'
EQ -> mcons ((i,) <$> h a b) $ merge m1' m2'
unionAssocWith :: (Ord i)
=> (a -> Maybe c)
-> (b -> Maybe c)
-> (a -> b -> Maybe c)
-> [(i,a)] -> [(i,b)] -> [(i,c)]
unionAssocWith f g h = zipAssocWith (map_ f) (map_ g) f g h
where
map_ f = mapMaybe (\ (i, a) -> (i,) <$> f a)
zipMatrices :: forall a b c i . (Ord i)
=> (a -> c)
-> (b -> c)
-> (a -> b -> c)
-> (c -> Bool)
-> Matrix i a -> Matrix i b -> Matrix i c
zipMatrices f g h zero m1 m2 = Matrix (supSize m1 m2) $
unionAssocWith (drop0 . f) (drop0 . g) (\ a -> drop0 . h a) (unM m1) (unM m2)
where
drop0 = filterMaybe (not . zero)
add :: (Ord i, HasZero a) => (a -> a -> a) -> Matrix i a -> Matrix i a -> Matrix i a
add plus = zipMatrices id id plus (== zeroElement)
intersectWith :: (Ord i) => (a -> a -> a) -> Matrix i a -> Matrix i a -> Matrix i a
intersectWith f m1 m2 = Matrix (infSize m1 m2) $ interAssocWith f (unM m1) (unM m2)
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'
interAssocWith2 :: Ord i => (a -> a -> a) -> [(i,a)] -> [(i,a)] -> [(i,a)]
interAssocWith2 f = zipAssocWith (const []) (const []) (const Nothing) (const Nothing) (\ a -> Just . f a)
prop_interAssocWith_correct2 :: [(Int,Int)] -> [(Int,Int)] -> Bool
prop_interAssocWith_correct2 xs ys =
interAssocWith (*) xs ys == interAssocWith2 (*) xs ys
where
l = List.sortBy (compare `on` fst) xs
l' = List.sortBy (compare `on` fst) ys
mul :: (Ix i, Eq a)
=> Semiring a -> Matrix i a -> Matrix i a -> Matrix i a
mul semiring m1 m2 = Matrix (Size { rows = rows (size m1), cols = cols (size m2) }) $
[ (MIx i j, b)
| (i,v) <- toSparseRows m1
, (j,w) <- toSparseRows $ transpose m2
, let b = inner v w
, b /= zero
]
where
zero = Semiring.zero semiring
plus = Semiring.add semiring
times = Semiring.mul semiring
inner v w = List.foldl' plus zero $
map snd $ interAssocWith times v w
instance (Ord i, PartialOrd a) => PartialOrd (Matrix i a) where
comparable m n
| size m /= size n = POAny
| otherwise = Fold.fold $
zipMatrices onlym onlyn both trivial m n
where
onlym o = POGT
onlyn o = POLT
both = comparable
trivial = (==mempty)
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
instance (Integral i, HasZero b, Show i, Show b) => Show (Matrix i b) where
showsPrec _ m =
showString "Agda.Termination.SparseMatrix.fromLists " . shows (size m) .
showString " " . shows (toLists m)
instance (Integral i, HasZero b, Pretty b) =>
Pretty (Matrix i b) where
pretty = vcat
. map text
. lines
. Boxes.render
. Boxes.hsep 1 Boxes.right
. map ( Boxes.vcat Boxes.right
. map ( Boxes.alignHoriz Boxes.right 2
. Boxes.text . render . pretty
)
)
. toLists
. transpose
instance 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
sizeInvariant :: (Ord i, Num i) => Size i -> Bool
sizeInvariant sz = rows sz >= 0 && cols sz >= 0
prop_Arbitrary_Size :: Size Integer -> Bool
prop_Arbitrary_Size = sizeInvariant
prop_size :: TM -> Bool
prop_size m = sizeInvariant (size m)
instance Integral i => Arbitrary (MIx i) where
arbitrary = MIx <$> positive <*> positive
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
matrixInvariant :: (Num i, Ix i, HasZero b) => Matrix i b -> Bool
matrixInvariant (Matrix size@Size{ rows, cols} m) =
sizeInvariant size
&& all inBounds m
&& all nonZero m
&& strictlySorted (MIx 0 0) m
where
inBounds (MIx i j, _) = 1 <= i && i <= rows
&& 1 <= j && j <= cols
nonZero (_, b) = b /= zeroElement
strictlySorted :: (Ord i) => i -> [(i, b)] -> Bool
strictlySorted i [] = True
strictlySorted i ((i', _) : l) = i < i' && strictlySorted i' l
matrixUsingRowGen :: (Integral i, 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 :: (Integral i, Arbitrary b, HasZero b)
=> Size i -> Gen (Matrix i b)
matrix sz = matrixUsingRowGen sz (\n -> vectorOf (fromIntegral n) arbitrary)
prop_matrix :: Size Int -> Property
prop_matrix sz = forAll (matrix sz :: Gen TM) $ \ m -> size m == sz
instance (Integral i, Arbitrary b, HasZero b)
=> Arbitrary (Matrix i b) where
arbitrary = matrix =<< arbitrary
instance (Integral i, CoArbitrary b, HasZero b) => CoArbitrary (Matrix i b) where
coarbitrary m = coarbitrary (toLists m)
type TM = Matrix Int Int
prop_Arbitrary_Matrix :: TM -> Bool
prop_Arbitrary_Matrix = matrixInvariant
prop_fromIndexList :: TM -> Bool
prop_fromIndexList m@(Matrix size vs) = fromIndexList size vs == m
prop_fromLists_toLists :: TM -> Bool
prop_fromLists_toLists m = fromLists (size m) (toLists m) == m
prop_isSingleton :: Int -> Bool
prop_isSingleton b = Just b == (isSingleton (fromLists (Size 1 1) [[b]] :: TM))
prop_diagonal :: TM -> Bool
prop_diagonal m@(Matrix (Size r c) _) =
length (diagonal m) == min r c
prop_transpose :: TM -> Bool
prop_transpose m = matrixInvariant m' && m == transpose m'
where m' = transpose m
zipMatrices' :: forall a b c i . (Ord i)
=> (a -> c)
-> (b -> c)
-> (a -> b -> c)
-> (c -> Bool)
-> Matrix i a -> Matrix i b -> Matrix i c
zipMatrices' f g h zero m1 m2 = Matrix (supSize m1 m2) (merge (unM m1) (unM m2))
where
merge :: [(MIx i,a)] -> [(MIx i,b)] -> [(MIx i,c)]
merge [] m2 = filter (not . zero . snd) $ map (mapSnd g) m2
merge m1 [] = filter (not . zero . snd) $ map (mapSnd f) m1
merge m1@((i,a):m1') m2@((j,b):m2') =
case compare i j of
LT -> if zero c then r else (i,c) : r where c = f a ; r = merge m1' m2
GT -> if zero c then r else (j,c) : r where c = g b ; r = merge m1 m2'
EQ -> if zero c then r else (i,c) : r where c = h a b ; r = merge m1' m2'
prop_zipMatrices_correct :: TM -> TM -> Bool
prop_zipMatrices_correct m1 m2 =
zipMatrices id id (+) (== 0) m1 m2 == zipMatrices' id id (+) (== 0) m1 m2
prop_add :: Size Int -> Property
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
prop_add_correct :: TM -> TM -> Bool
prop_add_correct m1 m2 = toLists (add (+) m1 m2) == toLists (add' (+) m1 m2)
where
add' :: (Ord i) => (a -> a -> a) -> Matrix i a -> Matrix i a -> Matrix i a
add' plus m1 m2 = Matrix (supSize m1 m2) $ mergeAssocWith plus (unM m1) (unM m2)
where
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'
interAssocWith' :: (Eq i) => (a -> a -> a) -> [(i,a)] -> [(i,a)] -> [(i,a)]
interAssocWith' f l l' = [ (i, f a b) | (i,a) <- l, (j,b) <- l', i == j ]
prop_interAssocWith_correct :: [(Int,Int)] -> [(Int,Int)] -> Bool
prop_interAssocWith_correct xs ys =
interAssocWith (*) l l' == interAssocWith' (*) l l'
where
l = List.sortBy (compare `on` fst) xs
l' = List.sortBy (compare `on` fst) ys
prop_mul :: Size Int -> Property
prop_mul sz =
mapSize (`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.intSemiring
tests :: IO Bool
tests = runTests "Agda.Termination.SparseMatrix"
[ 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_toSparseRows
, quickCheck' prop_fromLists_toLists
, quickCheck' prop_isSingleton
, quickCheck' prop_zipMatrices_correct
, quickCheck' prop_add
, quickCheck' prop_add_correct
, quickCheck' prop_mul
, quickCheck' prop_diagonal
, quickCheck' prop_addColumn
, quickCheck' prop_addRow
]