diff --git a/src/Theory/LinearArithmatic/Constraint.hs b/src/Theory/LinearArithmatic/Constraint.hs index 7558050..3c08c0a 100644 --- a/src/Theory/LinearArithmatic/Constraint.hs +++ b/src/Theory/LinearArithmatic/Constraint.hs @@ -90,7 +90,7 @@ modifyConstant f (AffineExpr constant coeffs) = AffineExpr (f constant) coeffs Returns `Nothing`, if the given variable does not appear in the expression or the coefficient is zero. -} -solveFor :: Constraint -> Var -> Maybe Constraint +solveFor :: Constraint -> Var -> Maybe (Var, ConstraintRelation, AffineExpr) solveFor (AffineExpr constant coeffs, rel) x = do coeff_of_x <- M.lookup x coeffs guard (coeff_of_x /= 0) @@ -99,7 +99,7 @@ solveFor (AffineExpr constant coeffs, rel) x = do return -- flip the relation: x >= -10/3y + 1/3 - $ (, rel') + $ (x, rel', ) -- also divide constant term: x <= -10/3y + 1/3 $ AffineExpr (- constant / coeff_of_x) -- divide coefficients: x <= -10/3y - 1 @@ -156,3 +156,11 @@ isModel assignment (affine_expr, rel) = (Just value, GreaterEquals) -> value >= 0 (Just value, GreaterThan ) -> value > 0 (Nothing , _ ) -> False + +-- | True iff constraint contains no free variables +isGround :: Constraint -> Bool +isGround = M.null . getCoeffMap . fst + +-- | True iff constraint contains at least one free variable. +isOpen :: Constraint -> Bool +isOpen = not . isGround diff --git a/src/Theory/LinearArithmatic/FourierMotzkin.hs b/src/Theory/LinearArithmatic/FourierMotzkin.hs index cfcf9d8..62384ee 100644 --- a/src/Theory/LinearArithmatic/FourierMotzkin.hs +++ b/src/Theory/LinearArithmatic/FourierMotzkin.hs @@ -20,19 +20,19 @@ import Data.Foldable (asum) import Control.Applicative ((<|>)) import Utils (fixpoint) -partitionByBound :: [Constraint] -> Var -> ([Constraint], [Constraint], [Constraint]) -partitionByBound constraints var = +partitionByBound :: [Constraint] -> Var -> ([AffineExpr], [AffineExpr], [Constraint]) +partitionByBound constraints var = let - go :: Constraint -> ([Constraint], [Constraint], [Constraint]) + go :: Constraint -> ([AffineExpr], [AffineExpr], [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) -> + Just (_, GreaterEquals, lower_bound) -> + ([lower_bound], [], []) + Just (_, LessEquals, upper_bound) -> + ([], [upper_bound], []) + Just not_supported_yet -> error "TODO: support all constraint types in Fourier Motzkin" combine (as1, as2, as3) (bs1, bs2, bs3) = @@ -67,10 +67,10 @@ eliminate constraints = listToMaybe $ do 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 + AffineExpr ub_const ub_coeffs <- upper_bounds + AffineExpr lb_const lb_coeffs <- lower_bounds let left_hand_side = AffineExpr (lb_const - ub_const) (M.unionWith (+) lb_coeffs $ negate <$> ub_coeffs) - return (left_hand_side, LessEquals) + return $ normalize $ (left_hand_side, LessEquals) return (var, constraints_without_var ++ constraints_with_var_eliminated) @@ -79,44 +79,42 @@ fourierMotzkin = go . fmap normalize where go :: [Constraint] -> Maybe Assignment go constraints = do - let (constraints_no_vars, constraints_with_vars) = - partition (S.null . varsIn) constraints + let (constraints_ground, constraints_open) = partition isGround constraints - -- otherwise `constraints_no_vars` contains a constraint like 1 <= 0 - guard $ all (M.empty `isModel`) constraints_no_vars + -- otherwise `constraints_ground` contains a constraint like 1 <= 0 + guard $ all (M.empty `isModel`) constraints_ground - if null constraints_with_vars then + if null constraints_open 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 + case eliminate constraints_open of + Nothing -> + return $ fixpoint (`refine` constraints_open) M.empty Just (next_var, constraints_without_var) -> do partial_assignment <- go constraints_without_var - let constraints' = first (substitute partial_assignment) <$> constraints_with_vars + let constraints' = first (substitute partial_assignment) <$> constraints_open return $ refine (M.insert next_var 0 partial_assignment) constraints' refine_with :: Constraint -> Var -> Assignment -> Assignment - refine_with (expr, rel) var assignment = + refine_with (expr, rel) var assignment = let - current_value = assignment M.! var - assignment' = M.delete var assignment + old_value = M.findWithDefault 0 var assignment + assignment_no_value = M.delete var assignment + assignment_old_value = M.insert var old_value assignment in - case (substitute assignment' expr, rel) `solveFor` var of - Nothing -> assignment - Just (AffineExpr new_value _, LessEquals) -> - if new_value < current_value then + case (substitute assignment_no_value expr, rel) `solveFor` var of + Nothing -> assignment_old_value + Just (_, LessEquals, AffineExpr new_value _) -> + if new_value < old_value then M.insert var new_value assignment else - assignment - Just (AffineExpr new_value _, GreaterEquals) -> - if new_value > current_value then + assignment_old_value + Just (_, GreaterEquals, AffineExpr new_value _) -> + if new_value > old_value then M.insert var new_value assignment else - assignment + assignment_old_value refine :: Assignment -> [Constraint] -> Assignment refine = foldr accum @@ -125,12 +123,11 @@ fourierMotzkin = go . fmap normalize 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 = - [ (AffineExpr 0 $ M.singleton 0 (-1), LessEquals) - , (AffineExpr 1 $ M.singleton 0 1, LessEquals) ] +example = + [ (AffineExpr (-1) $ M.singleton 3 1, GreaterEquals ) + , (AffineExpr 0 $ M.fromList [ (0, 65), (3, -26)], LessEquals ) + , (AffineExpr 0 $ M.singleton 0 (-1), GreaterEquals ) + , (AffineExpr 0 $ M.fromList [ (0, 5), (1, 1), (3, -2) ], GreaterEquals ) + ] diff --git a/src/Theory/LinearArithmatic/Simplex.hs b/src/Theory/LinearArithmatic/Simplex.hs index d1aecaa..8857cf0 100644 --- a/src/Theory/LinearArithmatic/Simplex.hs +++ b/src/Theory/LinearArithmatic/Simplex.hs @@ -153,7 +153,7 @@ rewrite (x, expr_x) (y, expr_y) = solveFor' :: Equation -> Var -> Maybe Equation solveFor' (x, expr_x) y = do let constraint = (AffineExpr 0 $ M.insert x (-1) expr_x, Equals) - (AffineExpr _ expr_y, _) <- solveFor constraint y + (_, _, AffineExpr _ expr_y) <- solveFor constraint y return (y, expr_y) {-| diff --git a/src/Theory/NonLinearRealArithmatic/Subtropical.hs b/src/Theory/NonLinearRealArithmatic/Subtropical.hs index f890e07..1dd85e0 100644 --- a/src/Theory/NonLinearRealArithmatic/Subtropical.hs +++ b/src/Theory/NonLinearRealArithmatic/Subtropical.hs @@ -11,26 +11,27 @@ import Theory.NonLinearRealArithmatic.Expr (Expr (..), Var, BinaryOp (..), subst import Theory.NonLinearRealArithmatic.Polynomial import Control.Arrow ((&&&)) --- | --- The frame of a polynomial is a set of points, obtained from the --- exponents of the individual monomials. E.g. for a polynomial over --- variables x,y like --- --- y + 2xy^3 - 3x^2y^2 - x^3 - 4x^4y^4 --- --- we get the following points --- --- (0,1), (1,3), (2,2), (4,4) --- --- The points are then partitioned by the sign of the coefficient. --- --- pos: (0,1) --- neg: (1,3), (2,2), (4,4) --- --- Computing the frame is the basis for identifiying a term that --- dominates the polynomial for sufficently large variables values. --- That in turn is sufficient to find solutions to inequality --- constraints. +{-| + The frame of a polynomial is a set of points, obtained from the + exponents of the individual monomials. E.g. for a polynomial over + variables x,y like + + y + 2xy^3 - 3x^2y^2 - x^3 - 4x^4y^4 + + we get the following points + + (0,1), (1,3), (2,2), (4,4) + + The points are then partitioned by the sign of the coefficient. + + pos: (0,1) + neg: (1,3), (2,2), (4,4) + + Computing the frame is the basis for identifiying a term that + dominates the polynomial for sufficently large variables values. + That in turn is sufficient to find solutions to inequality + constraints. +-} frame :: (Ord a, Num a) => Polynomial a -> ([Monomial], [Monomial]) frame polynomial = undefined -- TODO where @@ -144,7 +145,6 @@ subtropical (Equals, polynomial) = go :: Assignment a -> Polynomial a -> Maybe (Assignment a) go neg_sol polynomial = intermediateRoot polynomial neg_sol <$> positiveSolution polynomial - in case eval one polynomial `compare` 0 of LT -> go one polynomial diff --git a/test/TestLinearArithmatic.hs b/test/TestLinearArithmatic.hs index 10ebdee..b784fc7 100644 --- a/test/TestLinearArithmatic.hs +++ b/test/TestLinearArithmatic.hs @@ -15,26 +15,24 @@ import Theory.LinearArithmatic.FourierMotzkin (fourierMotzkin) import Data.Maybe (isJust) -- TODO: generate more representative constraint sets -genConstraints :: Int -> Gen [Constraint] -genConstraints max_constraints = Gen.list (Range.linear 1 max_constraints) $ do +genConstraints :: Int -> Int -> Gen [Constraint] +genConstraints max_vars 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) + var <- Gen.int (Range.linear 0 max_vars) return (var, coeff) - -- TODO: extend Simplex to all constraint relation types + -- TODO: extend to all constraint types rel <- Gen.element [LessEquals, GreaterEquals] constant <- toRational <$> Gen.int (Range.linearFrom 0 (-100) 100) - -- constant <- Gen.float (Range.linearFrac (-100.0) 100.0) -- 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 - constraints <- forAll (genConstraints 50) + constraints <- forAll (genConstraints 20 50) case simplex constraints of Nothing -> success Just assignment -> do @@ -49,7 +47,7 @@ prop_simplex_sound = property $ do prop_fourier_motzkin_sound :: Property prop_fourier_motzkin_sound = property $ do - constraints <- forAll (genConstraints 10) + constraints <- forAll (genConstraints 5 5) case fourierMotzkin constraints of Nothing -> success Just assignment -> do @@ -59,5 +57,5 @@ prop_fourier_motzkin_sound = property $ do -- 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) + constraints <- forAll (genConstraints 5 5) isJust (fourierMotzkin constraints) === isJust (simplex constraints)