fix in Fourier-Motzkin

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

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 )
+  ]

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)
 
 {-|

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

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)