diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,3 +1,45 @@ -# Changelog for decision-diagrams +# Changelog for `decision-diagrams` package -## Unreleased changes +## 0.2.0.0 + +### Changes + +* Make `Leaf :: Bool -> BDD a` as a basic constructor instead of + `F`/`T` (in case of BDD) and `Empty`/`Base` (in case of ZDD), and + remove `F`/`T`. + +* `ZDD.toList` now returns sorted list + +* Change signature of `fold` and `fold'` of BDD + * Before: `b -> b -> (Int -> b -> b -> b) -> BDD a -> b` + * After: `(Int -> b -> b -> b) -> (Bool -> b) -> BDD a -> b` + +* Change signature of `fold` and `fold'` of ZDD (ditto) + +* Add `HasCallStack` to some functions that are expected to raise excpetions + +### Additions + +* Introduce signature functor type (`Sig`) + +* Add new operations: + * BDD: + * fixed point operators `lfp` and `gfp` + * satisfiability related functions: `anySat`, `allSat`, `anySatComplete`, `allSatComplete`, `countSat`, `uniformSatM` + * pseudo-boolean constraint functions: `pbAtLeast`, `pbAtMost`, `pbExactly`, `pbExactlyIntegral` + * ZDD: + * `combinations` + * pseudo-boolean constraint functions: `subsetsAtLeast`, `subsetsAtMost`, `subsetsExactly`, `subsetsExactlyIntegral` + * Both BDD and ZDD + * `numNodes` + * `unfoldHashable` and `unfoldOrd` + +### Bug fixes + +* Fix laziness of `fold` and `fold'` + +### Other changes + +* Introduced `doctest` + +* Add other-extensions fields to `package.yaml` diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,8 +1,31 @@ # decision-diagrams +Hackage: +[![Hackage](https://img.shields.io/hackage/v/decision-diagrams.svg)](https://hackage.haskell.org/package/decision-diagrams) +[![Hackage Deps](https://img.shields.io/hackage-deps/v/decision-diagrams.svg)](https://packdeps.haskellers.com/feed?needle=decision-diagrams) +[![Hackage CI](https://matrix.hackage.haskell.org/api/v2/packages/decision-diagrams/badge)](https://matrix.hackage.haskell.org/#/package/decision-diagrams) + +Dev: [![Build Status](https://github.com/msakai/haskell-decision-diagrams/actions/workflows/build.yaml/badge.svg)](https://github.com/msakai/haskell-decision-diagrams/actions/workflows/build.yaml) [![Coverage Status](https://coveralls.io/repos/msakai/haskell-decision-diagrams/badge.svg)](https://coveralls.io/r/msakai/haskell-decision-diagrams) -Binary Decision Diagrams (BDD) and Zero-suppressed Binary Decision Diagrams (ZDD) implementation in Haskell. +[Binary Decision Diagrams (BDD)](https://en.wikipedia.org/wiki/Binary_decision_diagram) and [Zero-suppressed Binary Decision Diagrams (ZDD)](https://en.wikipedia.org/wiki/Zero-suppressed_decision_diagram) implementation in Haskell. -Hash-consing is implemented using [intern](https://hackage.haskell.org/package/intern) package. +BDD is a data stucture suitable for representing boolean functions (can be thought as a compressed representation of truth tables) and many operations on boolean functions can be performed efficiently. ZDD is a variant of BDD suitable for representing (sparse) families of sets compactly. + +BDD/ZDD uses hash-consing for compact representation and efficient comparison, and we use [intern](https://hackage.haskell.org/package/intern) package for implementing hash-consing. + +## Comparison with other BDD packages for Haskell + +|Package name|Repository|BDD|ZDD|Style|Implementation|Hash-consing / Fast equality test| +|------------|----------|---|---|-----|--------------|---------------------------------| +|[decision-diagrams](https://hackage.haskell.org/package/decision-diagrams) (this package)|[GitHub](https://github.com/msakai/haskell-decision-diagrams/)|✔️|✔️|pure|pure Haskell|✔️| +|[zsdd](https://hackage.haskell.org/package/zsdd)|[GitHub](https://github.com/eddiejones2108/decision-diagrams) (deleted?)|✔️|✔️|monadic|pure Haskell|✔️| +|[obdd](https://hackage.haskell.org/package/obdd)|[GitHub](https://github.com/jwaldmann/haskell-obdd)|✔️|-|pure|pure Haskell|-| +|[HasCacBDD](https://hackage.haskell.org/package/HasCacBDD)|[GitHub](https://github.com/m4lvin/HasCacBDD)|✔️|-|pure|FFI|✔️| +|[hBDD](https://hackage.haskell.org/package/hBDD) ([hBDD-CUDD](https://hackage.haskell.org/package/hBDD-CUDD), [hBDD-CMUBDD](https://hackage.haskell.org/package/hBDD-CMUBDD))|[GitHub](https://github.com/peteg/hBDD)|✔️|-|pure|FFI|✔️| +|[cudd](https://hackage.haskell.org/package/cudd)|[GitHub](https://github.com/adamwalker/haskell_cudd)|✔️|-|both pure\*1 and monadic|FFI|✔️| + +\*1: `cudd`'s pure interface is different from normal Haskell data types (like ones in the [containers](https://hackage.haskell.org/package/containers) package, for example) because it requires `DDManager` argument. + +Please feel free to make a pull request for addition or correction to the comparision. diff --git a/decision-diagrams.cabal b/decision-diagrams.cabal --- a/decision-diagrams.cabal +++ b/decision-diagrams.cabal @@ -5,10 +5,10 @@ -- see: https://github.com/sol/hpack name: decision-diagrams -version: 0.1.0.0 +version: 0.2.0.0 synopsis: Binary Decision Diagrams (BDD) and Zero-suppressed Binary Decision Diagrams (ZDD) description: Please see the README on GitHub at -category: Data, Logic +category: Data, Data Structures, Logic homepage: https://github.com/msakai/haskell-decision-diagrams#readme bug-reports: https://github.com/msakai/haskell-decision-diagrams/issues author: Masahiro Sakai @@ -34,19 +34,54 @@ Data.DecisionDiagram.BDD.Internal.Node hs-source-dirs: src + other-extensions: + BangPatterns + FlexibleContexts + FlexibleInstances + RankNTypes + ScopedTypeVariables + CPP + DeriveGeneric + DeriveTraversable + GeneralizedNewtypeDeriving + PatternSynonyms + TypeFamilies + UndecidableInstances + ViewPatterns build-depends: - base >=4.7 && <5 + base >=4.11.0.0 && <5 , containers >=0.5.11.0 && <0.7 - , hashable >=1.2.7.0 && <1.4 - , hashtables >=1.2.3.1 && <1.3 + , hashable >=1.2.7.0 && <1.5 + , hashtables >=1.2.3.1 && <1.4 , intern >=0.9.1.2 && <1.0.0.0 , mwc-random >=0.13.6.0 && <0.16 , primitive >=0.6.3.0 && <0.8 , random >=1.1 && <1.3 , reflection >=2.1.4 && <2.2 , unordered-containers >=0.2.9.0 && <0.3 + , vector >=0.12.0.2 && <0.13 default-language: Haskell2010 +test-suite decision-diagrams-doctest + type: exitcode-stdio-1.0 + main-is: test/doctests.hs + other-modules: + Paths_decision_diagrams + other-extensions: + BangPatterns + FlexibleContexts + FlexibleInstances + RankNTypes + ScopedTypeVariables + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + base >=4.11.0.0 && <5 + , containers >=0.5.11.0 && <0.7 + , decision-diagrams + , doctest + , mwc-random >=0.13.6.0 && <0.16 + default-language: Haskell2010 + test-suite decision-diagrams-test type: exitcode-stdio-1.0 main-is: TestSuite.hs @@ -57,16 +92,26 @@ Paths_decision_diagrams hs-source-dirs: test + other-extensions: + BangPatterns + FlexibleContexts + FlexibleInstances + RankNTypes + ScopedTypeVariables + TemplateHaskell ghc-options: -threaded -rtsopts -with-rtsopts=-N build-depends: QuickCheck >=2.11.3 && <2.15 - , base >=4.7 && <5 + , base >=4.11.0.0 && <5 , containers >=0.5.11.0 && <0.7 , decision-diagrams + , deepseq >=1.4.3.0 && <1.5 , mwc-random >=0.13.6.0 && <0.16 + , quickcheck-instances >=0.3.19 && <0.4 , statistics >=0.14.0.2 && <0.16 , tasty >=1.1.0.4 && <1.5 , tasty-hunit >=0.10.0.1 && <0.11 , tasty-quickcheck ==0.10.* , tasty-th >=0.1.7 && <0.2 + , vector >=0.12.0.2 && <0.13 default-language: Haskell2010 diff --git a/src/Data/DecisionDiagram/BDD.hs b/src/Data/DecisionDiagram/BDD.hs --- a/src/Data/DecisionDiagram/BDD.hs +++ b/src/Data/DecisionDiagram/BDD.hs @@ -1,5 +1,6 @@ {-# OPTIONS_GHC -Wall #-} {-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE PatternSynonyms #-} {-# LANGUAGE RankNTypes #-} @@ -28,7 +29,7 @@ module Data.DecisionDiagram.BDD ( -- * The BDD type - BDD (F, T, Branch) + BDD (Leaf, Branch) -- * Item ordering , ItemOrder (..) @@ -53,6 +54,12 @@ , andB , orB + -- * Pseudo-boolean constraints + , pbAtLeast + , pbAtMost + , pbExactly + , pbExactlyIntegral + -- * Quantification , forAll , exists @@ -64,6 +71,7 @@ -- * Query , support , evaluate + , numNodes -- * Restriction / Cofactor , restrict @@ -74,13 +82,33 @@ , subst , substSet + -- * Satisfiability + , anySat + , allSat + , anySatComplete + , allSatComplete + , countSat + , uniformSatM + + -- * (Co)algebraic structure + , Sig (..) + , inSig + , outSig + -- * Fold , fold , fold' + -- * Unfold + , unfoldHashable + , unfoldOrd + + -- * Fixpoints + , lfp + , gfp + -- * Conversion from/to graphs , Graph - , Node (..) , toGraph , toGraph' , fromGraph @@ -89,10 +117,16 @@ import Control.Exception (assert) import Control.Monad +#if !MIN_VERSION_mwc_random(0,15,0) +import Control.Monad.Primitive +#endif import Control.Monad.ST +import Control.Monad.ST.Unsafe +import Data.Bits (Bits (shiftL)) +import qualified Data.Foldable as Foldable import Data.Function (on) -import Data.Functor.Identity import Data.Hashable +import qualified Data.HashMap.Lazy as HashMap import qualified Data.HashTable.Class as H import qualified Data.HashTable.ST.Cuckoo as C import Data.IntMap (IntMap) @@ -100,11 +134,24 @@ import Data.IntSet (IntSet) import qualified Data.IntSet as IntSet import Data.List (sortBy) +import Data.Map.Lazy (Map) +import qualified Data.Map.Lazy as Map import Data.Proxy -import Data.STRef +import Data.Ratio +import qualified Data.Vector as V +import GHC.Stack +import Numeric.Natural +#if MIN_VERSION_mwc_random(0,15,0) +import System.Random.MWC (Uniform (..)) +import System.Random.Stateful (StatefulGen (..)) +#else +import System.Random.MWC (Gen, Variate (..)) +#endif +import System.Random.MWC.Distributions (bernoulli) import Text.Read import Data.DecisionDiagram.BDD.Internal.ItemOrder +import Data.DecisionDiagram.BDD.Internal.Node (Sig (..), Graph) import qualified Data.DecisionDiagram.BDD.Internal.Node as Node infixr 3 .&&. @@ -124,11 +171,14 @@ deriving (Eq, Hashable) pattern F :: BDD a -pattern F = BDD Node.F +pattern F = Leaf False pattern T :: BDD a -pattern T = BDD Node.T +pattern T = Leaf True +pattern Leaf :: Bool -> BDD a +pattern Leaf b = BDD (Node.Leaf b) + -- | Smart constructor that takes the BDD reduction rules into account pattern Branch :: Int -> BDD a -> BDD a -> BDD a pattern Branch x lo hi <- BDD (Node.Branch x (BDD -> lo) (BDD -> hi)) where @@ -137,6 +187,7 @@ | otherwise = BDD (Node.Branch x lo hi) {-# COMPLETE T, F, Branch #-} +{-# COMPLETE Leaf, Branch #-} nodeId :: BDD a -> Int nodeId (BDD node) = Node.nodeId node @@ -155,15 +206,11 @@ EQ -> BDDCase2EQ ptop p0 p1 q0 q1 bddCase2 _ (Branch ptop p0 p1) _ = BDDCase2LT ptop p0 p1 bddCase2 _ _ (Branch qtop q0 q1) = BDDCase2GT qtop q0 q1 -bddCase2 _ T T = BDDCase2EQ2 True True -bddCase2 _ T F = BDDCase2EQ2 True False -bddCase2 _ F T = BDDCase2EQ2 False True -bddCase2 _ F F = BDDCase2EQ2 False False +bddCase2 _ (Leaf b1) (Leaf b2) = BDDCase2EQ2 b1 b2 level :: BDD a -> Level a -level T = Terminal -level F = Terminal level (Branch x _ _) = NonTerminal x +level (Leaf _) = Terminal -- ------------------------------------------------------------------------ @@ -197,8 +244,7 @@ notB :: BDD a -> BDD a notB bdd = runST $ do h <- C.newSized defaultTableSize - let f T = return F - f F = return T + let f (Leaf b) = return (Leaf (not b)) f n@(Branch ind lo hi) = do m <- H.lookup h n case m of @@ -293,10 +339,7 @@ (.<=>.) :: forall a. ItemOrder a => BDD a -> BDD a -> BDD a (.<=>.) = apply True f where - f T T = Just T - f T F = Just F - f F T = Just F - f F F = Just T + f (Leaf b1) (Leaf b2) = Just (Leaf (b1 == b2)) f a b | a == b = Just T f _ _ = Nothing @@ -304,8 +347,7 @@ ite :: forall a. ItemOrder a => BDD a -> BDD a -> BDD a -> BDD a ite c' t' e' = runST $ do h <- C.newSized defaultTableSize - let f T t _ = return t - f F _ e = return e + let f (Leaf b) t e = if b then return t else return e f _ t e | t == e = return t f c t e = do case minimum [level c, level t, level e] of @@ -338,6 +380,88 @@ -- ------------------------------------------------------------------------ +-- | Pseudo-boolean constraint, /w1*x1 + w2*x2 + … ≥ k/. +pbAtLeast :: forall a w. (ItemOrder a, Real w) => IntMap w -> w -> BDD a +pbAtLeast xs k0 = unfoldOrd f (0, k0) + where + xs' :: V.Vector (Int, w) + xs' = V.fromList $ sortBy (compareItem (Proxy :: Proxy a) `on` fst) $ IntMap.toList xs + ys :: V.Vector (w, w) + ys = V.scanr (\(_, w) (lb,ub) -> if w >= 0 then (lb, ub+w) else (lb+w, ub)) (0,0) xs' + + f :: (Int, w) -> Sig (Int, w) + f (!i, !k) + | not (k <= ub) = SLeaf False + | i == V.length xs' && 0 >= k = SLeaf True + | lb >= k = SLeaf True -- all remaining variables are don't-care + | otherwise = SBranch x (i+1, k) (i+1, k-w) + where + (lb,ub) = ys V.! i + (x, w) = xs' V.! i + +-- | Pseudo-boolean constraint, /w1*x1 + w2*x2 + … ≤ k/. +pbAtMost :: forall a w. (ItemOrder a, Real w) => IntMap w -> w -> BDD a +pbAtMost xs k0 = unfoldOrd f (0, k0) + where + xs' :: V.Vector (Int, w) + xs' = V.fromList $ sortBy (compareItem (Proxy :: Proxy a) `on` fst) $ IntMap.toList xs + ys :: V.Vector (w, w) + ys = V.scanr (\(_, w) (lb,ub) -> if w >= 0 then (lb, ub+w) else (lb+w, ub)) (0,0) xs' + + f :: (Int, w) -> Sig (Int, w) + f (!i, !k) + | not (lb <= k) = SLeaf False + | i == V.length xs' && 0 <= k = SLeaf True + | ub <= k = SLeaf True -- all remaining variables are don't-care + | otherwise = SBranch x (i+1, k) (i+1, k-w) + where + (lb,ub) = ys V.! i + (x, w) = xs' V.! i + +-- | Pseudo-boolean constraint, /w1*x1 + w2*x2 + … = k/. +-- +-- If weight type is 'Integral', 'pbExactlyIntegral' is more efficient. +pbExactly :: forall a w. (ItemOrder a, Real w) => IntMap w -> w -> BDD a +pbExactly xs k0 = unfoldOrd f (0, k0) + where + xs' :: V.Vector (Int, w) + xs' = V.fromList $ sortBy (compareItem (Proxy :: Proxy a) `on` fst) $ IntMap.toList xs + ys :: V.Vector (w, w) + ys = V.scanr (\(_, w) (lb,ub) -> if w >= 0 then (lb, ub+w) else (lb+w, ub)) (0,0) xs' + + f :: (Int, w) -> Sig (Int, w) + f (!i, !k) + | not (lb <= k && k <= ub) = SLeaf False + | i == V.length xs' && 0 == k = SLeaf True + | otherwise = SBranch x (i+1, k) (i+1, k-w) + where + (lb,ub) = ys V.! i + (x, w) = xs' V.! i + +-- | Similar to 'pbExactly' but more efficient. +pbExactlyIntegral :: forall a w. (ItemOrder a, Real w, Integral w) => IntMap w -> w -> BDD a +pbExactlyIntegral xs k0 = unfoldOrd f (0, k0) + where + xs' :: V.Vector (Int, w) + xs' = V.fromList $ sortBy (compareItem (Proxy :: Proxy a) `on` fst) $ IntMap.toList xs + ys :: V.Vector (w, w) + ys = V.scanr (\(_, w) (lb,ub) -> if w >= 0 then (lb, ub+w) else (lb+w, ub)) (0,0) xs' + ds :: V.Vector w + ds = V.scanr1 (\w d -> if w /= 0 then gcd w d else d) (V.map snd xs') + + f :: (Int, w) -> Sig (Int, w) + f (!i, !k) + | not (lb <= k && k <= ub) = SLeaf False + | i == V.length xs' && 0 == k = SLeaf True + | d /= 0 && k `mod` d /= 0 = SLeaf False + | otherwise = SBranch x (i+1, k) (i+1, k-w) + where + (lb,ub) = ys V.! i + (x, w) = xs' V.! i + d = ds V.! i + +-- ------------------------------------------------------------------------ + -- | Universal quantification (∀) forAll :: forall a. ItemOrder a => Int -> BDD a -> BDD a forAll x bdd = runST $ do @@ -465,48 +589,59 @@ -- | Fold over the graph structure of the BDD. -- --- It takes values for substituting 'false' ('F') and 'true' ('T'), --- and a function for substiting non-terminal node ('Branch'). -fold :: b -> b -> (Int -> b -> b -> b) -> BDD a -> b -fold ff tt br bdd = runST $ do - h <- C.newSized defaultTableSize - let f F = return ff - f T = return tt - f p@(Branch top lo hi) = do - m <- H.lookup h p - case m of - Just ret -> return ret - Nothing -> do - r0 <- f lo - r1 <- f hi - let ret = br top r0 r1 - H.insert h p ret - return ret - f bdd +-- It takes two functions that substitute 'Branch' and 'Leaf' respectively. +-- +-- Note that its type is isomorphic to @('Sig' b -> b) -> BDD a -> b@. +fold :: (Int -> b -> b -> b) -> (Bool -> b) -> BDD a -> b +fold br lf (BDD node) = Node.fold br lf node -- | Strict version of 'fold' -fold' :: b -> b -> (Int -> b -> b -> b) -> BDD a -> b -fold' ff tt br bdd = runST $ do - op <- mkFold'Op ff tt br - op bdd +fold' :: (Int -> b -> b -> b) -> (Bool -> b) -> BDD a -> b +fold' br lf (BDD node) = Node.fold' br lf node -mkFold'Op :: b -> b -> (Int -> b -> b -> b) -> ST s (BDD a -> ST s b) -mkFold'Op !ff !tt br = do +mkFold'Op :: (Int -> b -> b -> b) -> (Bool -> b) -> ST s (BDD a -> ST s b) +mkFold'Op br lf = do + op <- Node.mkFold'Op br lf + return $ \(BDD node) -> op node + +-- ------------------------------------------------------------------------ + +-- | Top-down construction of BDD, memoising internal states using 'Hashable' instance. +unfoldHashable :: forall a b. (ItemOrder a, Eq b, Hashable b) => (b -> Sig b) -> b -> BDD a +unfoldHashable f b = runST $ do h <- C.newSized defaultTableSize - let f F = return ff - f T = return tt - f p@(Branch top lo hi) = do - m <- H.lookup h p - case m of - Just ret -> return ret + let g [] = return () + g (x : xs) = do + r <- H.lookup h x + case r of + Just _ -> g xs Nothing -> do - r0 <- f lo - r1 <- f hi - let ret = br top r0 r1 - seq ret $ H.insert h p ret - return ret - return f + let fx = f x + H.insert h x fx + g (xs ++ Foldable.toList fx) + g [b] + xs <- H.toList h + let h2 = HashMap.fromList [(x, inSig (fmap (h2 HashMap.!) s)) | (x,s) <- xs] + return $ h2 HashMap.! b +-- | Top-down construction of BDD, memoising internal states using 'Ord' instance. +unfoldOrd :: forall a b. (ItemOrder a, Ord b) => (b -> Sig b) -> b -> BDD a +unfoldOrd f b = m2 Map.! b + where + m1 :: Map b (Sig b) + m1 = g Map.empty [b] + + m2 :: Map b (BDD a) + m2 = Map.map (inSig . fmap (m2 Map.!)) m1 + + g m [] = m + g m (x : xs) = + case Map.lookup x m of + Just _ -> g m xs + Nothing -> + let fx = f x + in g (Map.insert x fx m) (xs ++ Foldable.toList fx) + -- ------------------------------------------------------------------------ -- | All the variables that this BDD depends on. @@ -516,9 +651,10 @@ op bdd mkSupportOp :: ST s (BDD a -> ST s IntSet) -mkSupportOp = mkFold'Op IntSet.empty IntSet.empty f +mkSupportOp = mkFold'Op f g where f x lo hi = IntSet.insert x (lo `IntSet.union` hi) + g _ = IntSet.empty -- | Evaluate a boolean function represented as BDD under the valuation -- given by @(Int -> Bool)@, i.e. it lifts a valuation function from @@ -526,20 +662,31 @@ evaluate :: (Int -> Bool) -> BDD a -> Bool evaluate f = g where - g F = False - g T = True + g (Leaf b) = b g (Branch x lo hi) | f x = g hi | otherwise = g lo +-- | Count the number of nodes in a BDD viewed as a rooted directed acyclic graph. +-- +-- See also 'toGraph'. +numNodes :: BDD a -> Int +numNodes (BDD node) = Node.numNodes node + -- ------------------------------------------------------------------------ --- | Compute \(F_x \) or \(F_{\neg x} \). +-- | Compute \(F|_{x_i} \) or \(F|_{\neg x_i} \). +-- +-- \[ +-- F|_{x_i}(\ldots, x_{i-1}, x_{i+1}, \ldots) = F(\ldots, x_{i-1}, \mathrm{True}, x_{i+1}, \ldots) +-- \] +-- \[ +-- F|_{\neg x_i}(\ldots, x_{i-1}, x_{i+1}, \ldots) = F(\ldots, x_{i-1}, \mathrm{False}, x_{i+1}, \ldots) +-- \] restrict :: forall a. ItemOrder a => Int -> Bool -> BDD a -> BDD a restrict x val bdd = runST $ do h <- C.newSized defaultTableSize - let f T = return T - f F = return F + let f n@(Leaf _) = return n f n@(Branch ind lo hi) = do m <- H.lookup h n case m of @@ -553,13 +700,12 @@ return ret f bdd --- | Compute \(F_{\{x_i = v_i\}_i} \). +-- | Compute \(F|_{\{x_i = v_i\}_i} \). restrictSet :: forall a. ItemOrder a => IntMap Bool -> BDD a -> BDD a restrictSet val bdd = runST $ do h <- C.newSized defaultTableSize let f [] n = return n - f _ T = return T - f _ F = return F + f _ n@(Leaf _) = return n f xxs@((x,v) : xs) n@(Branch ind lo hi) = do m <- H.lookup h n case m of @@ -581,8 +727,7 @@ h <- C.newSized defaultTableSize let f T n = return n f F _ = return T -- Is this correct? - f _ F = return F - f _ T = return T + f _ n@(Leaf _) = return n f n1 n2 | n1 == n2 = return T f n1 n2 = do m <- H.lookup h (n1, n2) @@ -606,8 +751,9 @@ -- ------------------------------------------------------------------------ --- | @subst x N M@ computes substitution \(M_{x = N}\). +-- | @subst x N M@ computes substitution M[x ↦ N]. -- +-- Note the order of the arguments. -- This operation is also known as /Composition/. subst :: forall a. ItemOrder a => Int -> BDD a -> BDD a -> BDD a subst x n m = runST $ do @@ -643,7 +789,9 @@ return ret f m m n --- | Simultaneous substitution +-- | Simultaneous substitution. +-- +-- Note that this is not the same as repeated application of 'subst'. substSet :: forall a. ItemOrder a => IntMap (BDD a) -> BDD a -> BDD a substSet s m = runST $ do supportOp <- mkSupportOp @@ -705,87 +853,220 @@ fixed = IntMap.mapMaybe asBool conditions asBool :: BDD a -> Maybe Bool - asBool a = - case a of - T -> Just True - F -> Just False - _ -> Nothing + asBool (Leaf b) = Just b + asBool _ = Nothing -- ------------------------------------------------------------------------ -type Graph = IntMap Node +-- | Least fixed point. +-- +-- Monotonicity of the operator is assumed but not checked. +lfp :: ItemOrder a => (BDD a -> BDD a) -> BDD a +lfp f = go false + where + go curr + | curr == next = curr + | otherwise = go next + where + next = f curr -data Node - = NodeF - | NodeT - | NodeBranch !Int Int Int - deriving (Eq, Show, Read) +-- | Greatest fixed point. +-- +-- Monotonicity of the operator is assumed but not checked. +gfp :: ItemOrder a => (BDD a -> BDD a) -> BDD a +gfp f = go true + where + go curr + | curr == next = curr + | otherwise = go next + where + next = f curr --- | Convert a BDD into a pointed graph -toGraph :: BDD a -> (Graph, Int) -toGraph bdd = - case toGraph' (Identity bdd) of - (g, Identity v) -> (g, v) +-- ------------------------------------------------------------------------ --- | Convert multiple BDDs into a graph -toGraph' :: Traversable t => t (BDD a) -> (Graph, t Int) -toGraph' bs = runST $ do +findSatM :: MonadPlus m => BDD a -> m (IntMap Bool) +findSatM = fold f g + where + f x lo hi = mplus (liftM (IntMap.insert x False) lo) (liftM (IntMap.insert x True) hi) + g b = if b then return IntMap.empty else mzero + +-- | Find one satisfying partial assignment +anySat :: BDD a -> Maybe (IntMap Bool) +anySat = findSatM + +-- | Enumerate all satisfying partial assignments +allSat :: BDD a -> [IntMap Bool] +allSat = findSatM + +findSatCompleteM :: forall a m. (MonadPlus m, ItemOrder a, HasCallStack) => IntSet -> BDD a -> m (IntMap Bool) +findSatCompleteM xs0 bdd = runST $ do h <- C.newSized defaultTableSize - H.insert h F 0 - H.insert h T 1 - counter <- newSTRef 2 - ref <- newSTRef $ IntMap.fromList [(0, NodeF), (1, NodeT)] + let f _ (Leaf False) = return $ mzero + f xs (Leaf True) = return $ foldM (\m x -> msum [return (IntMap.insert x v m) | v <- [False, True]]) IntMap.empty xs + f xs n@(Branch x lo hi) = do + case span (\x2 -> compareItem (Proxy :: Proxy a) x2 x == LT) xs of + (ys, (x':xs')) | x == x' -> do + r <- H.lookup h n + ps <- case r of + Just ret -> return ret + Nothing -> do + r0 <- f xs' lo + r1 <- unsafeInterleaveST $ f xs' hi + let ret = liftM (IntMap.insert x False) r0 `mplus` liftM (IntMap.insert x True) r1 + H.insert h n ret + return ret + return $ do + p <- ps + foldM (\m y -> msum [return (IntMap.insert y v m) | v <- [False, True]]) p ys + _ -> error ("findSatCompleteM: " ++ show x ++ " should not occur") + f (sortBy (compareItem (Proxy :: Proxy a)) (IntSet.toList xs0)) bdd - let f F = return 0 - f T = return 1 - f p@(Branch x lo hi) = do - m <- H.lookup h p - case m of - Just ret -> return ret - Nothing -> do - r0 <- f lo - r1 <- f hi - n <- readSTRef counter - writeSTRef counter $! n+1 - H.insert h p n - modifySTRef' ref (IntMap.insert n (NodeBranch x r0 r1)) - return n +-- | Find one satisfying (complete) assignment over a given set of variables +-- +-- The set of variables must be a superset of 'support'. +anySatComplete :: ItemOrder a => IntSet -> BDD a -> Maybe (IntMap Bool) +anySatComplete = findSatCompleteM - vs <- mapM f bs - g <- readSTRef ref - return (g, vs) +-- | Enumerate all satisfying (complete) assignment over a given set of variables +-- +-- The set of variables must be a superset of 'support'. +allSatComplete :: ItemOrder a => IntSet -> BDD a -> [IntMap Bool] +allSatComplete = findSatCompleteM --- | Convert a pointed graph into a BDD -fromGraph :: (Graph, Int) -> BDD a -fromGraph (g, v) = - case IntMap.lookup v (fromGraph' g) of - Nothing -> error ("Data.DecisionDiagram.BDD.fromGraph: invalid node id " ++ show v) - Just bdd -> bdd +{-# SPECIALIZE countSat :: ItemOrder a => IntSet -> BDD a -> Int #-} +{-# SPECIALIZE countSat :: ItemOrder a => IntSet -> BDD a -> Integer #-} +{-# SPECIALIZE countSat :: ItemOrder a => IntSet -> BDD a -> Natural #-} +-- | Count the number of satisfying (complete) assignment over a given set of variables. +-- +-- The set of variables must be a superset of 'support'. +-- +-- It is polymorphic in return type, but it is recommended to use 'Integer' or 'Natural' +-- because the size can be larger than fixed integer types such as @Int64@. +-- +-- >>> countSat (IntSet.fromList [1..128]) (true :: BDD AscOrder) +-- 340282366920938463463374607431768211456 +-- >>> import Data.Int +-- >>> maxBound :: Int64 +-- 9223372036854775807 +countSat :: forall a b. (ItemOrder a, Num b, Bits b, HasCallStack) => IntSet -> BDD a -> b +countSat xs bdd = runST $ do + h <- C.newSized defaultTableSize + let f _ (Leaf False) = return $ 0 + f ys (Leaf True) = return $! 1 `shiftL` length ys + f ys node@(Branch x lo hi) = do + case span (\x2 -> compareItem (Proxy :: Proxy a) x2 x == LT) ys of + (zs, y' : ys') | x == y' -> do + m <- H.lookup h node + n <- case m of + Just n -> return n + Nothing -> do + n <- liftM2 (+) (f ys' lo) (f ys' hi) + H.insert h node n + return n + return $! n `shiftL` length zs + (_, _) -> error ("countSat: " ++ show x ++ " should not occur") + f (sortBy (compareItem (Proxy :: Proxy a)) (IntSet.toList xs)) bdd --- | Convert nodes of a graph into BDDs -fromGraph' :: Graph -> IntMap (BDD a) -fromGraph' g = ret +-- | Sample an assignment from uniform distribution over complete satisfiable assignments ('allSatComplete') of the BDD. +-- +-- The function constructs a table internally and the table is shared across +-- multiple use of the resulting action (@m IntSet@). +-- Therefore, the code +-- +-- @ +-- let g = uniformSatM xs bdd gen +-- s1 <- g +-- s2 <- g +-- @ +-- +-- is more efficient than +-- +-- @ +-- s1 <- uniformSatM xs bdd gen +-- s2 <- uniformSatM xs bdd gen +-- @ +-- . +#if MIN_VERSION_mwc_random(0,15,0) +uniformSatM :: forall a g m. (ItemOrder a, StatefulGen g m, HasCallStack) => IntSet -> BDD a -> g -> m (IntMap Bool) +#else +uniformSatM :: forall a m. (ItemOrder a, PrimMonad m, HasCallStack) => IntSet -> BDD a -> Gen (PrimState m) -> m (IntMap Bool) +#endif +uniformSatM xs0 bdd0 = func IntMap.empty where - ret = IntMap.map f g - f NodeF = F - f NodeT = T - f (NodeBranch x lo hi) = - case (IntMap.lookup lo ret, IntMap.lookup hi ret) of - (Nothing, _) -> error ("Data.DecisionDiagram.BDD.fromGraph': invalid node id " ++ show lo) - (_, Nothing) -> error ("Data.DecisionDiagram.BDD.fromGraph': invalid node id " ++ show hi) - (Just lo', Just hi') -> Branch x lo' hi' + func = runST $ do + h <- C.newSized defaultTableSize + let f xs bdd = + case span (\x2 -> NonTerminal x2 < level bdd) xs of + (ys, xxs') -> do + xs' <- case (bdd, xxs') of + (Branch x _ _, x' : xs') | x == x' -> return xs' + (Branch x _ _, _) -> error ("uniformSatM: " ++ show x ++ " should not occur") + (Leaf _, []) -> return [] + (Leaf _, _:_) -> error ("uniformSatM: should not happen") + (s, func0) <- g xs' bdd + let func' !m !gen = do +#if MIN_VERSION_mwc_random(0,15,0) + vals <- replicateM (length ys) (uniformM gen) +#else + vals <- replicateM (length ys) (uniform gen) +#endif + func0 (m `IntMap.union` IntMap.fromList (zip ys vals)) gen + return (s `shiftL` length ys, func') + g _ (Leaf True) = return (1 :: Integer, \a _gen -> return a) + g _ (Leaf False) = return (0 :: Integer, \_a _gen -> error "uniformSatM: should not happen") + g xs bdd@(Branch x lo hi) = do + m <- H.lookup h bdd + case m of + Just ret -> return ret + Nothing -> do + (n0, func0) <- f xs lo + (n1, func1) <- f xs hi + let s = n0 + n1 + r :: Double + r = realToFrac (n1 % s) + seq r $ return () + let func' !a !gen = do + b <- bernoulli r gen + if b then + func1 (IntMap.insert x True a) gen + else + func0 (IntMap.insert x False a) gen + H.insert h bdd (s, func') + return (s, func') + liftM snd $ f (sortBy (compareItem (Proxy :: Proxy a)) (IntSet.toList xs0)) bdd0 -- ------------------------------------------------------------------------ --- https://ja.wikipedia.org/wiki/%E4%BA%8C%E5%88%86%E6%B1%BA%E5%AE%9A%E5%9B%B3 -_test_bdd :: BDD AscOrder -_test_bdd = (notB x1 .&&. notB x2 .&&. notB x3) .||. (x1 .&&. x2) .||. (x2 .&&. x3) - where - x1 = var 1 - x2 = var 2 - x3 = var 3 -{- -BDD (Node 880 (UBranch 1 (Node 611 (UBranch 2 (Node 836 UT) (Node 215 UF))) (Node 806 (UBranch 2 (Node 842 (UBranch 3 (Node 836 UT) (Node 215 UF))) (Node 464 (UBranch 3 (Node 215 UF) (Node 836 UT))))))) --} +-- | 'Sig'-algebra stucture of 'BDD', \(\mathrm{in}_\mathrm{Sig}\). +inSig :: Sig (BDD a) -> BDD a +inSig (SLeaf b) = Leaf b +inSig (SBranch x lo hi) = Branch x lo hi + +-- | 'Sig'-coalgebra stucture of 'BDD', \(\mathrm{out}_\mathrm{Sig}\). +outSig :: BDD a -> Sig (BDD a) +outSig (Leaf b) = SLeaf b +outSig (Branch x lo hi) = SBranch x lo hi + +-- ------------------------------------------------------------------------ + +-- | Convert a BDD into a pointed graph +-- +-- Nodes @0@ and @1@ are reserved for @SLeaf False@ and @SLeaf True@ +-- even if they are not actually used. Therefore the result may be +-- larger than 'numNodes' if the leaf nodes are not used. +toGraph :: BDD a -> (Graph Sig, Int) +toGraph (BDD node) = Node.toGraph node + +-- | Convert multiple BDDs into a graph +toGraph' :: Traversable t => t (BDD a) -> (Graph Sig, t Int) +toGraph' bs = Node.toGraph' (fmap (\(BDD node) -> node) bs) + +-- | Convert a pointed graph into a BDD +fromGraph :: HasCallStack => (Graph Sig, Int) -> BDD a +fromGraph = Node.foldGraph inSig + +-- | Convert nodes of a graph into BDDs +fromGraph' :: HasCallStack => Graph Sig -> IntMap (BDD a) +fromGraph' = Node.foldGraphNodes inSig -- ------------------------------------------------------------------------ diff --git a/src/Data/DecisionDiagram/BDD/Internal/ItemOrder.hs b/src/Data/DecisionDiagram/BDD/Internal/ItemOrder.hs --- a/src/Data/DecisionDiagram/BDD/Internal/ItemOrder.hs +++ b/src/Data/DecisionDiagram/BDD/Internal/ItemOrder.hs @@ -25,6 +25,9 @@ , withDescOrder , withCustomOrder + -- * Ordered item + , OrderedItem (..) + -- * Level , Level (..) ) where @@ -64,6 +67,14 @@ withCustomOrder :: forall r. (Int -> Int -> Ordering) -> (forall a. ItemOrder a => Proxy a -> r) -> r withCustomOrder cmp k = reify cmp (\(_ :: Proxy s) -> k (Proxy :: Proxy (CustomOrder s))) + +-- ------------------------------------------------------------------------ + +newtype OrderedItem a = OrderedItem Int + deriving (Eq, Show) + +instance ItemOrder a => Ord (OrderedItem a) where + compare (OrderedItem x) (OrderedItem y) = compareItem (Proxy :: Proxy a) x y -- ------------------------------------------------------------------------ diff --git a/src/Data/DecisionDiagram/BDD/Internal/Node.hs b/src/Data/DecisionDiagram/BDD/Internal/Node.hs --- a/src/Data/DecisionDiagram/BDD/Internal/Node.hs +++ b/src/Data/DecisionDiagram/BDD/Internal/Node.hs @@ -1,5 +1,6 @@ {-# OPTIONS_GHC -Wall #-} {-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE DeriveTraversable #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE PatternSynonyms #-} @@ -22,13 +23,40 @@ module Data.DecisionDiagram.BDD.Internal.Node ( -- * Low level node type - Node (T, F, Branch) + Node (Leaf, Branch) , nodeId + + , numNodes + + -- * Fold + , fold + , fold' + , mkFold'Op + + -- * (Co)algebraic structure + , Sig (..) + + -- * Graph + , Graph + , toGraph + , toGraph' + , foldGraph + , foldGraphNodes ) where +import Control.Monad +import Control.Monad.ST +import Control.Monad.ST.Unsafe +import Data.Functor.Identity import Data.Hashable +import qualified Data.HashTable.Class as H +import qualified Data.HashTable.ST.Cuckoo as C import Data.Interned +import Data.IntMap.Lazy (IntMap) +import qualified Data.IntMap.Lazy as IntMap +import Data.STRef import GHC.Generics +import GHC.Stack -- ------------------------------------------------------------------------ @@ -50,11 +78,22 @@ pattern F <- (unintern -> UF) where F = intern UF +pattern Leaf :: Bool -> Node +pattern Leaf b <- (asBool -> Just b) where + Leaf True = T + Leaf False = F + +asBool :: Node -> Maybe Bool +asBool T = Just True +asBool F = Just False +asBool _ = Nothing + pattern Branch :: Int -> Node -> Node -> Node pattern Branch ind lo hi <- (unintern -> UBranch ind lo hi) where Branch ind lo hi = intern (UBranch ind lo hi) {-# COMPLETE T, F, Branch #-} +{-# COMPLETE Leaf, Branch #-} data UNode = UT @@ -86,5 +125,143 @@ nodeId :: Node -> Id nodeId (Node id_ _) = id_ + +-- ------------------------------------------------------------------------ + +defaultTableSize :: Int +defaultTableSize = 256 + +-- | Counts the number of nodes when viewed as a rooted directed acyclic graph +numNodes :: Node -> Int +numNodes node0 = runST $ do + h <- C.newSized defaultTableSize + let f node = do + m <- H.lookup h node + case m of + Just _ -> return () + Nothing -> do + H.insert h node () + case node of + Branch _ lo hi -> f lo >> f hi + _ -> return () + f node0 + liftM length $ H.toList h + +-- ------------------------------------------------------------------------ + +-- | Signature functor of binary decision trees, BDD, and ZDD. +data Sig a + = SLeaf !Bool + | SBranch !Int a a + deriving (Eq, Ord, Show, Read, Generic, Functor, Foldable, Traversable) + +instance Hashable a => Hashable (Sig a) + +-- ------------------------------------------------------------------------ + +-- | Fold over the graph structure of Node. +-- +-- It takes two functions that substitute 'Branch' and 'Leaf' respectively. +-- +-- Note that its type is isomorphic to @('Sig' a -> a) -> 'Node' -> a@. +fold :: (Int -> a -> a -> a) -> (Bool -> a) -> Node -> a +fold br lf bdd = runST $ do + h <- C.newSized defaultTableSize + let f (Leaf b) = return (lf b) + f p@(Branch top lo hi) = do + m <- H.lookup h p + case m of + Just ret -> return ret + Nothing -> do + r0 <- unsafeInterleaveST $ f lo + r1 <- unsafeInterleaveST $ f hi + let ret = br top r0 r1 + H.insert h p ret -- Note that H.insert is value-strict + return ret + f bdd + +-- | Strict version of 'fold' +fold' :: (Int -> a -> a -> a) -> (Bool -> a) -> Node -> a +fold' br lf bdd = runST $ do + op <- mkFold'Op br lf + op bdd + +mkFold'Op :: (Int -> a -> a -> a) -> (Bool -> a) -> ST s (Node -> ST s a) +mkFold'Op br lf = do + h <- C.newSized defaultTableSize + let f (Leaf b) = return $! lf b + f p@(Branch top lo hi) = do + m <- H.lookup h p + case m of + Just ret -> return ret + Nothing -> do + r0 <- f lo + r1 <- f hi + let ret = br top r0 r1 + H.insert h p ret -- Note that H.insert is value-strict + return ret + return f + +-- ------------------------------------------------------------------------ + +-- | Graph where nodes are decorated using a functor @f@. +-- +-- The occurrences of the parameter of @f@ represent out-going edges. +type Graph f = IntMap (f Int) + +-- | Convert a node into a pointed graph +-- +-- Nodes @0@ and @1@ are reserved for @SLeaf False@ and @SLeaf True@ even if +-- they are not actually used. Therefore the result may be larger than +-- 'numNodes' if the leaf nodes are not used. +toGraph :: Node -> (Graph Sig, Int) +toGraph bdd = + case toGraph' (Identity bdd) of + (g, Identity v) -> (g, v) + +-- | Convert multiple nodes into a graph +toGraph' :: Traversable t => t Node -> (Graph Sig, t Int) +toGraph' bs = runST $ do + h <- C.newSized defaultTableSize + H.insert h (Leaf False) 0 + H.insert h (Leaf True) 1 + counter <- newSTRef 2 + ref <- newSTRef $ IntMap.fromList [(0, SLeaf False), (1, SLeaf True)] + + let f (Leaf False) = return 0 + f (Leaf True) = return 1 + f p@(Branch x lo hi) = do + m <- H.lookup h p + case m of + Just ret -> return ret + Nothing -> do + r0 <- f lo + r1 <- f hi + n <- readSTRef counter + writeSTRef counter $! n+1 + H.insert h p n + modifySTRef' ref (IntMap.insert n (SBranch x r0 r1)) + return n + + vs <- mapM f bs + g <- readSTRef ref + return (g, vs) + +-- | Fold over pointed graph +foldGraph :: (Functor f, HasCallStack) => (f a -> a) -> (Graph f, Int) -> a +foldGraph f (g, v) = + case IntMap.lookup v (foldGraphNodes f g) of + Just x -> x + Nothing -> error ("foldGraphNodes: invalid node id " ++ show v) + +-- | Fold over graph nodes +foldGraphNodes :: (Functor f, HasCallStack) => (f a -> a) -> Graph f -> IntMap a +foldGraphNodes f m = ret + where + ret = IntMap.map (f . fmap h) m + h v = + case IntMap.lookup v ret of + Just x -> x + Nothing -> error ("foldGraphNodes: invalid node id " ++ show v) -- ------------------------------------------------------------------------ diff --git a/src/Data/DecisionDiagram/ZDD.hs b/src/Data/DecisionDiagram/ZDD.hs --- a/src/Data/DecisionDiagram/ZDD.hs +++ b/src/Data/DecisionDiagram/ZDD.hs @@ -30,7 +30,9 @@ module Data.DecisionDiagram.ZDD ( -- * ZDD type - ZDD (Empty, Base, Branch) + ZDD (Leaf, Branch) + , pattern Empty + , pattern Base -- * Item ordering , ItemOrder (..) @@ -46,9 +48,16 @@ , base , singleton , subsets + , combinations , fromListOfIntSets , fromSetOfIntSets + -- ** Pseudo-boolean constraints + , subsetsAtLeast + , subsetsAtMost + , subsetsExactly + , subsetsExactlyIntegral + -- * Insertion , insert @@ -63,6 +72,7 @@ , isSubsetOf , isProperSubsetOf , disjoint + , numNodes -- * Combine , union @@ -81,10 +91,21 @@ , mapDelete , change + -- * (Co)algebraic structure + , Sig (..) + , pattern SEmpty + , pattern SBase + , inSig + , outSig + -- * Fold , fold , fold' + -- * Unfold + , unfoldHashable + , unfoldOrd + -- * Minimal hitting sets , minimalHittingSets , minimalHittingSetsToda @@ -107,7 +128,6 @@ -- ** Conversion from/to graphs , Graph - , Node (..) , toGraph , toGraph' , fromGraph @@ -121,7 +141,8 @@ import Control.Monad.Primitive #endif import Control.Monad.ST -import Data.Functor.Identity +import qualified Data.Foldable as Foldable +import Data.Function (on) import Data.Hashable import Data.HashMap.Lazy (HashMap) import qualified Data.HashMap.Lazy as HashMap @@ -132,13 +153,16 @@ import Data.IntSet (IntSet) import qualified Data.IntSet as IntSet import Data.List (foldl', sortBy) +import Data.Map.Lazy (Map) +import qualified Data.Map.Lazy as Map import Data.Maybe import Data.Proxy import Data.Ratio import Data.Set (Set) import qualified Data.Set as Set -import Data.STRef +import qualified Data.Vector as V import qualified GHC.Exts as Exts +import GHC.Stack import Numeric.Natural #if MIN_VERSION_mwc_random(0,15,0) import System.Random.Stateful (StatefulGen (..)) @@ -149,6 +173,7 @@ import Text.Read import Data.DecisionDiagram.BDD.Internal.ItemOrder +import Data.DecisionDiagram.BDD.Internal.Node (Sig (..), Graph) import qualified Data.DecisionDiagram.BDD.Internal.Node as Node import qualified Data.DecisionDiagram.BDD as BDD @@ -163,12 +188,17 @@ newtype ZDD a = ZDD Node.Node deriving (Eq, Hashable) +-- | Synonym of @'Leaf' False@ pattern Empty :: ZDD a -pattern Empty = ZDD Node.F +pattern Empty = Leaf False +-- | Synonym of @'Leaf' True@ pattern Base :: ZDD a -pattern Base = ZDD Node.T +pattern Base = Leaf True +pattern Leaf :: Bool -> ZDD a +pattern Leaf b = ZDD (Node.Leaf b) + -- | Smart constructor that takes the ZDD reduction rules into account pattern Branch :: Int -> ZDD a -> ZDD a -> ZDD a pattern Branch x lo hi <- ZDD (Node.Branch x (ZDD -> lo) (ZDD -> hi)) where @@ -176,7 +206,13 @@ Branch x (ZDD lo) (ZDD hi) = ZDD (Node.Branch x lo hi) {-# COMPLETE Empty, Base, Branch #-} +{-# COMPLETE Leaf, Branch #-} +-- Hack for avoiding spurious incomplete patterns warning on the above Branch pattern definition. +#if __GLASGOW_HASKELL__ < 810 +{-# COMPLETE ZDD #-} +#endif + nodeId :: ZDD a -> Int nodeId (ZDD node) = Node.nodeId node @@ -202,7 +238,7 @@ f :: IntSet -> [Int] f = sortBy (compareItem (Proxy :: Proxy a)) . IntSet.toList - toList = fold' [] [IntSet.empty] (\top lo hi -> lo <> map (IntSet.insert top) hi) + toList = toListOfIntSets -- ------------------------------------------------------------------------ @@ -225,15 +261,24 @@ zddCase2 _ Empty Base = ZDDCase2EQ2 False True zddCase2 _ Empty Empty = ZDDCase2EQ2 False False --- | The empty set (∅). +-- | The empty family (∅). +-- +-- >>> toSetOfIntSets (empty :: ZDD AscOrder) +-- fromList [] empty :: ZDD a empty = Empty --- | The set containing only the empty set ({∅}). +-- | The family containing only the empty set ({∅}). +-- +-- >>> toSetOfIntSets (base :: ZDD AscOrder) +-- fromList [fromList []] base :: ZDD a base = Base -- | Create a ZDD that contains only a given set. +-- +-- >>> toSetOfIntSets (singleton (IntSet.fromList [1,2,3]) :: ZDD AscOrder) +-- fromList [fromList [1,2,3]] singleton :: forall a. ItemOrder a => IntSet -> ZDD a singleton xs = insert xs empty @@ -243,7 +288,107 @@ where f zdd x = Branch x zdd zdd +-- | Set of all k-combination of a set +combinations :: forall a. (ItemOrder a, HasCallStack) => IntSet -> Int -> ZDD a +combinations xs k + | k < 0 = error "Data.DecisionDiagram.ZDD.combinations: negative size" + | otherwise = unfoldOrd f (0, k) + where + table = V.fromList $ sortBy (compareItem (Proxy :: Proxy a)) $ IntSet.toList xs + n = V.length table + + f :: (Int, Int) -> Sig (Int, Int) + f (!_, !0) = SLeaf True + f (!i, !k') + | i + k' > n = SLeaf False + | otherwise = SBranch (table V.! i) (i+1, k') (i+1, k'-1) + +-- | Set of all subsets whose sum of weights is at least k. +subsetsAtLeast :: forall a w. (ItemOrder a, Real w) => IntMap w -> w -> ZDD a +subsetsAtLeast xs k0 = unfoldOrd f (0, k0) + where + xs' :: V.Vector (Int, w) + xs' = V.fromList $ sortBy (compareItem (Proxy :: Proxy a) `on` fst) $ IntMap.toList xs + ys :: V.Vector (w, w) + ys = V.scanr (\(_, w) (lb,ub) -> if w >= 0 then (lb, ub+w) else (lb+w, ub)) (0,0) xs' + + f :: (Int, w) -> Sig (Int, w) + f (!i, !k) + | not (k <= ub) = SLeaf False + | i == V.length xs' && 0 >= k = SLeaf True + | lb >= k = SBranch x (i+1, lb) (i+1, lb) -- all remaining variables are don't-care + | otherwise = SBranch x (i+1, k) (i+1, k-w) + where + (lb,ub) = ys V.! i + (x, w) = xs' V.! i + +-- | Set of all subsets whose sum of weights is at most k. +subsetsAtMost :: forall a w. (ItemOrder a, Real w) => IntMap w -> w -> ZDD a +subsetsAtMost xs k0 = unfoldOrd f (0, k0) + where + xs' :: V.Vector (Int, w) + xs' = V.fromList $ sortBy (compareItem (Proxy :: Proxy a) `on` fst) $ IntMap.toList xs + ys :: V.Vector (w, w) + ys = V.scanr (\(_, w) (lb,ub) -> if w >= 0 then (lb, ub+w) else (lb+w, ub)) (0,0) xs' + + f :: (Int, w) -> Sig (Int, w) + f (!i, !k) + | not (lb <= k) = SLeaf False + | i == V.length xs' && 0 <= k = SLeaf True + | ub <= k = SBranch x (i+1, ub) (i+1, ub) -- all remaining variables are don't-care + | otherwise = SBranch x (i+1, k) (i+1, k-w) + where + (lb,ub) = ys V.! i + (x, w) = xs' V.! i + +-- | Set of all subsets whose sum of weights is exactly k. +-- +-- Note that 'combinations' is a special case where all weights are 1. +-- +-- If weight type is 'Integral', 'subsetsExactlyIntegral' is more efficient. +subsetsExactly :: forall a w. (ItemOrder a, Real w) => IntMap w -> w -> ZDD a +subsetsExactly xs k0 = unfoldOrd f (0, k0) + where + xs' :: V.Vector (Int, w) + xs' = V.fromList $ sortBy (compareItem (Proxy :: Proxy a) `on` fst) $ IntMap.toList xs + ys :: V.Vector (w, w) + ys = V.scanr (\(_, w) (lb,ub) -> if w >= 0 then (lb, ub+w) else (lb+w, ub)) (0,0) xs' + + f :: (Int, w) -> Sig (Int, w) + f (!i, !k) + | not (lb <= k && k <= ub) = SLeaf False + | i == V.length xs' && 0 == k = SLeaf True + | otherwise = SBranch x (i+1, k) (i+1, k-w) + where + (lb,ub) = ys V.! i + (x, w) = xs' V.! i + +-- | Similar to 'subsetsExactly' but more efficient. +subsetsExactlyIntegral :: forall a w. (ItemOrder a, Real w, Integral w) => IntMap w -> w -> ZDD a +subsetsExactlyIntegral xs k0 = unfoldOrd f (0, k0) + where + xs' :: V.Vector (Int, w) + xs' = V.fromList $ sortBy (compareItem (Proxy :: Proxy a) `on` fst) $ IntMap.toList xs + ys :: V.Vector (w, w) + ys = V.scanr (\(_, w) (lb,ub) -> if w >= 0 then (lb, ub+w) else (lb+w, ub)) (0,0) xs' + ds :: V.Vector w + ds = V.scanr1 (\w d -> if w /= 0 then gcd w d else d) (V.map snd xs') + + f :: (Int, w) -> Sig (Int, w) + f (!i, !k) + | not (lb <= k && k <= ub) = SLeaf False + | i == V.length xs' && 0 == k = SLeaf True + | d /= 0 && k `mod` d /= 0 = SLeaf False + | otherwise = SBranch x (i+1, k) (i+1, k-w) + where + (lb,ub) = ys V.! i + (x, w) = xs' V.! i + d = ds V.! i + -- | Select subsets that contain a particular element and then remove the element from them +-- +-- >>> toSetOfIntSets $ subset1 2 (fromListOfIntSets (map IntSet.fromList [[1,2,3], [1,3], [2,4]]) :: ZDD AscOrder) +-- fromList [fromList [1,3],fromList [4]] subset1 :: forall a. ItemOrder a => Int -> ZDD a -> ZDD a subset1 var zdd = runST $ do h <- C.newSized defaultTableSize @@ -263,6 +408,9 @@ f zdd -- | Subsets that does not contain a particular element +-- +-- >>> toSetOfIntSets $ subset0 2 (fromListOfIntSets (map IntSet.fromList [[1,2,3], [1,3], [2,4], [3,4]]) :: ZDD AscOrder) +-- fromList [fromList [1,3],fromList [3,4]] subset0 :: forall a. ItemOrder a => Int -> ZDD a -> ZDD a subset0 var zdd = runST $ do h <- C.newSized defaultTableSize @@ -282,11 +430,13 @@ f zdd -- | Insert a set into the ZDD. +-- +-- >>> toSetOfIntSets (insert (IntSet.fromList [1,2,3]) (fromListOfIntSets (map IntSet.fromList [[1,3], [2,4]])) :: ZDD AscOrder) +-- fromList [fromList [1,2,3],fromList [1,3],fromList [2,4]] insert :: forall a. ItemOrder a => IntSet -> ZDD a -> ZDD a insert xs = f (sortBy (compareItem (Proxy :: Proxy a)) (IntSet.toList xs)) where - f [] Empty = Base - f [] Base = Base + f [] (Leaf _) = Base f [] (Branch top p0 p1) = Branch top (f [] p0) p1 f (y : ys) Empty = Branch y Empty (f ys Empty) f (y : ys) Base = Branch y Base (f ys Empty) @@ -297,14 +447,15 @@ EQ -> Branch top p0 (f ys p1) -- | Delete a set from the ZDD. +-- +-- >>> toSetOfIntSets (delete (IntSet.fromList [1,3]) (fromListOfIntSets (map IntSet.fromList [[1,2,3], [1,3], [2,4]])) :: ZDD AscOrder) +-- fromList [fromList [1,2,3],fromList [2,4]] delete :: forall a. ItemOrder a => IntSet -> ZDD a -> ZDD a delete xs = f (sortBy (compareItem (Proxy :: Proxy a)) (IntSet.toList xs)) where - f [] Empty = Empty - f [] Base = Empty + f [] (Leaf _) = Empty f [] (Branch top p0 p1) = Branch top (f [] p0) p1 - f (_ : _) Empty = Empty - f (_ : _) Base = Base + f (_ : _) l@(Leaf _) = l f yys@(y : ys) p@(Branch top p0 p1) = case compareItem (Proxy :: Proxy a) y top of LT -> p @@ -312,6 +463,9 @@ EQ -> Branch top p0 (f ys p1) -- | Insert an item into each element set of ZDD. +-- +-- >>> toSetOfIntSets (mapInsert 2 (fromListOfIntSets (map IntSet.fromList [[1,2,3], [1,3], [1,4]])) :: ZDD AscOrder) +-- fromList [fromList [1,2,3],fromList [1,2,4]] mapInsert :: forall a. ItemOrder a => Int -> ZDD a -> ZDD a mapInsert var zdd = runST $ do unionOp <- mkUnionOp @@ -332,12 +486,14 @@ f zdd -- | Delete an item from each element set of ZDD. +-- +-- >>> toSetOfIntSets (mapDelete 2 (fromListOfIntSets (map IntSet.fromList [[1,2,3], [1,3], [1,2,4]])) :: ZDD AscOrder) +-- fromList [fromList [1,3],fromList [1,4]] mapDelete :: forall a. ItemOrder a => Int -> ZDD a -> ZDD a mapDelete var zdd = runST $ do unionOp <- mkUnionOp h <- C.newSized defaultTableSize - let f Base = return Base - f Empty = return Empty + let f l@(Leaf _) = return l f p@(Branch top p0 p1) = do m <- H.lookup h p case m of @@ -352,6 +508,9 @@ f zdd -- | @change x p@ returns {if x∈s then s∖{x} else s∪{x} | s∈P} +-- +-- >>> toSetOfIntSets (change 2 (fromListOfIntSets (map IntSet.fromList [[1,2,3], [1,3], [1,2,4]])) :: ZDD AscOrder) +-- fromList [fromList [1,2,3],fromList [1,3],fromList [1,4]] change :: forall a. ItemOrder a => Int -> ZDD a -> ZDD a change var zdd = runST $ do h <- C.newSized defaultTableSize @@ -463,14 +622,17 @@ -- | Given a family P and Q, it computes {S∈P | ∀X∈Q. X⊈S} -- -- Sometimes it is denoted as /P ↘ Q/. +-- +-- >>> toSetOfIntSets (fromListOfIntSets (map IntSet.fromList [[1,2,3], [1,3], [3,4]]) `nonSuperset` singleton (IntSet.fromList [1,3]) :: ZDD AscOrder) +-- fromList [fromList [3,4]] nonSuperset :: forall a. ItemOrder a => ZDD a -> ZDD a -> ZDD a nonSuperset zdd1 zdd2 = runST $ do - op <- mkNonSueprsetOp + op <- mkNonSupersetOp op zdd1 zdd2 -mkNonSueprsetOp :: forall a s. ItemOrder a => ST s (ZDD a -> ZDD a -> ST s (ZDD a)) -mkNonSueprsetOp = do - intersectionOp <- mkIntersectionOp +mkNonSupersetOp :: forall a s. ItemOrder a => ST s (ZDD a -> ZDD a -> ST s (ZDD a)) +mkNonSupersetOp = do + intersectionOp <- mkIntersectionOp h <- C.newSized defaultTableSize let f Empty _ = return Empty f _ Base = return Empty @@ -496,7 +658,7 @@ minimalHittingSetsKnuth' :: forall a. ItemOrder a => Bool -> ZDD a -> ZDD a minimalHittingSetsKnuth' imai zdd = runST $ do unionOp <- mkUnionOp - diffOp <- if imai then mkDifferenceOp else mkNonSueprsetOp + diffOp <- if imai then mkDifferenceOp else mkNonSupersetOp h <- C.newSized defaultTableSize let f Empty = return Base f Base = return Empty @@ -540,14 +702,13 @@ minimalHittingSetsToda = minimal . hittingSetsBDD hittingSetsBDD :: forall a. ItemOrder a => ZDD a -> BDD.BDD a -hittingSetsBDD = fold' BDD.true BDD.false (\top h0 h1 -> h0 BDD..&&. BDD.Branch top h1 BDD.true) +hittingSetsBDD = fold' (\top h0 h1 -> h0 BDD..&&. BDD.Branch top h1 BDD.true) (\b -> BDD.Leaf (not b)) minimal :: forall a. ItemOrder a => BDD.BDD a -> ZDD a minimal bdd = runST $ do diffOp <- mkDifferenceOp h <- C.newSized defaultTableSize - let f BDD.F = return Empty - f BDD.T = return Base + let f (BDD.Leaf b) = return (Leaf b) f p@(BDD.Branch x lo hi) = do m <- H.lookup h p case m of @@ -561,6 +722,9 @@ f bdd -- | See 'minimalHittingSetsToda'. +-- +-- >>> toSetOfIntSets (minimalHittingSets (fromListOfIntSets (map IntSet.fromList [[1], [2,3,5], [2,3,6], [2,4,5], [2,4,6]]) :: ZDD AscOrder)) +-- fromList [fromList [1,2],fromList [1,3,4],fromList [1,5,6]] minimalHittingSets :: forall a. ItemOrder a => ZDD a -> ZDD a minimalHittingSets = minimalHittingSetsToda @@ -584,7 +748,7 @@ notMember :: forall a. (ItemOrder a) => IntSet -> ZDD a -> Bool notMember xs = not . member xs --- | Is this the empty set? +-- | Is this the empty family? null :: ZDD a -> Bool null = (empty ==) @@ -592,24 +756,45 @@ {-# SPECIALIZE size :: ZDD a -> Integer #-} {-# SPECIALIZE size :: ZDD a -> Natural #-} -- | The number of sets in the family. +-- +-- Any 'Integral' type can be used as a result type, but it is recommended to use +-- 'Integer' or 'Natural' because the size can be larger than @Int64@ for example: +-- +-- >>> size (subsets (IntSet.fromList [1..128]) :: ZDD AscOrder) :: Integer +-- 340282366920938463463374607431768211456 +-- >>> import Data.Int +-- >>> maxBound :: Int64 +-- 9223372036854775807 +-- size :: (Integral b) => ZDD a -> b -size = fold' 0 1 (\_ n0 n1 -> n0 + n1) +size = fold' (\_ n0 n1 -> n0 + n1) (\b -> if b then 1 else 0) --- | @(s1 `isSubsetOf` s2)@ indicates whether @s1@ is a subset of @s2@. +-- | @(s1 \`isSubsetOf\` s2)@ indicates whether @s1@ is a subset of @s2@. isSubsetOf :: ItemOrder a => ZDD a -> ZDD a -> Bool isSubsetOf a b = union a b == b --- | @(s1 `isProperSubsetOf` s2)@ indicates whether @s1@ is a proper subset of @s2@. +-- | @(s1 \`isProperSubsetOf\` s2)@ indicates whether @s1@ is a proper subset of @s2@. isProperSubsetOf :: ItemOrder a => ZDD a -> ZDD a -> Bool isProperSubsetOf a b = a `isSubsetOf` b && a /= b --- | Check whether two sets are disjoint (i.e., their intersection is empty). +-- | Check whether two families are disjoint (i.e., their intersection is empty). disjoint :: ItemOrder a => ZDD a -> ZDD a -> Bool disjoint a b = null (a `intersection` b) ---- | Unions of all member sets +-- | Count the number of nodes in a ZDD viewed as a rooted directed acyclic graph. +-- +-- Please do not confuse it with 'size'. +-- +-- See also 'toGraph'. +numNodes :: ZDD a -> Int +numNodes (ZDD node) = Node.numNodes node + +-- | Unions of all member sets +-- +-- >>> flatten (fromListOfIntSets (map IntSet.fromList [[1,2,3], [1,3], [3,4]]) :: ZDD AscOrder) +-- fromList [1,2,3,4] flatten :: ItemOrder a => ZDD a -> IntSet -flatten = fold' IntSet.empty IntSet.empty (\top lo hi -> IntSet.insert top (lo `IntSet.union` hi)) +flatten = fold' (\top lo hi -> IntSet.insert top (lo `IntSet.union` hi)) (const IntSet.empty) -- | Create a ZDD from a set of 'IntSet' fromSetOfIntSets :: forall a. ItemOrder a => Set IntSet -> ZDD a @@ -617,7 +802,7 @@ -- | Convert the family to a set of 'IntSet'. toSetOfIntSets :: ZDD a -> Set IntSet -toSetOfIntSets = fold' Set.empty (Set.singleton IntSet.empty) (\top lo hi -> lo <> Set.map (IntSet.insert top) hi) +toSetOfIntSets = fold' (\top lo hi -> lo <> Set.map (IntSet.insert top) hi) (\b -> if b then Set.singleton IntSet.empty else Set.empty) -- | Create a ZDD from a list of 'IntSet' fromListOfIntSets :: forall a. ItemOrder a => [IntSet] -> ZDD a @@ -628,7 +813,11 @@ -- | Convert the family to a list of 'IntSet'. toListOfIntSets :: ZDD a -> [IntSet] -toListOfIntSets = fold [] [IntSet.empty] (\top lo hi -> lo <> map (IntSet.insert top) hi) +toListOfIntSets = g . fold' f (\b -> (b,[])) + where + f top (b, xss) hi = (b, map (IntSet.insert top) (g hi) <> xss) + g (True, xss) = IntSet.empty : xss + g (False, xss) = xss fromListOfSortedList :: forall a. ItemOrder a => [[Int]] -> ZDD a fromListOfSortedList = unions . map f @@ -638,43 +827,54 @@ -- | Fold over the graph structure of the ZDD. -- --- It takes values for substituting 'empty' and 'base', --- and a function for substiting non-terminal node. -fold :: b -> b -> (Int -> b -> b -> b) -> ZDD a -> b -fold ff tt br zdd = runST $ do - h <- C.newSized defaultTableSize - let f Empty = return ff - f Base = return tt - f p@(Branch top p0 p1) = do - m <- H.lookup h p - case m of - Just ret -> return ret - Nothing -> do - r0 <- f p0 - r1 <- f p1 - let ret = br top r0 r1 - H.insert h p ret - return ret - f zdd +-- It takes two functions that substitute 'Branch' and 'Leaf' respectively. +-- +-- Note that its type is isomorphic to @('Sig' b -> b) -> ZDD a -> b@. +fold :: (Int -> b -> b -> b) -> (Bool -> b) -> ZDD a -> b +fold br lf (ZDD node) = Node.fold br lf node -- | Strict version of 'fold' -fold' :: b -> b -> (Int -> b -> b -> b) -> ZDD a -> b -fold' !ff !tt br zdd = runST $ do +fold' :: (Int -> b -> b -> b) -> (Bool -> b) -> ZDD a -> b +fold' br lf (ZDD node) = Node.fold' br lf node + +-- ------------------------------------------------------------------------ + +-- | Top-down construction of ZDD, memoising internal states using 'Hashable' instance. +unfoldHashable :: forall a b. (ItemOrder a, Eq b, Hashable b) => (b -> Sig b) -> b -> ZDD a +unfoldHashable f b = runST $ do h <- C.newSized defaultTableSize - let f Empty = return ff - f Base = return tt - f p@(Branch top p0 p1) = do - m <- H.lookup h p - case m of - Just ret -> return ret + let g [] = return () + g (x : xs) = do + r <- H.lookup h x + case r of + Just _ -> g xs Nothing -> do - r0 <- f p0 - r1 <- f p1 - let ret = br top r0 r1 - seq ret $ H.insert h p ret - return ret - f zdd + let fx = f x + H.insert h x fx + g (xs ++ Foldable.toList fx) + g [b] + xs <- H.toList h + let h2 = HashMap.fromList [(x, inSig (fmap (h2 HashMap.!) s)) | (x,s) <- xs] + return $ h2 HashMap.! b +-- | Top-down construction of ZDD, memoising internal states using 'Ord' instance. +unfoldOrd :: forall a b. (ItemOrder a, Ord b) => (b -> Sig b) -> b -> ZDD a +unfoldOrd f b = m2 Map.! b + where + m1 :: Map b (Sig b) + m1 = g Map.empty [b] + + m2 :: Map b (ZDD a) + m2 = Map.map (inSig . fmap (m2 Map.!)) m1 + + g m [] = m + g m (x : xs) = + case Map.lookup x m of + Just _ -> g m xs + Nothing -> + let fx = f x + in g (Map.insert x fx m) (xs ++ Foldable.toList fx) + -- ------------------------------------------------------------------------ -- | Sample a set from uniform distribution over elements of the ZDD. @@ -697,9 +897,9 @@ -- @ -- . #if MIN_VERSION_mwc_random(0,15,0) -uniformM :: forall a g m. (ItemOrder a, StatefulGen g m) => ZDD a -> g -> m IntSet +uniformM :: forall a g m. (ItemOrder a, StatefulGen g m, HasCallStack) => ZDD a -> g -> m IntSet #else -uniformM :: forall a m. (ItemOrder a, PrimMonad m) => ZDD a -> Gen (PrimState m) -> m IntSet +uniformM :: forall a m. (ItemOrder a, PrimMonad m, HasCallStack) => ZDD a -> Gen (PrimState m) -> m IntSet #endif uniformM Empty = error "Data.DecisionDiagram.ZDD.uniformM: empty ZDD" uniformM zdd = func @@ -743,10 +943,10 @@ -- \[ -- \min_{X\in S} \sum_{x\in X} w(x) -- \] -findMinSum :: forall a w. (ItemOrder a, Num w, Ord w) => (Int -> w) -> ZDD a -> (w, IntSet) +findMinSum :: forall a w. (ItemOrder a, Num w, Ord w, HasCallStack) => (Int -> w) -> ZDD a -> (w, IntSet) findMinSum weight = fromMaybe (error "Data.DecisionDiagram.ZDD.findMinSum: empty ZDD") . - fold' Nothing (Just (0, IntSet.empty)) f + fold' f (\b -> if b then Just (0, IntSet.empty) else Nothing) where f _ _ Nothing = undefined f x z1 (Just (w2, s2)) = @@ -762,10 +962,13 @@ -- \[ -- \max_{X\in S} \sum_{x\in X} w(x) -- \] -findMaxSum :: forall a w. (ItemOrder a, Num w, Ord w) => (Int -> w) -> ZDD a -> (w, IntSet) +-- +-- >>> findMaxSum (IntMap.fromList [(1,2),(2,4),(3,-3)] IntMap.!) (fromListOfIntSets (map IntSet.fromList [[1], [2], [3], [1,2,3]]) :: ZDD AscOrder) +-- (4,fromList [2]) +findMaxSum :: forall a w. (ItemOrder a, Num w, Ord w, HasCallStack) => (Int -> w) -> ZDD a -> (w, IntSet) findMaxSum weight = fromMaybe (error "Data.DecisionDiagram.ZDD.findMinSum: empty ZDD") . - fold' Nothing (Just (0, IntSet.empty)) f + fold' f (\b -> if b then Just (0, IntSet.empty) else Nothing) where f _ _ Nothing = undefined f x z1 (Just (w2, s2)) = @@ -778,66 +981,44 @@ -- ------------------------------------------------------------------------ -type Graph = IntMap Node +-- | Synonym of @'SLeaf' False@ +pattern SEmpty :: Sig a +pattern SEmpty = SLeaf False -data Node - = NodeEmpty - | NodeBase - | NodeBranch !Int Int Int - deriving (Eq, Show, Read) +-- | Synonym of @'SLeaf' True@ +pattern SBase :: Sig a +pattern SBase = SLeaf True --- | Convert a ZDD into a pointed graph -toGraph :: ZDD a -> (Graph, Int) -toGraph bdd = - case toGraph' (Identity bdd) of - (g, Identity v) -> (g, v) +-- | 'Sig'-algebra stucture of 'ZDD', \(\mathrm{in}_\mathrm{Sig}\). +inSig :: Sig (ZDD a) -> ZDD a +inSig (SLeaf b) = Leaf b +inSig (SBranch x lo hi) = Branch x lo hi --- | Convert multiple ZDDs into a graph -toGraph' :: Traversable t => t (ZDD a) -> (Graph, t Int) -toGraph' bs = runST $ do - h <- C.newSized defaultTableSize - H.insert h Empty 0 - H.insert h Base 1 - counter <- newSTRef 2 - ref <- newSTRef $ IntMap.fromList [(0, NodeEmpty), (1, NodeBase)] +-- | 'Sig'-coalgebra stucture of 'ZDD', \(\mathrm{out}_\mathrm{Sig}\). +outSig :: ZDD a -> Sig (ZDD a) +outSig (Leaf b) = SLeaf b +outSig (Branch x lo hi) = SBranch x lo hi - let f Empty = return 0 - f Base = return 1 - f p@(Branch x lo hi) = do - m <- H.lookup h p - case m of - Just ret -> return ret - Nothing -> do - r0 <- f lo - r1 <- f hi - n <- readSTRef counter - writeSTRef counter $! n+1 - H.insert h p n - modifySTRef' ref (IntMap.insert n (NodeBranch x r0 r1)) - return n +-- ------------------------------------------------------------------------ - vs <- mapM f bs - g <- readSTRef ref - return (g, vs) +-- | Convert a ZDD into a pointed graph +-- +-- Nodes @0@ and @1@ are reserved for @SLeaf False@ and @SLeaf True@ even if +-- they are not actually used. Therefore the result may be larger than +-- 'numNodes' if the leaf nodes are not used. +toGraph :: ZDD a -> (Graph Sig, Int) +toGraph (ZDD node) = Node.toGraph node +-- | Convert multiple ZDDs into a graph +toGraph' :: Traversable t => t (ZDD a) -> (Graph Sig, t Int) +toGraph' bs = Node.toGraph' (fmap (\(ZDD node) -> node) bs) + -- | Convert a pointed graph into a ZDD -fromGraph :: (Graph, Int) -> ZDD a -fromGraph (g, v) = - case IntMap.lookup v (fromGraph' g) of - Nothing -> error ("Data.DecisionDiagram.ZDD.fromGraph: invalid node id " ++ show v) - Just bdd -> bdd +fromGraph :: HasCallStack => (Graph Sig, Int) -> ZDD a +fromGraph = Node.foldGraph inSig -- | Convert nodes of a graph into ZDDs -fromGraph' :: Graph -> IntMap (ZDD a) -fromGraph' g = ret - where - ret = IntMap.map f g - f NodeEmpty = Empty - f NodeBase = Base - f (NodeBranch x lo hi) = - case (IntMap.lookup lo ret, IntMap.lookup hi ret) of - (Nothing, _) -> error ("Data.DecisionDiagram.ZDD.fromGraph': invalid node id " ++ show lo) - (_, Nothing) -> error ("Data.DecisionDiagram.ZDD.fromGraph': invalid node id " ++ show hi) - (Just lo', Just hi') -> Branch x lo' hi' +fromGraph' :: HasCallStack => Graph Sig -> IntMap (ZDD a) +fromGraph' = Node.foldGraphNodes inSig -- ------------------------------------------------------------------------ diff --git a/test/TestBDD.hs b/test/TestBDD.hs --- a/test/TestBDD.hs +++ b/test/TestBDD.hs @@ -4,12 +4,24 @@ module TestBDD (bddTestGroup) where import Control.Monad -import qualified Data.IntMap as IntMap +import Control.Monad.ST +import Data.IntMap.Lazy (IntMap) +import qualified Data.IntMap.Lazy as IntMap import Data.IntSet (IntSet) import qualified Data.IntSet as IntSet +import Data.IORef import Data.List +import qualified Data.Map.Lazy as Map import Data.Proxy +import qualified Data.Set as Set +import Data.Vector (Vector) +import Data.Word +import Statistics.Distribution +import Statistics.Distribution.ChiSquared (chiSquared) +import System.IO.Unsafe +import qualified System.Random.MWC as Rand import Test.QuickCheck.Function (apply) +import Test.QuickCheck.Instances.Vector () import Test.Tasty import Test.Tasty.HUnit import Test.Tasty.QuickCheck @@ -23,10 +35,9 @@ -- ------------------------------------------------------------------------ instance BDD.ItemOrder a => Arbitrary (BDD a) where - arbitrary = arbitraryBDDOver =<< liftM IntSet.fromList arbitrary + arbitrary = arbitraryBDDOver =<< arbitrary - shrink (BDD.F) = [] - shrink (BDD.T) = [] + shrink (BDD.Leaf _) = [] shrink (BDD.Branch x p0 p1) = [p0, p1] ++ @@ -50,6 +61,25 @@ ] sized $ f (sortBy (BDD.compareItem (Proxy :: Proxy a)) $ IntSet.toList xs) +arbitrarySatisfyingAssignment :: forall a. BDD.ItemOrder a => BDD a -> IntSet -> Gen (IntMap Bool) +arbitrarySatisfyingAssignment bdd xs = do + m1 <- arbitrarySatisfyingPartialAssignment bdd + let ys = xs `IntSet.difference` IntMap.keysSet m1 + m2 <- liftM(IntMap.fromAscList) $ forM (IntSet.toAscList ys) $ \y -> do + v <- arbitrary + return (y,v) + return $ m1 `IntMap.union` m2 + +arbitrarySatisfyingPartialAssignment :: forall a. BDD.ItemOrder a => BDD a -> Gen (IntMap Bool) +arbitrarySatisfyingPartialAssignment = f + where + f (BDD.Leaf True) = return IntMap.empty + f (BDD.Leaf False) = undefined + f (BDD.Branch x lo hi) = oneof $ + [liftM (IntMap.insert x False) (f lo) | lo /= BDD.Leaf False] + ++ + [liftM (IntMap.insert x True) (f hi) | hi /= BDD.Leaf False] + -- ------------------------------------------------------------------------ -- conjunction -- ------------------------------------------------------------------------ @@ -305,6 +335,83 @@ (d BDD..||. BDD.ite c t e) === BDD.ite c (d BDD..||. t) (d BDD..||. e) -- ------------------------------------------------------------------------ +-- Pseudo-Boolean +-- ------------------------------------------------------------------------ + +prop_pbAtLeast :: Property +prop_pbAtLeast = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntMap $ \(xs :: IntMap Integer) -> + forAll arbitrary $ \k -> + let a :: BDD o + a = BDD.pbAtLeast xs k + in counterexample (show a) $ + if a == BDD.Leaf False then + property (k > sum [max 0 w | (_,w) <- IntMap.toList xs]) + else + forAll (arbitrarySatisfyingAssignment a (IntMap.keysSet xs)) $ \ys -> + (IntMap.keysSet ys `IntSet.isSubsetOf` IntMap.keysSet xs) + .&&. + sum [xs IntMap.! y | (y,b) <- IntMap.toList ys, b] >= k + +prop_pbAtMost :: Property +prop_pbAtMost = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntMap $ \(xs :: IntMap Integer) -> + forAll arbitrary $ \k -> + let a :: BDD o + a = BDD.pbAtMost xs k + in counterexample (show a) $ + if a == BDD.Leaf False then + property (k < sum [min 0 w | (_,w) <- IntMap.toList xs]) + else + forAll (arbitrarySatisfyingAssignment a (IntMap.keysSet xs)) $ \ys -> + (IntMap.keysSet ys `IntSet.isSubsetOf` IntMap.keysSet xs) + .&&. + sum [xs IntMap.! y | (y,b) <- IntMap.toList ys, b] <= k + +prop_pbExactly :: Property +prop_pbExactly = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntMap $ \(xs :: IntMap Integer) -> + forAll arbitrary $ \k -> + let a :: BDD o + a = BDD.pbExactly xs k + in counterexample (show a) $ + if a == BDD.Leaf False then + property True + else + forAll (arbitrarySatisfyingAssignment a (IntMap.keysSet xs)) $ \ys -> + (IntMap.keysSet ys `IntSet.isSubsetOf` IntMap.keysSet xs) + .&&. + sum [xs IntMap.! y | (y,b) <- IntMap.toList ys, b] === k + +prop_pbExactly_2 :: Property +prop_pbExactly_2 = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntMap $ \(xs :: IntMap Integer) -> + forAll (gen xs) $ \(m, k) -> + let a :: BDD o + a = BDD.pbExactly xs k + in counterexample (show a) $ BDD.evaluate (m IntMap.!) a + where + gen :: IntMap Integer -> Gen (IntMap Bool, Integer) + gen xs = do + ys <- sublistOf (IntMap.toList xs) + let ys' = IntSet.fromList [y | (y,_) <- ys] + return + ( IntMap.mapWithKey (\x _ -> x `IntSet.member` ys') xs + , sum [w | (_,w) <- ys] + ) + +prop_pbExactlyIntegral :: Property +prop_pbExactlyIntegral = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntMap $ \(xs :: IntMap Integer) -> + forAll arbitrary $ \k -> + (BDD.pbExactlyIntegral xs k :: BDD o) === BDD.pbExactly xs k + +-- ------------------------------------------------------------------------ -- Quantification -- ------------------------------------------------------------------------ @@ -416,6 +523,55 @@ -- ------------------------------------------------------------------------ +case_fold_laziness :: Assertion +case_fold_laziness = do + let bdd :: BDD BDD.AscOrder + bdd = BDD.Branch 0 (BDD.Branch 1 (BDD.Leaf False) (BDD.Leaf True)) (BDD.Branch 2 (BDD.Leaf False) (BDD.Leaf True)) + f x lo _hi = + if x == 2 then + error "unused value should not be evaluated" + else + lo + seq (BDD.fold f id bdd) $ return () + +case_fold'_strictness :: Assertion +case_fold'_strictness = do + ref <- newIORef False + let bdd :: BDD BDD.AscOrder + bdd = BDD.Branch 0 (BDD.Branch 1 (BDD.Leaf False) (BDD.Leaf True)) (BDD.Branch 2 (BDD.Leaf False) (BDD.Leaf True)) + f x lo _hi = unsafePerformIO $ do + when (x==2) $ writeIORef ref True + return lo + seq (BDD.fold' f id bdd) $ do + flag <- readIORef ref + assertBool "unused value should be evaluated" flag + +prop_fold_inSig :: Property +prop_fold_inSig = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrary $ \(bdd :: BDD o) -> + BDD.fold (\x lo hi -> BDD.inSig (BDD.SBranch x lo hi)) (BDD.inSig . BDD.SLeaf) bdd + === + bdd + +prop_fold'_inSig :: Property +prop_fold'_inSig = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrary $ \(bdd :: BDD o) -> + BDD.fold' (\x lo hi -> BDD.inSig (BDD.SBranch x lo hi)) (BDD.inSig . BDD.SLeaf) bdd + === + bdd + +-- ------------------------------------------------------------------------ + +prop_unfoldHashable_outSig :: Property +prop_unfoldHashable_outSig = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrary $ \(bdd :: BDD o) -> + BDD.unfoldHashable BDD.outSig bdd === bdd + +-- ------------------------------------------------------------------------ + case_support_false :: Assertion case_support_false = BDD.support BDD.false @?= IntSet.empty @@ -487,6 +643,15 @@ -- ------------------------------------------------------------------------ +case_numNodes :: Assertion +case_numNodes = do + let bdd, bdd1 :: BDD BDD.AscOrder + bdd = BDD.Branch 0 (BDD.Branch 1 BDD.false bdd1) (BDD.Branch 2 BDD.false bdd1) + bdd1 = BDD.Branch 3 BDD.true BDD.false + BDD.numNodes bdd @?= 6 + +-- ------------------------------------------------------------------------ + prop_restrict :: Property prop_restrict = forAllItemOrder $ \(_ :: Proxy o) -> @@ -613,8 +778,8 @@ case_restrictLaw_case_0 = (law BDD..&&. BDD.restrictLaw law a) @?= (law BDD..&&. a) where a, law :: BDD BDD.AscOrder - a = BDD.Branch 2 BDD.F BDD.T - law = BDD.Branch 1 (BDD.Branch 2 BDD.T BDD.F) (BDD.Branch 2 BDD.F BDD.T) + a = BDD.Branch 2 (BDD.Leaf False) (BDD.Leaf True) + law = BDD.Branch 1 (BDD.Branch 2 (BDD.Leaf True) (BDD.Leaf False)) (BDD.Branch 2 (BDD.Leaf False) (BDD.Leaf True)) prop_restrictLaw_true :: Property prop_restrictLaw_true = @@ -657,13 +822,13 @@ case_restrictLaw_case_1 :: Assertion case_restrictLaw_case_1 = do -- BDD.restrictLaw val a @?= BDD.restrictLaw val2 (BDD.restrictLaw val1 a) - BDD.restrictLaw val a @?= BDD.Branch 2 BDD.F BDD.T - BDD.restrictLaw val2 (BDD.restrictLaw val1 a) @?= BDD.Branch 1 BDD.T (Branch 2 BDD.F BDD.T) + BDD.restrictLaw val a @?= BDD.Branch 2 (BDD.Leaf False) (BDD.Leaf True) + BDD.restrictLaw val2 (BDD.restrictLaw val1 a) @?= BDD.Branch 1 (BDD.Leaf True) (Branch 2 (BDD.Leaf False) (BDD.Leaf True)) where a :: BDD BDD.AscOrder - a = Branch 2 BDD.F BDD.T -- x2 - val1 = BDD.Branch 1 BDD.F BDD.T -- x1 - val2 = BDD.Branch 1 (BDD.Branch 2 BDD.F BDD.T) BDD.T -- x1 ∨ x2 + a = Branch 2 (BDD.Leaf False) (BDD.Leaf True) -- x2 + val1 = BDD.Branch 1 (BDD.Leaf False) (BDD.Leaf True) -- x1 + val2 = BDD.Branch 1 (BDD.Branch 2 (BDD.Leaf False) (BDD.Leaf True)) (BDD.Leaf True) -- x1 ∨ x2 val = val1 BDD..&&. val2 -- x1 prop_restrictLaw_or_condition :: Property @@ -722,9 +887,9 @@ BDD.restrictLaw law a @?= b -- should be 'a'? where law, a :: BDD BDD.AscOrder - law = BDD.Branch 1 (BDD.Branch 2 BDD.F BDD.T) BDD.T -- x1 ∨ x2 - a = BDD.Branch 2 BDD.T BDD.F -- ¬x2 - b = BDD.Branch 1 BDD.F (BDD.Branch 2 BDD.T BDD.F) -- x1 ∧ ¬x2 + law = BDD.Branch 1 (BDD.Branch 2 (BDD.Leaf False) (BDD.Leaf True)) (BDD.Leaf True) -- x1 ∨ x2 + a = BDD.Branch 2 (BDD.Leaf True) (BDD.Leaf False) -- ¬x2 + b = BDD.Branch 1 (BDD.Leaf False) (BDD.Branch 2 (BDD.Leaf True) (BDD.Leaf False)) -- x1 ∧ ¬x2 case_restrictLaw_non_minimal_2 :: Assertion case_restrictLaw_non_minimal_2 = do @@ -732,9 +897,9 @@ BDD.restrictLaw law a @?= b -- should be 'a'? where law, a, b :: BDD BDD.AscOrder - law = BDD.Branch 1 BDD.T (BDD.Branch 2 BDD.F BDD.T) -- ¬x1 ∨ x2 - a = BDD.Branch 2 BDD.F BDD.T -- x2 - b = BDD.Branch 1 (BDD.Branch 2 BDD.F BDD.T) BDD.T -- x1 ∨ x2 + law = BDD.Branch 1 (BDD.Leaf True) (BDD.Branch 2 (BDD.Leaf False) (BDD.Leaf True)) -- ¬x1 ∨ x2 + a = BDD.Branch 2 (BDD.Leaf False) (BDD.Leaf True) -- x2 + b = BDD.Branch 1 (BDD.Branch 2 (BDD.Leaf False) (BDD.Leaf True)) (BDD.Leaf True) -- x1 ∨ x2 -- ------------------------------------------------------------------------ @@ -799,9 +964,9 @@ BDD.substSet (IntMap.singleton x m1) m @?= BDD.subst x m1 m where m :: BDD BDD.AscOrder - m = BDD.Branch 1 (BDD.Branch 2 BDD.T BDD.F) (BDD.Branch 2 BDD.F BDD.F) + m = BDD.Branch 1 (BDD.Branch 2 (BDD.Leaf True) (BDD.Leaf False)) (BDD.Branch 2 (BDD.Leaf False) (BDD.Leaf True)) x = 1 - m1 = BDD.Branch 1 BDD.T BDD.F + m1 = BDD.Branch 1 (BDD.Leaf True) (BDD.Leaf False) prop_substSet_same_vars :: Property prop_substSet_same_vars = @@ -825,8 +990,8 @@ prop_substSet_compose :: Property prop_substSet_compose = forAllItemOrder $ \(_ :: Proxy o) -> - forAll (liftM IntSet.fromList arbitrary) $ \xs -> - forAll (liftM IntSet.fromList arbitrary) $ \ys -> + forAll arbitrary $ \xs -> + forAll arbitrary $ \ys -> forAll (liftM IntMap.fromList $ mapM (\x -> (,) <$> pure x <*> arbitraryBDDOver ys) (IntSet.toList xs)) $ \s1 -> forAll (liftM IntMap.fromList $ mapM (\y -> (,) <$> pure y <*> arbitrary) (IntSet.toList ys)) $ \s2 -> forAll (arbitraryBDDOver xs) $ \(m :: BDD o) -> @@ -841,11 +1006,172 @@ -- ------------------------------------------------------------------------ +data MonotoneExpr a + = MVar a + | MAnd (MonotoneExpr a) (MonotoneExpr a) + | MOr (MonotoneExpr a) (MonotoneExpr a) + | MConst Bool + deriving (Show) + +arbitraryMonotoneExpr :: forall a. Gen a -> Gen (MonotoneExpr a) +arbitraryMonotoneExpr gen = sized f + where + f :: Int -> Gen (MonotoneExpr a) + f n = oneof $ + [ liftM MConst arbitrary + , liftM MVar gen + ] + ++ + concat + [ [liftM2 MAnd sub sub, liftM2 MOr sub sub] + | n > 0, let sub = f (n `div` 2) + ] + +evalMonotoneExpr :: ItemOrder a => (b -> BDD a) -> MonotoneExpr b -> BDD a +evalMonotoneExpr f = g + where + g (MVar a) = f a + g (MConst v) = BDD.Leaf v + g (MOr a b) = g a BDD..||. g b + g (MAnd a b) = g a BDD..&&. g b + +forAllMonotonicFunction :: forall o prop. (ItemOrder o, Testable prop) => IntSet -> ((BDD o -> BDD o) -> prop) -> Property +forAllMonotonicFunction xs k = + forAll (arbitraryMonotoneExpr (elements (Nothing : map Just (IntSet.toList xs)))) $ \e -> do + let f :: BDD o -> BDD o + f x = evalMonotoneExpr g e + where + g Nothing = x + g (Just v) = BDD.var v + in k f + +prop_lfp_is_fixed_point :: Property +prop_lfp_is_fixed_point = + forAll arbitrary $ \(xs :: IntSet) -> + forAllItemOrder $ \(_ :: Proxy o) -> + forAllMonotonicFunction xs $ \(f :: BDD o -> BDD o) -> do + let a = BDD.lfp f + in counterexample (show a) $ f a === a + + +prop_gfp_is_fixed_point :: Property +prop_gfp_is_fixed_point = + forAll arbitrary $ \(xs :: IntSet) -> + forAllItemOrder $ \(_ :: Proxy o) -> + forAllMonotonicFunction xs $ \(f :: BDD o -> BDD o) -> do + let a = BDD.gfp f + in counterexample (show a) $ f a === a + +prop_lfp_imply_gfp :: Property +prop_lfp_imply_gfp = + forAll arbitrary $ \(xs :: IntSet) -> + forAllItemOrder $ \(_ :: Proxy o) -> + forAllMonotonicFunction xs $ \(f :: BDD o -> BDD o) -> do + let a = BDD.lfp f + b = BDD.gfp f + in counterexample (show (a, b)) $ (a BDD..=>. b) === BDD.true + +-- ------------------------------------------------------------------------ + +prop_anySat :: Property +prop_anySat = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrary $ \(bdd :: BDD o) -> + case BDD.anySat bdd of + Just p -> counterexample (show p) $ BDD.evaluate (p IntMap.!) bdd + Nothing -> bdd === BDD.Leaf False + +prop_allSat :: Property +prop_allSat = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntSet $ \xs -> + forAll (arbitraryBDDOver xs) $ \(bdd :: BDD o) -> + let ps = BDD.allSat bdd + in null ps === (bdd == BDD.Leaf False) + .&&. + conjoin [counterexample (show p) $ BDD.evaluate (p IntMap.!) bdd | p <- ps] + +prop_anySatComplete :: Property +prop_anySatComplete = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrary $ \xs -> + forAll (arbitraryBDDOver xs) $ \(bdd :: BDD o) -> + case BDD.anySatComplete xs bdd of + Just p -> counterexample (show p) $ + IntMap.keysSet p === xs + .&&. + BDD.evaluate (p IntMap.!) bdd + Nothing -> bdd === BDD.Leaf False + +prop_allSatComplete :: Property +prop_allSatComplete = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntSet $ \xs -> + forAll (arbitraryBDDOver xs) $ \(bdd :: BDD o) -> + let ps = BDD.allSatComplete xs bdd + qs = [q | q <- foldM (\m x -> [IntMap.insert x v m | v <- [False, True]]) IntMap.empty (IntSet.toList xs) + , BDD.evaluate (q IntMap.!) bdd] + in conjoin [counterexample (show p) (IntMap.keysSet p === xs) | p <- ps] + .&&. + Set.fromList ps === Set.fromList qs + +prop_countSat_allSatComplete :: Property +prop_countSat_allSatComplete = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntSet $ \xs -> + forAll (arbitraryBDDOver xs) $ \(bdd :: BDD o) -> + let ps = BDD.allSatComplete xs bdd + n = BDD.countSat xs bdd + in counterexample (show n) $ + if bdd == BDD.Leaf False then + n === 0 + else + -- Note that the number of partial assignments is smaller than the number of total assignments + (n > 0) .&&. n === length ps + +prop_uniformSatM :: Property +prop_uniformSatM = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll (arbitrarySmallIntSet `suchThat` ((>= 2) . IntSet.size)) $ \xs -> + forAll (arbitraryBDDOver xs `suchThat` ((>= (2::Integer)) . BDD.countSat xs)) $ \(bdd :: BDD o) -> + forAll arbitrary $ \(seed :: Vector Word32) -> + let m :: Integer + m = BDD.countSat xs bdd + n = 1000 + samples = runST $ do + gen <- Rand.initialize seed + replicateM n $ BDD.uniformSatM xs bdd gen + hist_actual = Map.fromListWith (+) [(s, 1 :: Double) | s <- samples] + hist_expected = [(s, fromIntegral n / fromIntegral m :: Double) | s <- BDD.allSatComplete xs bdd] + chi_sq = sum [(Map.findWithDefault 0 s hist_actual - cnt) ** 2 / cnt | (s, cnt) <- hist_expected] + threshold = complQuantile (chiSquared (fromIntegral m - 1)) 0.0001 + in counterexample (show hist_actual ++ " /= " ++ show (Map.fromList hist_expected)) $ + and [BDD.evaluate (a IntMap.!) bdd | a <- Map.keys hist_actual] + .&&. + counterexample ("χ² = " ++ show chi_sq ++ " >= " ++ show threshold) (chi_sq < threshold) + +-- ------------------------------------------------------------------------ + prop_toGraph_fromGraph :: Property prop_toGraph_fromGraph = do forAllItemOrder $ \(_ :: Proxy o) -> forAll arbitrary $ \(a :: BDD o) -> BDD.fromGraph (BDD.toGraph a) === a + +-- ------------------------------------------------------------------------ + +arbitrarySmallIntSet :: Gen IntSet +arbitrarySmallIntSet = do + n <- choose (0, 12) + liftM IntSet.fromList $ replicateM n arbitrary + +arbitrarySmallIntMap :: Arbitrary a => Gen (IntMap a) +arbitrarySmallIntMap = do + n <- choose (0, 12) + liftM IntMap.fromList $ replicateM n $ do + k <- arbitrary + v <- arbitrary + return (k, v) -- ------------------------------------------------------------------------ diff --git a/test/TestZDD.hs b/test/TestZDD.hs --- a/test/TestZDD.hs +++ b/test/TestZDD.hs @@ -3,25 +3,34 @@ {-# LANGUAGE TemplateHaskell #-} module TestZDD (zddTestGroup) where +import Control.DeepSeq import Control.Monad +import Control.Monad.ST +import Data.IntMap (IntMap) +import qualified Data.IntMap as IntMap import Data.IntSet (IntSet) import qualified Data.IntSet as IntSet +import Data.IORef import Data.List import qualified Data.Map.Strict as Map import Data.Proxy import Data.Set (Set) import qualified Data.Set as Set +import Data.Vector (Vector) +import Data.Word import qualified GHC.Exts as Exts import Statistics.Distribution import Statistics.Distribution.ChiSquared (chiSquared) +import System.IO.Unsafe import qualified System.Random.MWC as Rand import Test.QuickCheck.Function (apply) -import qualified Test.QuickCheck.Monadic as QM +import Test.QuickCheck.Instances.Vector () import Test.Tasty import Test.Tasty.HUnit import Test.Tasty.QuickCheck import Test.Tasty.TH +import Data.DecisionDiagram.BDD.Internal.ItemOrder (OrderedItem (..)) import Data.DecisionDiagram.ZDD (ZDD (..), ItemOrder (..)) import qualified Data.DecisionDiagram.ZDD as ZDD @@ -31,7 +40,7 @@ instance ZDD.ItemOrder a => Arbitrary (ZDD a) where arbitrary = do - vars <- liftM (sortBy (ZDD.compareItem (Proxy :: Proxy a)) . IntSet.toList . IntSet.fromList) arbitrary + vars <- liftM (sortBy (ZDD.compareItem (Proxy :: Proxy a)) . IntSet.toList) arbitrary let f vs n = oneof $ [ return ZDD.empty , return ZDD.base @@ -55,6 +64,13 @@ | (p0', p1') <- shrink (p0, p1), p1' /= ZDD.empty ] +arbitraryMember :: ZDD.ItemOrder a => ZDD a -> Gen IntSet +arbitraryMember zdd = do + (seed :: Vector Word32) <- arbitrary + return $ runST $ do + gen <- Rand.initialize seed + ZDD.uniformM zdd gen + -- ------------------------------------------------------------------------ -- Union -- ------------------------------------------------------------------------ @@ -292,7 +308,7 @@ prop_singleton :: Property prop_singleton = forAllItemOrder $ \(_ :: Proxy o) -> - forAll (liftM IntSet.fromList arbitrary) $ \xs -> + forAll arbitrary $ \xs -> let a :: ZDD o a = ZDD.singleton xs in counterexample (show a) $ ZDD.toSetOfIntSets a === Set.singleton xs @@ -303,7 +319,7 @@ forAll arbitrary $ \xs -> let a :: ZDD o a = ZDD.subsets xs - in counterexample (show a) $ forAll (liftM IntSet.fromList (sublistOf (IntSet.toList xs))) $ \ys -> + in counterexample (show a) $ forAll (subsetOf xs) $ \ys -> ys `ZDD.member` a prop_subsets_member_empty :: Property @@ -330,6 +346,110 @@ a = ZDD.subsets xs in counterexample (show a) $ ZDD.size a === (2 :: Integer) ^ (IntSet.size xs) +prop_subsetsAtLeast :: Property +prop_subsetsAtLeast = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntMap $ \(xs :: IntMap Integer) -> + forAll arbitrary $ \k -> + let a :: ZDD o + a = ZDD.subsetsAtLeast xs k + in counterexample (show a) $ + if ZDD.null a then + property (k > sum [max 0 w | (_,w) <- IntMap.toList xs]) + else + forAll (arbitraryMember a) $ \ys -> + (ys `IntSet.isSubsetOf` IntMap.keysSet xs) + .&&. + sum [xs IntMap.! y | y <- IntSet.toList ys] >= k + +prop_subsetsAtMost :: Property +prop_subsetsAtMost = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntMap $ \(xs :: IntMap Integer) -> + forAll arbitrary $ \k -> + let a :: ZDD o + a = ZDD.subsetsAtMost xs k + in counterexample (show a) $ + if ZDD.null a then + property (k < sum [min 0 w | (_,w) <- IntMap.toList xs]) + else + forAll (arbitraryMember a) $ \ys -> + (ys `IntSet.isSubsetOf` IntMap.keysSet xs) + .&&. + sum [xs IntMap.! y | y <- IntSet.toList ys] <= k + +prop_subsetsExactly :: Property +prop_subsetsExactly = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntMap $ \(xs :: IntMap Integer) -> + forAll arbitrary $ \k -> + let a :: ZDD o + a = ZDD.subsetsExactly xs k + in counterexample (show a) $ + if ZDD.null a then + property True + else + forAll (arbitraryMember a) $ \ys -> + (ys `IntSet.isSubsetOf` IntMap.keysSet xs) + .&&. + sum [xs IntMap.! y | y <- IntSet.toList ys] === k + +prop_subsetsExactly_2 :: Property +prop_subsetsExactly_2 = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntMap $ \(xs :: IntMap Integer) -> + forAll (gen xs) $ \(ys, k) -> + let a :: ZDD o + a = ZDD.subsetsExactly xs k + in counterexample (show a) $ ys `ZDD.member` a + where + gen xs = do + ys <- sublistOf (IntMap.toList xs) + return (IntSet.fromList [y | (y,_) <- ys], sum [w | (_,w) <- ys]) + +prop_subsetsExactlyIntegral :: Property +prop_subsetsExactlyIntegral = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrarySmallIntMap $ \(xs :: IntMap Integer) -> + forAll arbitrary $ \k -> + (ZDD.subsetsExactlyIntegral xs k :: ZDD o) === ZDD.subsetsExactly xs k + +prop_combinations_are_combinations :: Property +prop_combinations_are_combinations = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrary $ \xs -> + forAll arbitrary $ \(NonNegative k) -> + let a :: ZDD o + a = ZDD.combinations xs k + in counterexample (show a) $ + not (ZDD.null a) + ==> + (forAll (arbitraryMember a) $ \ys -> (ys `IntSet.isSubsetOf` xs) .&&. (IntSet.size ys === k)) + +prop_combinations_size :: Property +prop_combinations_size = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrary $ \xs -> + forAll arbitrary $ \(NonNegative k) -> + let a :: ZDD o + a = ZDD.combinations xs k + n = toInteger $ IntSet.size xs + in counterexample (show a) $ ZDD.size a === (product [(n - toInteger k + 1)..n] `div` (product [1..toInteger k])) + +case_toList_lazyness :: Assertion +case_toList_lazyness = do + let xss :: ZDD ZDD.AscOrder + xss = ZDD.subsets (IntSet.fromList [1..128]) + deepseq (take 100 (Exts.toList xss)) $ return () + +prop_toList_sorted :: Property +prop_toList_sorted = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrary $ \(xss :: ZDD o) -> + let yss :: [[OrderedItem o]] + yss = map (sort . map OrderedItem . IntSet.toList) $ take 100 $ Exts.toList xss + in yss === sort yss + prop_toList_fromList :: Property prop_toList_fromList = forAllItemOrder $ \(_ :: Proxy o) -> @@ -349,7 +469,7 @@ prop_toSetOfIntSets_fromSetOfIntSets :: Property prop_toSetOfIntSets_fromSetOfIntSets = forAllItemOrder $ \(_ :: Proxy o) -> - forAll (liftM (Set.fromList . map IntSet.fromList) arbitrary) $ \xss -> + forAll arbitrary $ \xss -> let a :: ZDD o a = ZDD.fromSetOfIntSets xss in counterexample (show a) $ ZDD.toSetOfIntSets a === xss @@ -365,14 +485,14 @@ prop_insert = forAllItemOrder $ \(_ :: Proxy o) -> forAll arbitrary $ \(a :: ZDD o) -> - forAll (liftM IntSet.fromList arbitrary) $ \xs -> + forAll arbitrary $ \xs -> ZDD.toSetOfIntSets (ZDD.insert xs a) === Set.insert xs (ZDD.toSetOfIntSets a) prop_insert_idempotent :: Property prop_insert_idempotent = forAllItemOrder $ \(_ :: Proxy o) -> forAll arbitrary $ \(a :: ZDD o) -> - forAll (liftM IntSet.fromList arbitrary) $ \xs -> + forAll arbitrary $ \xs -> let b = ZDD.insert xs a in counterexample (show b) $ ZDD.insert xs b === b @@ -380,14 +500,14 @@ prop_delete = forAllItemOrder $ \(_ :: Proxy o) -> forAll arbitrary $ \(a :: ZDD o) -> - forAll (liftM IntSet.fromList $ sublistOf (IntSet.toList (ZDD.flatten a))) $ \xs -> + forAll (subsetOf (ZDD.flatten a)) $ \xs -> ZDD.toSetOfIntSets (ZDD.delete xs a) === Set.delete xs (ZDD.toSetOfIntSets a) prop_delete_idempotent :: Property prop_delete_idempotent = forAllItemOrder $ \(_ :: Proxy o) -> forAll arbitrary $ \(a :: ZDD o) -> - forAll (liftM IntSet.fromList $ sublistOf (IntSet.toList (ZDD.flatten a))) $ \xs -> + forAll (subsetOf (ZDD.flatten a)) $ \xs -> let b = ZDD.delete xs a in counterexample (show b) $ ZDD.delete xs b === b @@ -462,7 +582,7 @@ prop_member_2 = forAllItemOrder $ \(_ :: Proxy o) -> forAll arbitrary $ \(a :: ZDD o) -> - forAll (liftM IntSet.fromList $ sublistOf (IntSet.toList (ZDD.flatten a))) $ \s2 -> + forAll (subsetOf (ZDD.flatten a)) $ \s2 -> (s2 `ZDD.member` a) === (s2 `Set.member` ZDD.toSetOfIntSets a) prop_size :: Property @@ -507,6 +627,13 @@ forAll arbitrary $ \(a :: ZDD o, b) -> ZDD.disjoint a b === ZDD.null (a `ZDD.intersection` b) +case_numNodes :: Assertion +case_numNodes = do + let bdd, bdd1 :: ZDD ZDD.AscOrder + bdd = ZDD.Branch 0 (ZDD.Branch 1 ZDD.empty bdd1) (ZDD.Branch 2 ZDD.empty bdd1) + bdd1 = ZDD.Branch 3 ZDD.empty ZDD.base + ZDD.numNodes bdd @?= 6 + prop_flatten :: Property prop_flatten = forAllItemOrder $ \(_ :: Proxy o) -> @@ -517,20 +644,20 @@ prop_uniformM = forAllItemOrder $ \(_ :: Proxy o) -> forAll (arbitrary `suchThat` ((>= (2::Integer)) . ZDD.size)) $ \(a :: ZDD o) -> - QM.monadicIO $ do - gen <- QM.run Rand.create + forAll arbitrary $ \(seed :: Vector Word32) -> let m :: Integer m = ZDD.size a n = 1000 - samples <- QM.run $ replicateM n $ ZDD.uniformM a gen - let hist_actual = Map.fromListWith (+) [(s, 1) | s <- samples] + samples = runST $ do + gen <- Rand.initialize seed + replicateM n $ ZDD.uniformM a gen + hist_actual = Map.fromListWith (+) [(s, 1) | s <- samples] hist_expected = [(s, fromIntegral n / fromIntegral m) | s <- ZDD.toListOfIntSets a] chi_sq = sum [(Map.findWithDefault 0 s hist_actual - cnt) ** 2 / cnt | (s, cnt) <- hist_expected] - threshold = complQuantile (chiSquared (fromIntegral m - 1)) 0.001 - QM.monitor $ counterexample $ show hist_actual ++ " /= " ++ show (Map.fromList hist_expected) - QM.assert $ and [xs `ZDD.member` a | xs <- Map.keys hist_actual] - QM.monitor $ counterexample $ "χ² = " ++ show chi_sq ++ " >= " ++ show threshold - QM.assert $ chi_sq < threshold + threshold = complQuantile (chiSquared (fromIntegral m - 1)) 0.0001 + in counterexample (show hist_actual ++ " /= " ++ show (Map.fromList hist_expected)) $ + and [xs `ZDD.member` a | xs <- Map.keys hist_actual] .&&. + counterexample ("χ² = " ++ show chi_sq ++ " >= " ++ show threshold) (chi_sq < threshold) prop_findMinSum :: Property prop_findMinSum = @@ -562,6 +689,55 @@ -- ------------------------------------------------------------------------ +case_fold_laziness :: Assertion +case_fold_laziness = do + let bdd :: ZDD ZDD.AscOrder + bdd = ZDD.Branch 0 (ZDD.Branch 1 ZDD.Empty ZDD.Base) (ZDD.Branch 2 ZDD.Empty ZDD.Base) + f x lo _hi = + if x == 2 then + error "unused value should not be evaluated" + else + lo + seq (ZDD.fold f id bdd) $ return () + +case_fold'_strictness :: Assertion +case_fold'_strictness = do + ref <- newIORef False + let bdd :: ZDD ZDD.AscOrder + bdd = ZDD.Branch 0 (ZDD.Branch 1 ZDD.Empty ZDD.Base) (ZDD.Branch 2 ZDD.Empty ZDD.Base) + f x lo _hi = unsafePerformIO $ do + when (x==2) $ writeIORef ref True + return lo + seq (ZDD.fold' f id bdd) $ do + flag <- readIORef ref + assertBool "unused value should be evaluated" flag + +prop_fold_inSig :: Property +prop_fold_inSig = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrary $ \(zdd :: ZDD o) -> + ZDD.fold (\x lo hi -> ZDD.inSig (ZDD.SBranch x lo hi)) (ZDD.inSig . ZDD.SLeaf) zdd + === + zdd + +prop_fold'_inSig :: Property +prop_fold'_inSig = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrary $ \(zdd :: ZDD o) -> + ZDD.fold' (\x lo hi -> ZDD.inSig (ZDD.SBranch x lo hi)) (ZDD.inSig . ZDD.SLeaf) zdd + === + zdd + +-- ------------------------------------------------------------------------ + +prop_unfoldHashable_outSig :: Property +prop_unfoldHashable_outSig = + forAllItemOrder $ \(_ :: Proxy o) -> + forAll arbitrary $ \(zdd :: ZDD o) -> + ZDD.unfoldHashable ZDD.outSig zdd === zdd + +-- ------------------------------------------------------------------------ + prop_toGraph_fromGraph :: Property prop_toGraph_fromGraph = do forAllItemOrder $ \(_ :: Proxy o) -> @@ -674,6 +850,19 @@ ZDD.union tmp tmp2 @?= ZDD.union tmp ZDD.base -- 3. DIFF ZDD.difference tmp tmp2 @?= ZDD.fromListOfIntSets (map IntSet.fromList [[1], [2]]) + +-- ------------------------------------------------------------------------ + +subsetOf :: IntSet -> Gen IntSet +subsetOf = liftM IntSet.fromList . sublistOf . IntSet.toList + +arbitrarySmallIntMap :: Arbitrary a => Gen (IntMap a) +arbitrarySmallIntMap = do + n <- choose (0, 12) + liftM IntMap.fromList $ replicateM n $ do + k <- arbitrary + v <- arbitrary + return (k, v) -- ------------------------------------------------------------------------ diff --git a/test/doctests.hs b/test/doctests.hs new file mode 100644 --- /dev/null +++ b/test/doctests.hs @@ -0,0 +1,10 @@ +module Main (main) where + +import Test.DocTest + +main :: IO () +main = doctest + [ "-isrc" + , "src/Data/DecisionDiagram/BDD.hs" + , "src/Data/DecisionDiagram/ZDD.hs" + ]