Skip to content

Commit

Permalink
Construct model in Fourier-Motzkin
Browse files Browse the repository at this point in the history
  • Loading branch information
gruhn committed Apr 3, 2023
1 parent fe256c5 commit 823bf57
Show file tree
Hide file tree
Showing 7 changed files with 221 additions and 140 deletions.
1 change: 1 addition & 0 deletions SMT.cabal
Original file line number Diff line number Diff line change
Expand Up @@ -59,3 +59,4 @@ test-suite test
hs-source-dirs: test
default-language: Haskell2010
default-extensions: StrictData
, ScopedTypeVariables
82 changes: 67 additions & 15 deletions src/Theory/LinearArithmatic/Constraint.hs
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -78,29 +91,68 @@ 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)

let rel' = if coeff_of_x < 0 then flipRelation rel else rel

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
$ 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
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
204 changes: 119 additions & 85 deletions src/Theory/LinearArithmatic/FourierMotzkin.hs
Original file line number Diff line number Diff line change
@@ -1,102 +1,136 @@
{-|
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
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) ]
Loading

0 comments on commit 823bf57

Please sign in to comment.