From 823bf57183982d5850ef12c0d4a9aa462ef1e357 Mon Sep 17 00:00:00 2001 From: Niklas Gruhn Date: Sun, 2 Apr 2023 16:41:12 +0200 Subject: [PATCH] Construct model in Fourier-Motzkin --- SMT.cabal | 1 + src/Theory/LinearArithmatic/Constraint.hs | 82 +++++-- src/Theory/LinearArithmatic/FourierMotzkin.hs | 204 ++++++++++-------- src/Theory/LinearArithmatic/Simplex.hs | 30 +-- src/Utils.hs | 4 +- test/Test.hs | 4 +- test/TestLinearArithmatic.hs | 36 ++-- 7 files changed, 221 insertions(+), 140 deletions(-) diff --git a/SMT.cabal b/SMT.cabal index e1a3f0d..b2a9066 100644 --- a/SMT.cabal +++ b/SMT.cabal @@ -59,3 +59,4 @@ test-suite test hs-source-dirs: test default-language: Haskell2010 default-extensions: StrictData + , ScopedTypeVariables diff --git a/src/Theory/LinearArithmatic/Constraint.hs b/src/Theory/LinearArithmatic/Constraint.hs index 8c1a05b..7558050 100644 --- a/src/Theory/LinearArithmatic/Constraint.hs +++ b/src/Theory/LinearArithmatic/Constraint.hs @@ -6,6 +6,7 @@ import Control.Monad (guard) import qualified Data.IntMap.Merge.Lazy as MM import Data.List (intercalate) import Data.Ratio (denominator, numerator) +import Debug.Trace type Var = Int @@ -26,9 +27,18 @@ type Assignment = M.IntMap Rational -- | -- For example, constraint such as `3x + y + 3 <= 0` is represented as: -- --- (LessEquals, AffineExpr 3 (3x+y)) +-- (AffineExpr 3 (3x+y), LessEquals) -- -type Constraint = (ConstraintRelation, AffineExpr) +type Constraint = (AffineExpr, ConstraintRelation) + +{-| + Remove variables with zero coefficient. + TODO: ensure that constraints are always normalized. +-} +normalize :: Constraint -> Constraint +normalize (AffineExpr constant expr, rel) = + (AffineExpr constant $ M.filter (/= 0) expr, rel) + -- | Map from variable IDs to coefficients type LinearExpr = M.IntMap Rational @@ -60,10 +70,13 @@ instance Show AffineExpr where intercalate " + " (fmap show_var (M.toList coeff_map) ++ [show_signed constant]) varsIn :: Constraint -> S.IntSet -varsIn (_, AffineExpr _ coeff_map) = M.keysSet coeff_map +varsIn (AffineExpr _ coeff_map, _) = M.keysSet coeff_map appearsIn :: Var -> Constraint -> Bool -appearsIn var = M.member var . getCoeffMap . snd +appearsIn var = M.member var . getCoeffMap . fst + +modifyConstant :: (Rational -> Rational) -> AffineExpr -> AffineExpr +modifyConstant f (AffineExpr constant coeffs) = AffineExpr (f constant) coeffs {-| Solve a constraint for a variable. For example solving @@ -78,7 +91,7 @@ appearsIn var = M.member var . getCoeffMap . snd expression or the coefficient is zero. -} solveFor :: Constraint -> Var -> Maybe Constraint -solveFor (rel, AffineExpr constant coeffs) x = do +solveFor (AffineExpr constant coeffs, rel) x = do coeff_of_x <- M.lookup x coeffs guard (coeff_of_x /= 0) @@ -86,7 +99,7 @@ solveFor (rel, AffineExpr constant coeffs) x = do return -- flip the relation: x >= -10/3y + 1/3 - $ (rel', ) + $ (, rel') -- also divide constant term: x <= -10/3y + 1/3 $ AffineExpr (- constant / coeff_of_x) -- divide coefficients: x <= -10/3y - 1 @@ -94,13 +107,52 @@ solveFor (rel, AffineExpr constant coeffs) x = do -- get x on the other side: -3x <= 10y - 1 $ M.delete x coeffs -eval :: Assignment -> AffineExpr -> Rational -eval assignment (AffineExpr constant coeff_map) = - (constant +) $ sum $ - MM.merge - MM.dropMissing - MM.dropMissing - (MM.zipWithMatched (const (*))) - assignment - coeff_map +substitute :: Assignment -> AffineExpr -> AffineExpr +substitute partial_assignment (AffineExpr constant coeff_map) = + let + constant' = (constant +) $ sum $ + MM.merge + MM.dropMissing + MM.dropMissing + (MM.zipWithMatched (const (*))) + partial_assignment + coeff_map + + coeff_map' = coeff_map `M.withoutKeys` M.keysSet partial_assignment + in + AffineExpr constant' coeff_map' +{-| + Evaluate an expression by substituting all variables + using the given assignment. Returns `Nothing` if the + assignment is partial. +-} +eval :: Assignment -> AffineExpr -> Maybe Rational +eval assignment expr + | M.null coeff_map = Just constant + | otherwise = traceShow coeff_map Nothing + where + AffineExpr constant coeff_map = substitute assignment expr + + +-- linearDependent :: Constraint -> Constraint -> Bool +-- linearDependent vs ws = coeff /= 0 && V.map (* coeff) vs == ws +-- where +-- div' x y +-- | y == 0 = 0 +-- | otherwise = x/y +-- coeff = +-- fromMaybe 0 $ V.find (/=0) $ V.zipWith div' ws vs + +{-| + Checks if the assignment satisfies the constraints. +-} +isModel :: Assignment -> Constraint -> Bool +isModel assignment (affine_expr, rel) = + case (eval assignment affine_expr, rel) of + (Just value, LessThan ) -> value < 0 + (Just value, LessEquals ) -> value <= 0 + (Just value, Equals ) -> value == 0 + (Just value, GreaterEquals) -> value >= 0 + (Just value, GreaterThan ) -> value > 0 + (Nothing , _ ) -> False diff --git a/src/Theory/LinearArithmatic/FourierMotzkin.hs b/src/Theory/LinearArithmatic/FourierMotzkin.hs index ba1b236..cfcf9d8 100644 --- a/src/Theory/LinearArithmatic/FourierMotzkin.hs +++ b/src/Theory/LinearArithmatic/FourierMotzkin.hs @@ -1,3 +1,8 @@ +{-| + Fourier-Motzkin is a sound and complete method for solving linear + constraints. The method is simpler but less efficient then Simplex. + We use it only for testing the Simplex implementation. +-} module Theory.LinearArithmatic.FourierMotzkin (fourierMotzkin) where import Theory.LinearArithmatic.Constraint @@ -5,98 +10,127 @@ import qualified Data.IntMap.Strict as M import qualified Data.IntSet as S import Control.Monad (guard) import qualified Data.Vector as V -import Data.Maybe (fromMaybe, mapMaybe) +import Data.Maybe (fromMaybe, mapMaybe, maybeToList, listToMaybe, fromJust) import Data.List (partition) import Debug.Trace +import Control.Arrow (Arrow (first)) +import Control.Exception (assert) +import Data.Either (partitionEithers) +import Data.Foldable (asum) +import Control.Applicative ((<|>)) +import Utils (fixpoint) + +partitionByBound :: [Constraint] -> Var -> ([Constraint], [Constraint], [Constraint]) +partitionByBound constraints var = + let + go :: Constraint -> ([Constraint], [Constraint], [Constraint]) + go constraint = + case constraint `solveFor` var of + Nothing -> -- constraint does not include var + ([], [], [constraint]) + Just c@(expr, GreaterEquals) -> -- lower bound + ([c], [], []) + Just c@(expr, LessEquals) -> -- upper bound + ([], [c], []) + Just (expr, not_supported_yet) -> + error "TODO: support all constraint types in Fourier Motzkin" + + combine (as1, as2, as3) (bs1, bs2, bs3) = + (as1 ++ bs1, as2 ++ bs2, as3 ++ bs3) + in + foldr combine ([], [], []) $ go <$> constraints {-| - Trivial constraints have no variables and are satisfied, such as 0 < 1. + Identifies a variable that has both upper- and lower bounds, if one exists, as in + + x <= 2y - 1 (upper bound) + 3 <= x (lower bound) + + Then constructs new constraints, where the variable is eliminated, by pairing + each upper- with each lower bound: + + 3 <= x <= 2y - 1 ==> 3 <= 2y - 1 + + and rearranging to make right-hand-side zero. + + -2y + 2 <= 0 + -} -isTrivial :: Constraint -> Bool -isTrivial (rel, AffineExpr constant coeff_map) = - S.null (M.keysSet coeff_map) && case rel of - LessThan -> constant < 0 - LessEquals -> constant <= 0 - Equals -> constant == 0 - GreaterEquals -> constant >= 0 - GreaterThan -> constant > 0 - -eliminate :: Var -> [Constraint] -> [Constraint] -eliminate var constraints = - let - (constraints_with_var, constraints_without_var) = - partition (var `appearsIn`) constraints - - constraints_with_var_eliminated :: [Constraint] - constraints_with_var_eliminated = do - -- c1: x - 1 >= 0 - -- c2: 2x + 4 <= 0 - let constraints_with_var' = - mapMaybe (`solveFor` var) constraints_with_var - - -- c1: x <= 1 - (upper_bound_rel, AffineExpr ub_const ub_coeffs) <- constraints_with_var' - guard (upper_bound_rel == LessEquals) - - -- c2: x >= -2 - (lower_bound_rel, AffineExpr lb_const lb_coeffs) <- constraints_with_var' - guard (lower_bound_rel == GreaterEquals) - - -- -2 <= x <= 1 - -- -2 <= 1 - -- -2 - 1 <= 0 - return - ( LessEquals - , AffineExpr - (lb_const - ub_const) - (M.unionWith (-) lb_coeffs ub_coeffs) - ) - in - constraints_without_var <> constraints_with_var_eliminated +eliminate :: [Constraint] -> Maybe (Var, [Constraint]) +eliminate constraints = listToMaybe $ do + var <- S.toList $ S.unions $ varsIn <$> constraints -fourierMotzkin :: [Constraint] -> Bool -fourierMotzkin = go . fmap reject_zero_coeffs + let (lower_bounds, upper_bounds, constraints_without_var) = + partitionByBound constraints var + + guard $ not (null lower_bounds) && not (null upper_bounds) + + let constraints_with_var_eliminated :: [Constraint] + constraints_with_var_eliminated = do + AffineExpr ub_const ub_coeffs <- fst <$> upper_bounds + AffineExpr lb_const lb_coeffs <- fst <$> lower_bounds + let left_hand_side = AffineExpr (lb_const - ub_const) (M.unionWith (+) lb_coeffs $ negate <$> ub_coeffs) + return (left_hand_side, LessEquals) + + return (var, constraints_without_var ++ constraints_with_var_eliminated) + +fourierMotzkin :: [Constraint] -> Maybe Assignment +fourierMotzkin = go . fmap normalize where - reject_zero_coeffs :: Constraint -> Constraint - reject_zero_coeffs (rel, AffineExpr constant expr) = - (rel, AffineExpr constant $ M.filter (/= 0) expr) - - go :: [Constraint] -> Bool - go constraints - | not (all isTrivial constraints_no_vars) = False - | null constraints_with_vars = True - | otherwise = fourierMotzkin (eliminate next_var constraints_with_vars) + go :: [Constraint] -> Maybe Assignment + go constraints = do + let (constraints_no_vars, constraints_with_vars) = + partition (S.null . varsIn) constraints + + -- otherwise `constraints_no_vars` contains a constraint like 1 <= 0 + guard $ all (M.empty `isModel`) constraints_no_vars + + if null constraints_with_vars then + Just M.empty + else + case eliminate constraints_with_vars of + Nothing -> do + let unassigned_vars = S.unions $ varsIn <$> constraints_with_vars + initial_assignment = M.fromSet (const 0) unassigned_vars + return $ fixpoint (`refine` constraints_with_vars) initial_assignment + + Just (next_var, constraints_without_var) -> do + partial_assignment <- go constraints_without_var + let constraints' = first (substitute partial_assignment) <$> constraints_with_vars + return $ refine (M.insert next_var 0 partial_assignment) constraints' + + refine_with :: Constraint -> Var -> Assignment -> Assignment + refine_with (expr, rel) var assignment = + let + current_value = assignment M.! var + assignment' = M.delete var assignment + in + case (substitute assignment' expr, rel) `solveFor` var of + Nothing -> assignment + Just (AffineExpr new_value _, LessEquals) -> + if new_value < current_value then + M.insert var new_value assignment + else + assignment + Just (AffineExpr new_value _, GreaterEquals) -> + if new_value > current_value then + M.insert var new_value assignment + else + assignment + + refine :: Assignment -> [Constraint] -> Assignment + refine = foldr accum where - (constraints_no_vars, constraints_with_vars) = - partition (S.null . varsIn) constraints - - next_var = fst - $ M.findMin - $ getCoeffMap . snd - $ head constraints_with_vars - --- linearDependent :: Constraint -> Constraint -> Bool --- linearDependent vs ws = coeff /= 0 && V.map (* coeff) vs == ws --- where --- div' x y --- | y == 0 = 0 --- | otherwise = x/y --- coeff = --- fromMaybe 0 $ V.find (/=0) $ V.zipWith div' ws vs - --------------------------------------- - --- SAT -example1 = - [ ( GreaterEquals - , AffineExpr 1 (M.singleton 0 (-1)) - ) - , ( GreaterEquals - , AffineExpr 0 (M.singleton 0 7) - ) - ] + accum :: Constraint -> Assignment -> Assignment + accum constraint assignment = + foldr (refine_with constraint) assignment (S.toList $ varsIn constraint) + +--------------------------- + +example1 = + [ (AffineExpr 0 $ M.fromList [ (0, -1), (1, -1) ], LessEquals) + , (AffineExpr 1 $ M.fromList [ (0, 1) ], LessEquals) ] example2 = - [ (LessEquals, AffineExpr 0 $ M.singleton 0 (-1)) - , (GreaterEquals, AffineExpr 1 $ M.singleton 0 (-1)) - ] + [ (AffineExpr 0 $ M.singleton 0 (-1), LessEquals) + , (AffineExpr 1 $ M.singleton 0 1, LessEquals) ] diff --git a/src/Theory/LinearArithmatic/Simplex.hs b/src/Theory/LinearArithmatic/Simplex.hs index a5677ca..d1aecaa 100644 --- a/src/Theory/LinearArithmatic/Simplex.hs +++ b/src/Theory/LinearArithmatic/Simplex.hs @@ -35,7 +35,7 @@ initTableau constraints = fresh_vars = [max_original_var + 1 ..] (basis, bounds) = bimap M.fromList M.fromList $ unzip $ do - (slack_var, (relation, affine_expr)) <- zip fresh_vars constraints + (slack_var, (affine_expr, relation)) <- zip fresh_vars constraints let bound_type = case relation of @@ -152,8 +152,8 @@ rewrite (x, expr_x) (y, expr_y) = -} solveFor' :: Equation -> Var -> Maybe Equation solveFor' (x, expr_x) y = do - let constraint = (Equals, AffineExpr 0 $ M.insert x (-1) expr_x) - (_, AffineExpr _ expr_y) <- solveFor constraint y + let constraint = (AffineExpr 0 $ M.insert x (-1) expr_x, Equals) + (AffineExpr _ expr_y, _) <- solveFor constraint y return (y, expr_y) {-| @@ -193,10 +193,14 @@ pivot basic_var non_basic_var (Tableau basis bounds assignment) = new_value_non_basic_var = old_value_non_basic_var + (old_value_basic_var - new_value_basic_var) / non_basic_var_coeff assignment' = - M.insert non_basic_var new_value_non_basic_var $ - M.insert basic_var new_value_basic_var assignment - - assignment'' = M.union (eval assignment' . AffineExpr 0 <$> basis') assignment' + M.insert non_basic_var new_value_non_basic_var + $ M.insert basic_var new_value_basic_var assignment + + assignment'' = + M.union + ( from_just "Can't evaluate expression, because assignment is partial" + $ traverse (eval assignment' . AffineExpr 0) basis' ) + assignment' in Tableau basis' bounds assignment'' @@ -253,19 +257,19 @@ simplex constraints = do -- SAT example1 = simplex - [ (LessEquals, AffineExpr (-3) $ M.fromList [(0, 1), (1, 1)]) - , (GreaterEquals, AffineExpr (-1) $ M.fromList [(0, 1), (1, 1)]) - , (LessEquals, AffineExpr (-3) $ M.fromList [(0, 1), (1, -1)]) - , (GreaterEquals, AffineExpr (-1) $ M.fromList [(0, 1), (1, -1)]) + [ (AffineExpr (-3) $ M.fromList [(0, 1), (1, 1)], LessEquals) + , (AffineExpr (-1) $ M.fromList [(0, 1), (1, 1)], GreaterEquals) + , (AffineExpr (-3) $ M.fromList [(0, 1), (1, -1)], LessEquals) + , (AffineExpr (-1) $ M.fromList [(0, 1), (1, -1)], GreaterEquals) ] -- SAT example2 = - [ ( GreaterEquals - , AffineExpr + [ ( AffineExpr { getConstant = -100 , getCoeffMap = M.fromList [ ( 0 , -100) ] } + , GreaterEquals ) ] diff --git a/src/Utils.hs b/src/Utils.hs index 2fafa3e..bf5526a 100644 --- a/src/Utils.hs +++ b/src/Utils.hs @@ -3,6 +3,8 @@ import Data.Maybe (isJust, catMaybes) import Data.Foldable (toList, concatMap) import qualified Data.Set as S import Data.Set (Set) +import Control.Exception (assert) +import Data.List (uncons) fixpoint :: Eq a => (a -> a) -> a -> a fixpoint f a @@ -16,4 +18,4 @@ takeWhileJust :: [Maybe a] -> [a] takeWhileJust = catMaybes . takeWhile isJust distinct :: Ord a => [a] -> [a] -distinct = toList . S.fromList \ No newline at end of file +distinct = toList . S.fromList diff --git a/test/Test.hs b/test/Test.hs index 2ea072d..9b6c98c 100644 --- a/test/Test.hs +++ b/test/Test.hs @@ -21,8 +21,8 @@ main = defaultMain $ checkParallel <$> , ("No roots are lost", prop_no_roots_are_lost) ] , Group "Linear Arithmatic" - [ ("Simplex is sound", withTests 1000 $ prop_simplex_sound) - -- , ("Fourier-Motzkin is sound" prop_fourier_motzkin_sound) + [ ("Simplex is sound", withTests 100 $ prop_simplex_sound) + , ("Fourier-Motzkin is sound", prop_fourier_motzkin_sound) , ("Fourier-Motzkin equivalent to Simplex", prop_fourier_motzkin_equiv_simplex) ] ] diff --git a/test/TestLinearArithmatic.hs b/test/TestLinearArithmatic.hs index 5721aae..10ebdee 100644 --- a/test/TestLinearArithmatic.hs +++ b/test/TestLinearArithmatic.hs @@ -1,6 +1,6 @@ -{-# LANGUAGE ScopedTypeVariables #-} module TestLinearArithmatic ( prop_simplex_sound + , prop_fourier_motzkin_sound , prop_fourier_motzkin_equiv_simplex ) where @@ -14,25 +14,7 @@ import Theory.LinearArithmatic.Constraint import Theory.LinearArithmatic.FourierMotzkin (fourierMotzkin) import Data.Maybe (isJust) -isModel :: Assignment -> Constraint -> Bool -isModel assignment (rel, affine_expr) = - case rel of - LessThan -> eval assignment affine_expr < 0 - LessEquals -> eval assignment affine_expr <= 0 - Equals -> eval assignment affine_expr == 0 - GreaterEquals -> eval assignment affine_expr >= 0 - GreaterThan -> eval assignment affine_expr > 0 - -- TODO: generate more representative constraint sets --- newtype Constraint' a = Constraint' (Constraint a) - --- instance Arbitrary a => Arbitrary (Constraint' a) where --- arbitrary = do --- var <- chooseInt (0, 10) --- coeff <- arbitrary --- _ --- return $ Constraint' - genConstraints :: Int -> Gen [Constraint] genConstraints max_constraints = Gen.list (Range.linear 1 max_constraints) $ do linear_expr <- fmap M.fromList $ Gen.list (Range.linear 0 10) $ do @@ -47,7 +29,8 @@ genConstraints max_constraints = Gen.list (Range.linear 1 max_constraints) $ do constant <- toRational <$> Gen.int (Range.linearFrom 0 (-100) 100) -- constant <- Gen.float (Range.linearFrac (-100.0) 100.0) - return (rel, AffineExpr constant linear_expr) + -- TODO: make sure that constraints are always normalized by construction. + return $ normalize (AffineExpr constant linear_expr, rel) prop_simplex_sound :: Property prop_simplex_sound = property $ do @@ -64,12 +47,17 @@ prop_simplex_sound = property $ do -- TODO: -- invariant: assignment always satisfies `A*x = 0` --- TODO --- prop_fourier_motzkin_sound :: Property --- prop_fourier_motzkin_sound = _ +prop_fourier_motzkin_sound :: Property +prop_fourier_motzkin_sound = property $ do + constraints <- forAll (genConstraints 10) + case fourierMotzkin constraints of + Nothing -> success + Just assignment -> do + annotate $ show assignment + assert $ all (assignment `isModel`) constraints -- TODO: gets "stuck" on some inputs. Too slow? Memory leak? prop_fourier_motzkin_equiv_simplex :: Property prop_fourier_motzkin_equiv_simplex = property $ do constraints <- forAll (genConstraints 10) - fourierMotzkin constraints === isJust (simplex constraints) + isJust (fourierMotzkin constraints) === isJust (simplex constraints)