{-# OPTIONS_GHC -Wall #-} {-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-} ----------------------------------------------------------------------------- -- | -- Module : FourierMotzkin -- Copyright : (c) Masahiro Sakai 2011 -- License : BSD-style -- -- Maintainer : masahiro.sakai@gmail.com -- Stability : provisional -- Portability : non-portable (MultiParamTypeClasses, FunctionalDependencies) -- -- Naïve implementation of Fourier-Motzkin Variable Elimination -- -- Reference: -- -- * -- ----------------------------------------------------------------------------- module FourierMotzkin ( module Data.Expr , module Data.Formula , Lit (..) , eliminateQuantifiers , solve -- FIXME , termR , Rat , collectBounds , boundConditions , evalBounds , simplify , constraintsToDNF ) where import Control.Monad import Data.List import Data.Maybe import Data.Ratio import qualified Data.IntMap as IM import qualified Data.IntSet as IS import Data.Expr import Data.Formula import Data.Linear import qualified Data.LA as LA import qualified Data.Interval as Interval -- --------------------------------------------------------------------------- type ExprZ = LA.Expr Integer -- | (t,c) represents t/c, and c must be >0. type Rat = (ExprZ, Integer) evalRat :: Model Rational -> Rat -> Rational evalRat model (e, d) = LA.lift1 1 (model IM.!) (LA.mapCoeff fromIntegral e) / (fromIntegral d) -- | Literal data Lit = Nonneg ExprZ | Pos ExprZ deriving (Show, Eq, Ord) instance Variables Lit where vars (Pos t) = vars t vars (Nonneg t) = vars t instance Complement Lit where notB (Pos t) = Nonneg (lnegate t) notB (Nonneg t) = Pos (lnegate t) -- 制約集合の単純化 -- It returns Nothing when a inconsistency is detected. simplify :: [Lit] -> Maybe [Lit] simplify = fmap concat . mapM f where f :: Lit -> Maybe [Lit] f lit@(Pos e) = case LA.asConst e of Just x -> guard (x > 0) >> return [] Nothing -> return [lit] f lit@(Nonneg e) = case LA.asConst e of Just x -> guard (x >= 0) >> return [] Nothing -> return [lit] -- --------------------------------------------------------------------------- atomR :: RelOp -> Expr Rational -> Expr Rational -> Maybe (DNF Lit) atomR op a b = do a' <- termR a b' <- termR b return $ case op of Le -> DNF [[a' `leR` b']] Lt -> DNF [[a' `ltR` b']] Ge -> DNF [[a' `geR` b']] Gt -> DNF [[a' `gtR` b']] Eql -> DNF [[a' `leR` b', a' `geR` b']] NEq -> DNF [[a' `ltR` b'], [a' `gtR` b']] termR :: Expr Rational -> Maybe Rat termR (Const c) = return (LA.constant (numerator c), denominator c) termR (Var v) = return (LA.var v, 1) termR (a :+: b) = do (t1,c1) <- termR a (t2,c2) <- termR b return (c2 .*. t1 .+. c1 .*. t2, c1*c2) termR (a :*: b) = do (t1,c1) <- termR a (t2,c2) <- termR b msum [ do{ c <- LA.asConst t1; return (c .*. t2, c1*c2) } , do{ c <- LA.asConst t2; return (c .*. t1, c1*c2) } ] termR (a :/: b) = do (t1,c1) <- termR a (t2,c2) <- termR b c3 <- LA.asConst t2 guard $ c3 /= 0 return (c2 .*. t1, c1*c3) leR, ltR, geR, gtR :: Rat -> Rat -> Lit leR (e1,c) (e2,d) = Nonneg $ normalizeExprR $ c .*. e2 .-. d .*. e1 ltR (e1,c) (e2,d) = Pos $ normalizeExprR $ c .*. e2 .-. d .*. e1 geR = flip leR gtR = flip gtR normalizeExprR :: ExprZ -> ExprZ normalizeExprR e = LA.mapCoeff (`div` d) e where d = abs $ gcd' $ map fst $ LA.terms e -- --------------------------------------------------------------------------- {- (ls1,ls2,us1,us2) represents { x | ∀(M,c)∈ls1. M/c≤x, ∀(M,c)∈ls2. M/c [Lit] -> DNF Lit eliminate v xs = DNF [rest] .&&. boundConditions bnd where (bnd, rest) = collectBounds v xs collectBounds :: Var -> [Lit] -> (BoundsR, [Lit]) collectBounds v = foldr phi (([],[],[],[]),[]) where phi :: Lit -> (BoundsR, [Lit]) -> (BoundsR, [Lit]) phi lit@(Nonneg t) x = f False lit t x phi lit@(Pos t) x = f True lit t x f :: Bool -> Lit -> ExprZ -> (BoundsR, [Lit]) -> (BoundsR, [Lit]) f strict lit t (bnd@(ls1,ls2,us1,us2), xs) = case LA.extract v t of (c,t') -> case c `compare` 0 of EQ -> (bnd, lit : xs) GT -> if strict then ((ls1, (lnegate t', c) : ls2, us1, us2), xs) -- 0 < cx + M ⇔ -M/c < x else (((lnegate t', c) : ls1, ls2, us1, us2), xs) -- 0 ≤ cx + M ⇔ -M/c ≤ x LT -> if strict then ((ls1, ls2, us1, (t', negate c) : us2), xs) -- 0 < cx + M ⇔ x < M/-c else ((ls1, ls2, (t', negate c) : us1, us2), xs) -- 0 ≤ cx + M ⇔ x ≤ M/-c boundConditions :: BoundsR -> DNF Lit boundConditions (ls1, ls2, us1, us2) = DNF $ maybeToList $ simplify $ [ x `leR` y | x <- ls1, y <- us1 ] ++ [ x `ltR` y | x <- ls1, y <- us2 ] ++ [ x `ltR` y | x <- ls2, y <- us1 ] ++ [ x `ltR` y | x <- ls2, y <- us2 ] eliminateQuantifiers :: Formula Rational -> Maybe (DNF Lit) eliminateQuantifiers = f where f T = return true f F = return false f (Atom (Rel a op b)) = atomR op a b f (And a b) = liftM2 (.&&.) (f a) (f b) f (Or a b) = liftM2 (.||.) (f a) (f b) f (Not a) = f (pushNot a) f (Imply a b) = f (Or (Not a) b) f (Equiv a b) = f (And (Imply a b) (Imply b a)) f (Forall v a) = do dnf <- f (Exists v (pushNot a)) return (notB dnf) f (Exists v a) = do dnf <- f a return $ orB [eliminate v xs | xs <- unDNF dnf] solve :: Formula Rational -> SatResult Rational solve formula = case eliminateQuantifiers formula of Nothing -> Unknown Just dnf -> case msum [solve' vs xs | xs <- unDNF dnf] of Nothing -> Unsat Just m -> Sat m where vs = IS.toList (vars formula) solve' :: [Var] -> [Lit] -> Maybe (Model Rational) solve' vs xs = simplify xs >>= go vs where go [] [] = return IM.empty go [] _ = mzero go (v:vs) ys = msum (map f (unDNF (boundConditions bnd))) where (bnd, rest) = collectBounds v ys f zs = do model <- go vs (zs ++ rest) val <- Interval.pickup (evalBounds model bnd) return $ IM.insert v val model evalBounds :: Model Rational -> BoundsR -> Interval.Interval Rational evalBounds model (ls1,ls2,us1,us2) = foldl' Interval.intersection Interval.univ $ [ Interval.interval (Just (True, evalRat model x)) Nothing | x <- ls1 ] ++ [ Interval.interval (Just (False, evalRat model x)) Nothing | x <- ls2 ] ++ [ Interval.interval Nothing (Just (True, evalRat model x)) | x <- us1 ] ++ [ Interval.interval Nothing (Just (False, evalRat model x)) | x <- us2 ] -- --------------------------------------------------------------------------- constraintsToDNF :: [LA.Atom Rational] -> DNF Lit constraintsToDNF = andB . map constraintToDNF constraintToDNF :: LA.Atom Rational -> DNF Lit constraintToDNF (LA.Atom a op b) = DNF $ case op of Eql -> [[Nonneg c, Nonneg (lnegate c)]] NEq -> [[Pos c], [Pos (lnegate c)]] Ge -> [[Nonneg c]] Le -> [[Nonneg (lnegate c)]] Gt -> [[Pos c]] Lt -> [[Pos (lnegate c)]] where c = normalize (a .-. b) normalize :: LA.Expr Rational -> ExprZ normalize e = LA.mapCoeff (round . (*c)) e where c = fromIntegral $ foldl' lcm 1 ds ds = [denominator c | (c,v) <- LA.terms e] -- --------------------------------------------------------------------------- gcd' :: [Integer] -> Integer gcd' [] = 1 gcd' xs = foldl1' gcd xs -- --------------------------------------------------------------------------- {- 7x + 12y + 31z = 17 3x + 5y + 14z = 7 1 ≤ x ≤ 40 -50 ≤ y ≤ 50 satisfiable in R but unsatisfiable in Z -} test1 :: Formula Rational test1 = c1 .&&. c2 .&&. c3 .&&. c4 where x = Var 0 y = Var 1 z = Var 2 c1 = 7*x + 12*y + 31*z .==. 17 c2 = 3*x + 5*y + 14*z .==. 7 c3 = 1 .<=. x .&&. x .<=. 40 c4 = (-50) .<=. y .&&. y .<=. 50 test1' :: [LA.Atom Rational] test1' = [c1, c2] ++ c3 ++ c4 where x = LA.var 0 y = LA.var 1 z = LA.var 2 c1 = 7.*.x .+. 12.*.y .+. 31.*.z .==. LA.constant 17 c2 = 3.*.x .+. 5.*.y .+. 14.*.z .==. LA.constant 7 c3 = [LA.constant 1 .<=. x, x .<=. LA.constant 40] c4 = [LA.constant (-50) .<=. y, y .<=. LA.constant 50] {- 27 ≤ 11x+13y ≤ 45 -10 ≤ 7x-9y ≤ 4 satisfiable in R but unsatisfiable in Z -} test2 :: Formula Rational test2 = c1 .&&. c2 where x = Var 0 y = Var 1 t1 = 11*x + 13*y t2 = 7*x - 9*y c1 = 27 .<=. t1 .&&. t1 .<=. 45 c2 = (-10) .<=. t2 .&&. t2 .<=. 4 -- ---------------------------------------------------------------------------