Simplex: equality constraints and cutting planes

Extend simplex with equality constraint support and improve API
to allow easy incremental adding of new constraints (cutting planes).

src/Theory/LinearArithmatic/BranchAndBound.hs

@@ -1,104 +1,116 @@
-module Theory.LinearArithmatic.BranchAndBound 
+module Theory.LinearArithmatic.BranchAndBound
   ( branchAndBound
   , isIntegral
+  , steps
+  , isSat
   ) where
 
-import Theory.LinearArithmatic.Simplex ( simplexSteps, initTableau, Tableau(..), BoundType(..), isSatisfied )
+import qualified Theory.LinearArithmatic.Simplex as Simplex
+import Theory.LinearArithmatic.Simplex (Tableau, BoundType(..))
 import Theory.LinearArithmatic.Constraint
-import Data.Ratio (denominator, numerator)
+import Data.Ratio (denominator, numerator, (%))
 import qualified Data.IntSet as S
 import qualified Data.IntMap as M
 import Control.Applicative ((<|>))
 import Data.Maybe (fromJust)
+import Utils (takeUntil, firstJust)
 
 isIntegral :: Rational -> Bool
-isIntegral rational = 
+isIntegral rational =
   numerator rational `mod` denominator rational == 0
 
-findFractional :: S.IntSet -> BBState -> Maybe (Var, Rational)
-findFractional vars state =  M.lookupMin 
-  $ M.filter (not . isIntegral)
-  $ M.restrictKeys (getAssignment $ getTableau state) vars
-
-data BBState = BB 
-  { getFreshVar :: Var
-  , getTableau :: Tableau
-  } deriving Show
-
-initState :: [Constraint] -> Maybe BBState
-initState constraints = do 
-  tableau_0 <- initTableau constraints 
-
-  let fresh_var =
-        case M.lookupMax $ getAssignment tableau_0 of
-          Nothing           -> 0
-          Just (max_var, _) -> max_var + 1
-
-  return $ BB fresh_var tableau_0
-
-branch :: Var -> Rational -> BoundType -> BBState -> BBState
-branch var frac_value bound_type (BB fresh_var tableau) = 
-  let Tableau basis bounds assignment = tableau
-
-      int_value = 
-        case bound_type of
-          UpperBound -> fromIntegral $ floor frac_value
-          LowerBound -> fromIntegral $ ceiling frac_value
-
-      tableau_row =
-        case M.lookup var basis of
-          Just expr -> expr
-          Nothing   -> M.singleton var 1
-
-      fresh_var_value = fromJust $ eval assignment $ AffineExpr 0 tableau_row
-
-   in BB (fresh_var + 1) $ Tableau
-      (M.insert fresh_var tableau_row basis)
-      (M.insert fresh_var (bound_type, int_value) bounds)
-      (M.insert fresh_var fresh_var_value assignment)
-
-solveRelaxation :: BBState -> Maybe BBState
-solveRelaxation (BB fresh_var tableau)
-  | isSatisfied tableau' = Just $ BB fresh_var tableau'
+findFractional :: S.IntSet -> Tableau -> Maybe (Var, Rational)
+findFractional vars tableau =
+  M.lookupMin
+    $ M.filter (not . isIntegral)
+    $ M.restrictKeys (Simplex.getAssignment tableau) vars
+
+branchOn :: Var -> Rational -> Tableau -> (Tableau, Tableau)
+branchOn var frac_value tableau =
+  let
+    left_branch :: Tableau
+    left_branch = 
+      Simplex.addTableauRow
+        (AffineExpr (- fromIntegral (floor frac_value)) (M.singleton var 1))
+        UpperBound
+        tableau
+
+    right_branch :: Tableau
+    right_branch = 
+      Simplex.addTableauRow
+        (AffineExpr (- fromIntegral (ceiling frac_value)) (M.singleton var 1))
+        LowerBound
+        tableau
+  in
+    (left_branch, right_branch)
+
+solveRelaxation :: Tableau -> Maybe Tableau
+solveRelaxation tableau
+  | Simplex.isSatisfied tableau' = Just tableau'
   | otherwise = Nothing
   where
-    tableau' = last $ simplexSteps tableau 
-
+    tableau' = last $ Simplex.steps tableau
+
+data BBState = Sat Assignment | Pending Tableau
+
+steps :: [Constraint] -> S.IntSet -> [BBState] 
+steps constraints int_vars = 
+  let 
+    original_vars = varsInAll constraints
+
+    get_model :: Tableau -> Assignment
+    get_model tableau = M.restrictKeys (Simplex.getAssignment tableau) original_vars
+
+    go :: BBState -> [BBState]
+    go (Sat _)           = []
+    go (Pending tableau) = Pending tableau : 
+      case solveRelaxation tableau of
+        Nothing -> []
+        Just tableau_next -> 
+          case findFractional int_vars tableau_next of
+            Nothing                -> [ Sat (get_model tableau_next) ]
+            Just (var, frac_value) -> 
+              let 
+                (left_branch, right_branch) = branchOn var frac_value tableau_next
+              in
+                go (Pending left_branch) <|> go (Pending right_branch)
+  in
+    maybe [] (go . Pending) $ Simplex.initTableau constraints
+
+isSat :: BBState -> Maybe Assignment
+isSat (Sat model) = Just model
+isSat _           = Nothing
+
+{-|
+  TODO: 
+    * Currently incomplete, i.e. algorithm may not terminate on some inputs.
+      - especially equality constraints
+    * Constraint preprocessing: 
+      - remove redundant constraints
+      - tighten bounds
+
+-}
 branchAndBound :: [Constraint] -> S.IntSet -> Maybe Assignment
-branchAndBound constraints int_vars = do
-  let original_vars = varsInAll constraints
-
-      restrict :: Assignment -> Assignment
-      restrict assignment = M.restrictKeys assignment original_vars
-
-      go :: BBState -> Maybe BBState
-      go state = do
-        state' <- solveRelaxation state
-        case findFractional int_vars state' of
-          Nothing                -> Just state'
-          Just (var, frac_value) -> go left_branch <|> go right_branch
-            where
-              left_branch = branch var frac_value UpperBound state'
-              right_branch = branch var frac_value LowerBound state'
-
-  state_0 <- initState constraints 
-  state_n <- go state_0
-
-  case findFractional int_vars state_n of
-    Nothing -> Just $ restrict $ getAssignment $ getTableau state_n
-    Just _  -> Nothing
+branchAndBound constraints int_vars = 
+  firstJust isSat $ steps constraints int_vars
 
 --------------------
 
-example = 
-  [ (AffineExpr (-2)   $ M.fromList [ (0,1), (1,1) ], GreaterEquals )
-  , (AffineExpr (-1/2) $ M.fromList [ (0, 1), (1, -1)], LessEquals )
-  , (AffineExpr (-1) $ M.fromList [ (0, 1) ], LessEquals )
-  , (AffineExpr 1      $ M.fromList [ (0, 1), (1, -1)], GreaterEquals )
-  , (AffineExpr (-3/2) $ M.fromList [ (0,1) ], LessEquals )
-  , (AffineExpr (-3/2) $ M.fromList [ (1,1) ], GreaterEquals )
-  , (AffineExpr (-7/4) $ M.fromList [ (1,1) ], LessEquals )
+example =
+  [ (AffineExpr (-2)   $ M.fromList [ (0,1), (1,1) ]   , GreaterEquals )
+  , (AffineExpr (-1/2) $ M.fromList [ (0, 1), (1, -1)] , LessEquals    )
+  , (AffineExpr (-1)   $ M.fromList [ (0, 1) ]         , LessEquals    )
+  , (AffineExpr 1      $ M.fromList [ (0, 1), (1, -1)] , GreaterEquals )
+  , (AffineExpr (-3/2) $ M.fromList [ (0,1) ]          , LessEquals    )
+  , (AffineExpr (-3/2) $ M.fromList [ (1,1) ]          , GreaterEquals )
+  , (AffineExpr (-7/4) $ M.fromList [ (1,1) ]          , LessEquals    )
   ]
 
-example2 = 
-  [ (AffineExpr 1 $ M.fromList [ (0,-2), (1,1) ], LessEquals ) ]
+example2 =
+  [ (AffineExpr 1 $ M.fromList [ (0,-2), (1,1) ], LessEquals ) ]
+
+example3 = 
+  [ ( AffineExpr (-1) (M.fromList [(0,3),(1,-3)]), GreaterEquals )
+  , ( AffineExpr 0 (M.singleton 0 1), LessEquals)
+  , ( AffineExpr 1 (M.singleton 1 1), LessEquals)
+  ]

src/Theory/LinearArithmatic/Constraint.hs

@@ -57,6 +57,26 @@ instance Num AffineExpr where
   negate (AffineExpr constant coeffs) = 
     AffineExpr (negate constant) (M.map negate coeffs)
 
+-- | Construct `Equals` constraint.
+(.==) :: AffineExpr -> AffineExpr -> Constraint
+(.==) lhs_expr rhs_expr = (lhs_expr - rhs_expr, Equals)
+
+-- | Construct `LessEquals` constraint.
+(.<=) :: AffineExpr -> AffineExpr -> Constraint
+(.<=) lhs_expr rhs_expr = (lhs_expr - rhs_expr, LessEquals)
+
+-- | Construct `GreaterEquals` constraint.
+(.>=) :: AffineExpr -> AffineExpr -> Constraint
+(.>=) lhs_expr rhs_expr = (lhs_expr - rhs_expr, GreaterEquals)
+
+-- | Construct `GreaterThan` constraint.
+(.>) :: AffineExpr -> AffineExpr -> Constraint
+(.>) lhs_expr rhs_expr = (lhs_expr - rhs_expr, GreaterThan)
+
+-- | Construct `LessThan` constraint.
+(.<) :: AffineExpr -> AffineExpr -> Constraint
+(.<) lhs_expr rhs_expr = (lhs_expr - rhs_expr, LessThan)
+
 instance Show AffineExpr where
   show (AffineExpr constant coeff_map) =
     let 
@@ -74,10 +94,17 @@ instance Show AffineExpr where
           show_ratio ratio
 
       show_var (var, coeff) = 
-        show_signed coeff ++ "*x" ++ show var
+        if coeff == 1 then
+          "x" ++ show var
+        else 
+          show_signed coeff ++ "*x" ++ show var
 
+      terms_showed :: [String]
+      terms_showed = fmap show_var (M.toList coeff_map) ++ [ show_signed constant | constant /= 0 ]
      in 
-      intercalate " + " (fmap show_var (M.toList coeff_map) ++ [ show_signed constant | constant /= 0 ])
+      case terms_showed of
+        [] -> "0"
+        _  -> intercalate "+" terms_showed
 
 varsIn :: Constraint -> S.IntSet
 varsIn (AffineExpr _ coeff_map, _) = M.keysSet coeff_map

src/Theory/LinearArithmatic/FourierMotzkin.hs

@@ -38,7 +38,7 @@ partitionByBound constraints var =
     combine (as1, as2, as3) (bs1, bs2, bs3) =
       (as1 ++ bs1, as2 ++ bs2, as3 ++ bs3)
   in
-    foldr combine ([], [], []) $ go <$> constraints
+    foldr (combine . go) ([], [], []) constraints
 
 {-|
   Identifies a variable that has both upper- and lower bounds, if one exists, as in
@@ -70,12 +70,18 @@ eliminate constraints = listToMaybe $ do
         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 $ normalize $ (left_hand_side, LessEquals)
+        return $ normalize (left_hand_side, LessEquals)
 
   return (var, constraints_without_var ++ constraints_with_var_eliminated)
 
+-- | TODO: a bit ad-hoc way to deal with equlities by just turing them into two inequalities.
+-- I feel like there is a better way.
+replaceEqualities :: Constraint -> [Constraint] 
+replaceEqualities (expr, Equals) = [ (expr, LessEquals), (expr, GreaterEquals) ]
+replaceEqualities constraint     = [ constraint ]
+
 fourierMotzkin :: [Constraint] -> Maybe Assignment
-fourierMotzkin = go . fmap normalize
+fourierMotzkin = go . concatMap (replaceEqualities . normalize)
   where
     go :: [Constraint] -> Maybe Assignment
     go constraints = do