diff --git a/SMT.cabal b/SMT.cabal index 406e596..e1a3f0d 100644 --- a/SMT.cabal +++ b/SMT.cabal @@ -17,6 +17,7 @@ library , Theory.BitVectors , Theory.UninterpretedFunctions.Lazy , Theory.UninterpretedFunctions.Eager + , Theory.LinearArithmatic.Constraint , Theory.LinearArithmatic.FourierMotzkin , Theory.LinearArithmatic.Simplex , Theory.LinearArithmatic.BranchAndBound diff --git a/src/Theory/LinearArithmatic/Constraint.hs b/src/Theory/LinearArithmatic/Constraint.hs new file mode 100644 index 0000000..8c1a05b --- /dev/null +++ b/src/Theory/LinearArithmatic/Constraint.hs @@ -0,0 +1,106 @@ +module Theory.LinearArithmatic.Constraint where + +import qualified Data.IntMap.Strict as M +import qualified Data.IntSet as S +import Control.Monad (guard) +import qualified Data.IntMap.Merge.Lazy as MM +import Data.List (intercalate) +import Data.Ratio (denominator, numerator) + +type Var = Int + +data ConstraintRelation = LessThan | LessEquals | Equals | GreaterEquals | GreaterThan + deriving (Show, Eq) + +flipRelation :: ConstraintRelation -> ConstraintRelation +flipRelation = \case + LessThan -> GreaterThan + LessEquals -> GreaterEquals + Equals -> Equals + GreaterEquals -> LessEquals + GreaterThan -> LessThan + +-- | Map from variable IDs to assigned values +type Assignment = M.IntMap Rational + +-- | +-- For example, constraint such as `3x + y + 3 <= 0` is represented as: +-- +-- (LessEquals, AffineExpr 3 (3x+y)) +-- +type Constraint = (ConstraintRelation, AffineExpr) + +-- | Map from variable IDs to coefficients +type LinearExpr = M.IntMap Rational + +data AffineExpr = AffineExpr + { getConstant :: Rational + , getCoeffMap :: LinearExpr } + +instance Show AffineExpr where + show (AffineExpr constant coeff_map) = + let + show_ratio ratio = + show (numerator ratio) ++ + if denominator ratio == 1 then + "" + else + "/" ++ show (denominator ratio) + + show_signed ratio = + if ratio < 0 then + "(" ++ show_ratio ratio ++ ")" + else + show_ratio ratio + + show_var (var, coeff) = + show_signed coeff ++ "*x" ++ show var + + in + intercalate " + " (fmap show_var (M.toList coeff_map) ++ [show_signed constant]) + +varsIn :: Constraint -> S.IntSet +varsIn (_, AffineExpr _ coeff_map) = M.keysSet coeff_map + +appearsIn :: Var -> Constraint -> Bool +appearsIn var = M.member var . getCoeffMap . snd + +{-| + Solve a constraint for a variable. For example solving + + 0 <= 3x + 10y - 1 + + for x, yields + + x >= -10/3y + 1/3 + + Returns `Nothing`, if the given variable does not appear in the + expression or the coefficient is zero. +-} +solveFor :: Constraint -> Var -> Maybe Constraint +solveFor (rel, AffineExpr constant coeffs) x = do + coeff_of_x <- M.lookup x coeffs + guard (coeff_of_x /= 0) + + let rel' = if coeff_of_x < 0 then flipRelation rel else rel + + return + -- flip the relation: x >= -10/3y + 1/3 + $ (rel', ) + -- also divide constant term: x <= -10/3y + 1/3 + $ AffineExpr (- constant / coeff_of_x) + -- divide coefficients: x <= -10/3y - 1 + $ M.map (/ (-coeff_of_x)) + -- 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 + diff --git a/src/Theory/LinearArithmatic/FourierMotzkin.hs b/src/Theory/LinearArithmatic/FourierMotzkin.hs index e81355e..ba1b236 100644 --- a/src/Theory/LinearArithmatic/FourierMotzkin.hs +++ b/src/Theory/LinearArithmatic/FourierMotzkin.hs @@ -1,111 +1,102 @@ -{-# LANGUAGE OverloadedLists #-} -module Theory.LinearArithmatic.FourierMotzkin where +module Theory.LinearArithmatic.FourierMotzkin (fourierMotzkin) where -import Data.Ratio (Ratio) -import qualified Data.Vector as V -import Data.Maybe (fromMaybe) +import Theory.LinearArithmatic.Constraint +import qualified Data.IntMap.Strict as M +import qualified Data.IntSet as S import Control.Monad (guard) -import Data.List ((\\)) - --- [a,b,c,d] :<=> (a*x + b*y + c*z + d <= 0) -type Constraint = V.Vector (Ratio Integer) - -constraints :: [Constraint] -constraints = - [ [-2, 4,-5, 4 ] - , [-4, 5, 2, 1 ] - , [-5,-4, 0, 3 ] - , [ 3, 2,-4, 2 ] - , [ 1, 3, 4,-4 ] - , [ 4, 1, 1,-5 ] ] - -data Bound = Upper Constraint | Lower Constraint | UnBounded - deriving Show - -solveFor :: Int -> Constraint -> Bound -solveFor i vs - | x > 0 = Upper $ V.imap (\j y -> if j == i then 0 else y/x) vs - | x < 0 = Lower $ V.imap (\j y -> if j == i then 0 else -y/x) vs - | otherwise = UnBounded - where - x = vs V.! i - -partitionByBound :: [Bound] -> ([Constraint], [Constraint]) -partitionByBound [] = ([], []) -partitionByBound (b:bs) = - case b of - Lower l -> (l:ls, us) - where (ls, us) = partitionByBound bs - Upper u -> (ls, u:us) - where (ls, us) = partitionByBound bs - UnBounded -> partitionByBound bs - -hasNonZeroVars :: Constraint -> Bool -hasNonZeroVars vs = foldl go 0 vs > 1 - where - go sum x - | x == 0 = sum - | otherwise = sum + 1 - -fourierMotzkin :: [Constraint] -> [Constraint] -fourierMotzkin = go 0 - where - go :: Int -> [Constraint] -> [Constraint] - go i constraints = - if any hasNonZeroVars constraints then - let bounds = solveFor i <$> constraints - (lowers, uppers) = partitionByBound bounds - new_constraints = [ V.zipWith (+) l u | l <- lowers, u <- uppers ] - in new_constraints <> go (i+1) new_constraints - else - [] - -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 - -f :: [Constraint] -f = do - let cs' = filter hasNonZeroVars $ fourierMotzkin constraints - cs = cs' -- <> constraints - c1 <- cs - guard (not (any (linearDependent c1) (cs \\ [c1]))) - -- return (V.map realToFrac c1) -- , V.map realToFrac c2] - return c1 -- , V.map realToFrac c2] - ----------------------------------------------------- - --- simplexStep :: - -{- - -data Term = - Var (Ratio Integer) Int -- variable with rational coefficient - | Const (Ratio Integer) -- constant rational value - | Add Term Term -- sum of terms - -data Atom = Leq Term Term - -type Constraint' = V.Vector (Ratio Integer) - -variables :: Term -> [Int] -variables (Var _ v) = [v] -variables (Const _) = [] -variables (Add t1 t2) = - variables t1 <> variables t2 - -toConstraint :: Atom -> Constraint' -toConstraint (Leq t1 t2) = V.generate vector_size go - where - vector_size = 1 + maximum (variables t1 <> variables t2) - - go :: Int -> Ratio Integer - go i = undefined +import qualified Data.Vector as V +import Data.Maybe (fromMaybe, mapMaybe) +import Data.List (partition) +import Debug.Trace +{-| + Trivial constraints have no variables and are satisfied, such as 0 < 1. -} +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 + +fourierMotzkin :: [Constraint] -> Bool +fourierMotzkin = go . fmap reject_zero_coeffs + 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) + 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) + ) + ] + +example2 = + [ (LessEquals, AffineExpr 0 $ M.singleton 0 (-1)) + , (GreaterEquals, AffineExpr 1 $ M.singleton 0 (-1)) + ] diff --git a/src/Theory/LinearArithmatic/Simplex.hs b/src/Theory/LinearArithmatic/Simplex.hs index 4b8c7b0..a5677ca 100644 --- a/src/Theory/LinearArithmatic/Simplex.hs +++ b/src/Theory/LinearArithmatic/Simplex.hs @@ -3,7 +3,7 @@ to sets of linear constraints. All coefficients, bounds, and solutions are of type `Rational` to avoid having to deal with numeric issues at the cost of performance. -} -module Theory.LinearArithmatic.Simplex (simplex, Assignment, Constraint, ConstraintRelation (..), eval) where +module Theory.LinearArithmatic.Simplex (simplex) where import Control.Monad (guard) import Data.Bifunctor (bimap) @@ -11,61 +11,41 @@ import qualified Data.IntMap.Strict as M import qualified Data.IntMap.Merge.Strict as MM import qualified Data.IntSet as S import Data.Maybe (fromMaybe, mapMaybe) +import Theory.LinearArithmatic.Constraint import Debug.Trace -type Var = Int - -data ConstraintRelation = LessThan | LessEquals | Equals | GreaterEquals | GreaterThan - deriving (Show) - --- | --- For example, constraint such as `3x + y <= 3` is represented as: --- --- (3x+y, LessEquals, 3) -type Constraint = (LinearTerm, ConstraintRelation, Rational) - --- | Map from variable IDs to assigned values -type Assignment = M.IntMap Rational - --- | Map from variable IDs to coefficients -type LinearTerm = M.IntMap Rational - data BoundType = UpperBound | LowerBound deriving (Show) data Tableau = Tableau - { getBasis :: M.IntMap LinearTerm, + { getBasis :: M.IntMap LinearExpr, getBounds :: M.IntMap (BoundType, Rational), getAssignment :: Assignment } deriving (Show) -varsIn :: [Constraint] -> S.IntSet -varsIn constraints = S.unions $ do - (linear_term, _, _) <- constraints - return $ M.keysSet linear_term - --- | --- Transform constraints to "General From" and initialize Simplex tableau. +{-| + Transform constraints to "General From" and initialize Simplex tableau. +-} initTableau :: [Constraint] -> Tableau initTableau constraints = let - original_vars = varsIn constraints + original_vars = foldr (S.union . varsIn) S.empty constraints max_original_var = if S.null original_vars then -1 else S.findMax original_vars fresh_vars = [max_original_var + 1 ..] (basis, bounds) = bimap M.fromList M.fromList $ unzip $ do - (slack_var, (linear_term, relation, bound)) <- zip fresh_vars constraints + (slack_var, (relation, affine_expr)) <- zip fresh_vars constraints let bound_type = case relation of - LessEquals -> UpperBound + LessEquals -> UpperBound GreaterEquals -> LowerBound - other_rels -> error "TODO: extend Simplex to all constraint relations" + other_rels -> error "TODO: extend Simplex to all constraint relations" return - ( (slack_var, linear_term), - (slack_var, (bound_type, bound)) + ( (slack_var, getCoeffMap affine_expr) + , (slack_var, (bound_type, - getConstant affine_expr)) ) slack_vars = M.keys bounds @@ -108,12 +88,12 @@ pivotCandidates :: Tableau -> [(Var, Var)] pivotCandidates tableau@(Tableau basis bounds assignment) = do (basic_var, violation) <- violatedBasicVars tableau - let term = basis M.! basic_var + let expr = basis M.! basic_var non_basis = M.keysSet assignment S.\\ M.keysSet basis non_basic_var <- S.toAscList non_basis - let non_basic_var_coeff = M.findWithDefault 0 non_basic_var term + let non_basic_var_coeff = M.findWithDefault 0 non_basic_var expr non_basic_bound_type = fst <$> M.lookup non_basic_var bounds can_pivot = case (non_basic_bound_type, signum non_basic_var_coeff, violation) of @@ -146,97 +126,79 @@ pivotCandidates tableau@(Tableau basis bounds assignment) = do guard can_pivot return (basic_var, non_basic_var) -type Equation = (Var, LinearTerm) - --- | --- Solve an equation for a given variable. For example solving --- --- y = 3x + 10 --- --- for x, yields --- --- x = 1/3 y - 10/3 -solveFor :: Equation -> Var -> Maybe Equation -solveFor (y, right_hand_side) x = do - coeff_of_x <- M.lookup x right_hand_side - guard (coeff_of_x /= 0) - return $ - (x,) - -- divide by coefficient: x = -1/3y - 10/3 - $ - M.map (/ (-coeff_of_x)) - -- exchange x and y: -3x = -y + 10 - $ - M.insert y (-1) $ - M.delete x right_hand_side - --- | --- Given two equations, such as --- --- e1: x = w + 2z --- e2: y = 3x + 4z --- --- rewrite x in e2 using e1: --- --- y = 3(w + 2z) + 4z --- = 3w + 10z +type Equation = (Var, LinearExpr) + +{-| + Given two equations, such as + + e1: x = w + 2z + e2: y = 3x + 4z + + rewrite x in e2 using e1: + + y = 3(w + 2z) + 4z + = 3w + 10z +-} rewrite :: Equation -> Equation -> Equation -rewrite (x, term_x) (y, term_y) = - let coeff_of_x = M.findWithDefault 0 x term_y - term_x' = M.map (* coeff_of_x) term_x - term_y' = M.unionWith (+) (M.delete x term_y) term_x' - in (y, term_y') - -eval :: Assignment -> LinearTerm -> Rational -eval assignment term = - sum $ - MM.merge - MM.dropMissing - MM.dropMissing - (MM.zipWithMatched (const (*))) - assignment - term - --- | --- TODO: make sure: --- --- - only slack variables are pivoted --- - non-basic variables must violate bound --- - basic variable is suitable for pivoting +rewrite (x, expr_x) (y, expr_y) = + let coeff_of_x = M.findWithDefault 0 x expr_y + expr_x' = M.map (* coeff_of_x) expr_x + expr_y' = M.unionWith (+) (M.delete x expr_y) expr_x' + in (y, expr_y') + +{-| + Specialization of `solveFor` for linear equations, instead of affine constraints, i.e. + no constant terms and no inequalities. +-} +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 + return (y, expr_y) + +{-| + TODO: make sure: + + - only slack variables are pivoted + - non-basic variables must violate bound + - basic variable is suitable for pivoting +-} pivot :: Var -> Var -> Tableau -> Tableau pivot basic_var non_basic_var (Tableau basis bounds assignment) = - let from_just msg (Just a) = a - from_just msg Nothing = error msg + let + from_just msg (Just a) = a + from_just msg Nothing = error msg - term = - from_just "Given variable is not actually in the basis ==> invalid pivot pair" $ - M.lookup basic_var basis + expr = + from_just "Given variable is not actually in the basis ==> invalid pivot pair" $ + M.lookup basic_var basis - equation = - from_just "Can't solve for non-basic variable ==> invalid pivot pair" $ - solveFor (basic_var, term) non_basic_var + (_, expr') = + from_just "Can't solve for non-basic variable ==> invalid pivot pair" $ + solveFor' (basic_var, expr) non_basic_var - basis' = - uncurry M.insert equation $ - M.mapWithKey (\y term_y -> snd $ rewrite equation (y, term_y)) $ - M.delete basic_var basis + basis' = + M.insert non_basic_var expr' + $ M.mapWithKey (\y expr_y -> snd $ rewrite (non_basic_var, expr') (y, expr_y)) + $ M.delete basic_var basis - old_value_basic_var = assignment M.! basic_var - new_value_basic_var = - from_just "Basic variable doesn't have a bound so it's not actually violated" $ - snd <$> M.lookup basic_var bounds + old_value_basic_var = assignment M.! basic_var + new_value_basic_var = + from_just "Basic variable doesn't have a bound so it's not actually violated" $ + snd <$> M.lookup basic_var bounds - non_basic_var_coeff = term M.! non_basic_var + non_basic_var_coeff = expr M.! non_basic_var - old_value_non_basic_var = assignment M.! non_basic_var - new_value_non_basic_var = old_value_non_basic_var + (old_value_basic_var - new_value_basic_var) / non_basic_var_coeff + old_value_non_basic_var = assignment M.! non_basic_var + 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.insert non_basic_var new_value_non_basic_var $ + M.insert basic_var new_value_basic_var assignment - assignment'' = M.union (eval assignment' <$> basis') assignment' - in Tableau basis' bounds assignment'' + assignment'' = M.union (eval assignment' . AffineExpr 0 <$> basis') assignment' + in + Tableau basis' bounds assignment'' {-| Tableau rows with only zero entries correspond to constant constraints @@ -261,9 +223,10 @@ eliminateZeroRows tableau@(Tableau basis bounds assignment) = else Nothing --- | --- TODO: --- in case of UNSAT include explanation, i.e. minimal infeasible subset of constraints. +{-| + TODO: + in case of UNSAT include explanation, i.e. minimal infeasible subset of constraints. +-} simplex :: [Constraint] -> Maybe Assignment simplex constraints = do let go :: Tableau -> Tableau @@ -275,11 +238,9 @@ simplex constraints = do tableau_0 <- eliminateZeroRows $ initTableau constraints - -- traceShowM $ getBasis tableau_0 - let tableau_n = go tableau_0 - original_vars = varsIn constraints + original_vars = foldr (S.union . varsIn) S.empty constraints assignment = M.restrictKeys (getAssignment tableau_n) original_vars if null (violatedBasicVars tableau_n) then @@ -289,10 +250,22 @@ simplex constraints = do ----------------------------------------------------------- +-- SAT example1 = simplex - [ (M.fromList [(0, 1), (1, 1)], LessEquals, 3) - , (M.fromList [(0, 1), (1, 1)], GreaterEquals, 1) - , (M.fromList [(0, 1), (1, -1)], LessEquals, 3) - , (M.fromList [(0, 1), (1, -1)], GreaterEquals, 1) + [ (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)]) ] + +-- SAT +example2 = + [ ( GreaterEquals + , AffineExpr + { getConstant = -100 + , getCoeffMap = M.fromList [ ( 0 , -100) ] + } + ) + ] + diff --git a/test/Test.hs b/test/Test.hs index 853fbce..2ea072d 100644 --- a/test/Test.hs +++ b/test/Test.hs @@ -20,8 +20,10 @@ main = defaultMain $ checkParallel <$> [ ("Intervals never widen", prop_intervals_never_widen) , ("No roots are lost", prop_no_roots_are_lost) ] - , Group "Theory - Linear Arithmatic" - [ ("Simplex method is sound", withTests 1000 $ prop_simplex_sound) + , Group "Linear Arithmatic" + [ ("Simplex is sound", withTests 1000 $ prop_simplex_sound) + -- , ("Fourier-Motzkin is sound" prop_fourier_motzkin_sound) + , ("Fourier-Motzkin equivalent to Simplex", prop_fourier_motzkin_equiv_simplex) ] ] @@ -31,3 +33,4 @@ main = defaultMain $ checkParallel <$> -- main :: IO () -- main = do -- recheckAt (Seed 6028160336680363614 11864191702326251993) "1667:" prop_simplex_sound +-- recheckAt (Seed 7297858649592928895 18067415188796872511) "297:" Simplex is sound diff --git a/test/TestLinearArithmatic.hs b/test/TestLinearArithmatic.hs index 193afb4..5721aae 100644 --- a/test/TestLinearArithmatic.hs +++ b/test/TestLinearArithmatic.hs @@ -1,5 +1,8 @@ {-# LANGUAGE ScopedTypeVariables #-} -module TestLinearArithmatic (prop_simplex_sound) where +module TestLinearArithmatic + ( prop_simplex_sound + , prop_fourier_motzkin_equiv_simplex + ) where import Hedgehog hiding (eval) import qualified Hedgehog.Gen as Gen @@ -7,20 +10,18 @@ import qualified Hedgehog.Range as Range import Theory.LinearArithmatic.Simplex import qualified Data.IntMap as M +import Theory.LinearArithmatic.Constraint +import Theory.LinearArithmatic.FourierMotzkin (fourierMotzkin) +import Data.Maybe (isJust) -isModel :: Assignment -> [Constraint] -> Bool -isModel assignment constraints = - let - check :: Constraint -> Bool - check (linear_term, rel, bound) = - case rel of - LessThan -> eval assignment linear_term < bound - LessEquals -> eval assignment linear_term <= bound - Equals -> eval assignment linear_term == bound - GreaterEquals -> eval assignment linear_term >= bound - GreaterThan -> eval assignment linear_term > bound - in - all check constraints +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) @@ -32,25 +33,43 @@ isModel assignment constraints = -- _ -- return $ Constraint' -prop_simplex_sound :: Property -prop_simplex_sound = property $ do - constraints <- forAll $ Gen.list (Range.linear 1 50) $ do - linear_term <- fmap M.fromList $ Gen.list (Range.linear 0 10) $ do - -- coeff <- Gen.float (Range.linearFrac (-1000.0) 1000.0) - coeff <- toRational <$> Gen.int (Range.linear (-100) 100) - var <- Gen.int (Range.linear 0 20) - return (var, coeff) +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 + -- coeff <- Gen.float (Range.linearFrac (-1000.0) 1000.0) + coeff <- toRational <$> Gen.int (Range.linearFrom 0 (-100) 100) + var <- Gen.int (Range.linear 0 20) + return (var, coeff) - -- TODO: extend Simplex to all constraint relation types - rel <- Gen.element [LessEquals, GreaterEquals] + -- TODO: extend Simplex to all constraint relation types + rel <- Gen.element [LessEquals, GreaterEquals] - constant <- toRational <$> Gen.int (Range.linear (-100) 100) - -- constant <- Gen.float (Range.linearFrac (-100.0) 100.0) + constant <- toRational <$> Gen.int (Range.linearFrom 0 (-100) 100) + -- constant <- Gen.float (Range.linearFrac (-100.0) 100.0) - return (linear_term, rel, constant) + return (rel, AffineExpr constant linear_expr) +prop_simplex_sound :: Property +prop_simplex_sound = property $ do + constraints <- forAll (genConstraints 50) case simplex constraints of Nothing -> success Just assignment -> do annotate $ show assignment - assert $ assignment `isModel` constraints + assert $ all (assignment `isModel`) constraints + +-- TODO +-- invariant: non basic variables always satisfy their bounds + +-- TODO: +-- invariant: assignment always satisfies `A*x = 0` + +-- TODO +-- prop_fourier_motzkin_sound :: Property +-- prop_fourier_motzkin_sound = _ + +-- 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)