Skip to content

Commit

Permalink
work on Simplex
Browse files Browse the repository at this point in the history
  • Loading branch information
@gruhn
gruhn committed Mar 15, 2023
1 parent b2fdddb commit 4cd7735
Show file tree
Hide file tree
Showing 5 changed files with 234 additions and 95 deletions.
2 changes: 2 additions & 0 deletions SMT.cabal
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,12 +43,14 @@ library
default-extensions: ScopedTypeVariables
, LambdaCase
, StrictData
, TupleSections

test-suite test
type: exitcode-stdio-1.0
main-is: Test.hs
other-modules: GameUnequal
, TestNonLinearRealArithmatic
, TestLinearArithmatic
build-depends: base ^>=4.14.3.0
, SMT
, tasty
Expand Down
275 changes: 183 additions & 92 deletions src/Theory/LinearArithmatic/Simplex.hs
View file
Original file line number Diff line number Diff line change
@@ -1,29 +1,146 @@
{-# LANGUAGE TupleSections #-}
module Theory.LinearArithmatic.Simplex () where
module Theory.LinearArithmatic.Simplex ( simplex, Assignment, Constraint, ConstraintRelation(..), eval ) where

import qualified Data.IntMap as M
import qualified Data.IntMap.Merge.Lazy as MM
import qualified Data.IntSet as S
import Control.Monad (unless, guard)
import Data.Bifunctor (bimap)
import Debug.Trace
import Data.Maybe (fromMaybe)

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 a = (LinearTerm a, ConstraintRelation, a)

-- | Map from variable IDs to assigned values
type Assignment a = M.IntMap a

-- | Map from variable IDs to coefficients
type LinearTerm a = M.IntMap a

data BoundType = UpperBound | LowerBound
deriving Show

data Tableau a = Tableau
{ getNonBasis :: M.IntMap (LinearTerm a)
{ getBasis :: M.IntMap (LinearTerm a)
, getBounds :: M.IntMap (BoundType, a)
, getAssignment :: Assignment a
}
} deriving Show

varsIn :: [Constraint a] -> S.IntSet
varsIn constraints = S.unions $ do
(linear_term, _, _) <- constraints
return $ M.keysSet linear_term

{-|
Transform constraints to "General From" and initialize Simplex tableau.
-}
initTableau :: forall a. Num a => [Constraint a] -> Tableau a
initTableau constraints =
let
original_vars = varsIn 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

let bound_type =
case relation of
LessEquals -> UpperBound
GreaterEquals -> LowerBound
other_rels -> error "TODO: extend Simplex to all constraint relations"

return
( (slack_var, linear_term)
, (slack_var, (bound_type, bound)))

slack_vars = M.keys bounds

assignment :: Assignment a
assignment = M.fromList
$ zip (S.toList original_vars ++ slack_vars)
$ repeat 0

in
Tableau basis bounds assignment

data BoundViolation = MustIncrease | MustDecrease

violatedBasicVars :: Ord a => Tableau a -> [ (Var, BoundViolation) ]
violatedBasicVars (Tableau basis bounds assignment) = do
-- Following "Blend's Rule" by enumerating variables in ascending order.
basic_var <- M.keys basis

let current_value = assignment M.! basic_var

case M.lookup basic_var bounds of
Just (UpperBound, bound) -> do
guard (current_value > bound)
return (basic_var, MustDecrease)
Just (LowerBound, bound) -> do
guard (current_value < bound)
return (basic_var, MustIncrease)
Nothing -> []

pivotCandidates :: (Ord a, Num a) => Tableau a -> [ (Var, Var) ]
pivotCandidates tableau@(Tableau basis bounds assignment) = do
(basic_var, violation) <- violatedBasicVars tableau

let term = 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
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
-- If the coefficient of the non basic variable is 0, then the value of the variable
-- is not affected by pivoting, so it can't resolve the bound violation.
(_, 0, _) -> False

-- If the non-basic variable is not bounded and it's coefficient is non-zero,
-- then pivot is always possible.
(Nothing, 1, _) -> True
(Nothing, -1, _) -> True

-- The value of the basic variable must be decreased. The coefficient of the non-basic
-- variable is positive, so we need to decrease it's value. The non-basic variable is
-- at it's upper bound so decreasing it is possible.
(Just UpperBound, 1, MustDecrease) -> True

-- The value of the basic variable must be increased. The coefficient of the non-basic
-- variable is negative, so we need to decrease it's value. The non-basic variable is
-- at it's upper bound so decreasing it is possible.
(Just UpperBound, -1, MustIncrease) -> True

-- The value of the basic variable must be increased. The coefficient of the non-basic
-- variable is positive, so we need to increase it's value. The non-basic variable is
-- at it's lower bound so increasing it is possible.
(Just LowerBound, 1, MustIncrease) -> True

-- The value of the basic variable must be decreased. The coefficient of the non-basic
-- variable is negative, so we need to increase it's value. The non-basic variable is
-- at it's lower bound so increasing it is possible.
(Just LowerBound, -1, MustDecrease) -> True

-- In all other cases the bound of the non-basic variable would be violated by pivoting.
all_other_cases -> False

guard can_pivot
return (basic_var, non_basic_var)

type Equation a = (Var, LinearTerm a)

{-|
Expand All @@ -36,7 +153,7 @@ type Equation a = (Var, LinearTerm a)
x = 1/3 y - 10/3
-}
solveFor :: Equation a -> Var -> Maybe (Equation a)
solveFor :: (Fractional a, Eq a) => Equation a -> Var -> Maybe (Equation a)
solveFor (y, right_hand_side) x = do
coeff_of_x <- M.lookup x right_hand_side
guard (coeff_of_x /= 0)
Expand All @@ -58,7 +175,7 @@ solveFor (y, right_hand_side) x = do
y = 3(w + 2z) + 4z
= 3w + 10z
-}
rewrite :: Equation a -> Equation a -> Equation a
rewrite :: Num a => Equation a -> Equation a -> Equation a
rewrite (x, term_x) (y, term_y) =
let
coeff_of_x = M.findWithDefault 0 x term_y
Expand All @@ -67,8 +184,14 @@ rewrite (x, term_x) (y, term_y) =
in
(y, term_y')

eval :: Assignment a -> LinearTerm a -> a
eval assignment term = sum . MM.zipWithMatched (const (*)) assignment term
eval :: Num a => Assignment a -> LinearTerm a -> a
eval assignment term = sum $
MM.merge
MM.dropMissing
MM.dropMissing
(MM.zipWithMatched (const (*)))
assignment
term

{-|
TODO: make sure:
Expand All @@ -77,104 +200,72 @@ eval assignment term = sum . MM.zipWithMatched (const (*)) assignment term
- non-basic variables must violate bound
- basic variable is suitable for pivoting
-}
pivot' :: Var -> Var -> Tableau a -> Tableau a
pivot' basic_var non_basic_var (Tableau non_basis bounds assignment) =
pivot :: (Fractional a, Eq a) => Var -> Var -> Tableau a -> Tableau a
pivot basic_var non_basic_var (Tableau basis bounds assignment) =
let
from_just msg (Just a) = a
from_just msg Nothing = error msg

term = from_just "Given variable is not actually in the non-basis ==> invalid pivot pair"
$ M.lookup non_basic_var non_basis
term = 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 basic variable ==> invalid pivot pair"
$ solveFor (non_basic_var, term) basic_var
equation = from_just "Can't solve for non-basic variable ==> invalid pivot pair"
$ solveFor (basic_var, term) non_basic_var

non_basis' = M.fromList $ rewrite equation <$> M.toList non_basis
basis' = M.fromList $ rewrite equation <$> M.toList 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

basic_var_coeff = snd equation M.! basic_var
non_basic_var_coeff = term M.! non_basic_var

old_value_non_basic_var = non_basis M.! non_basic_bar
new_value_non_basic_var = old_value_non_basic_var + (old_value_basic_bar - new_value_basic_var) / 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.union (eval <$> non_basis')
$ M.insert non_basic_var new_value_non_basic_var
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 non_basis' bounds assignment'
Tableau basis' bounds assignment''

pivot :: forall a. (Num a, Ord a) => Tableau a -> Tableau a
pivot (Tableau non_basis bounds assignment) =
{-|
TODO:
in case of UNSAT include explanation, i.e. minimal infeasible subset of constraints.
-}
simplex :: forall a. (Num a, Ord a, Fractional a, Show a) => [Constraint a] -> Maybe (Assignment a)
simplex constraints =
let
basis = M.keysSet assignment S.\\ M.keysSet non_basis

violated_basic_vars :: [ (Var, BoundViolation) ]
violated_basic_vars = do
-- Following "Blend's Rule" by enumerating variables in ascending order.
basic_var <- S.toAscList basis

let current_value = assignment M.! basic_var

case M.lookup basic_var bounds of
Just (UpperBound bound) -> do
guard (current_value > bound)
return (basic_var, MustDecrease)
Just (LowerBound bound) -> do
guard (current_value < bound)
return (basic_var, MustIncrease)
Nothing -> []

pivot_candidates :: [ (Var, Var) ]
pivot_candidates = do
(basic_var, violation) <- violated_basic_vars
(non_basic_var, term) <- M.toAscList non_basis

let basic_var_coeff = M.findWithDefault 0 basic_var term
bound_type = fst <$> M.lookup non_basic_var bounds
can_pivot =
case (bound_type, signum basic_var_coeff, violation) of
-- If the coefficient of the basic variable is 0, then the value of the variable
-- is not affected by pivoting, so it can't resolve the bound violation.
(_, 0, _) -> False

-- If the coefficient of the basic variable is non-zero and the non-basic
-- variable is not bounded, pivot is always possible.
(Nothing, 1, _) -> True
(Nothing, -1, _) -> True

-- The non-basic variable is at it's upper bound so we can only decrease it,
-- without violating it. The value of the basic variable must be decreased,
-- and since the the coefficient is postive, decreasing it makes the non-basic
-- variable value smaller.
(Just UpperBound, 1, MustDecrease) -> True

-- The non-basic variable is at it's upper bound so we can only decrease it,
-- without violating it. The value of the basic variable must be increased,
-- but the coefficient is negative, so increasing it actually makes the non-basic
-- variable value smaller.
(Just UpperBound, -1, MustIncrease) -> True

-- The non-basic variable is at it's lower bound so we can only increase it,
-- without violating it. The value of the basic variable must be increased,
-- and since the the coefficient is postive, decreasing it makes the non-basic
-- variable value larger.
(Just LowerBound, 1, MustIncrease) -> True

-- The non-basic variable is at it's lower bound so we can only increase it,
-- without violating it. The value of the basic variable must be decreased,
-- but the coefficient is negative, so decreasing it actually makes the non-basic
-- variable value larger.
(Just LowerBound, -1, MustDecrease) -> True

-- In all other cases the bound of the non-basic variable would be violated by
-- pivoting.
all_other_cases -> False

guard can_pivot
return (basic_var, non_basic_var)

in undefined

go :: Tableau a -> Tableau a
go tableau =
case pivotCandidates $ traceShowId $ tableau of
[] -> tableau
((basic_var, non_basic_var) : _) ->
go (pivot basic_var non_basic_var tableau)

tableau_0 = initTableau constraints
tableau_n = go tableau_0

original_vars = varsIn constraints
assignment = M.restrictKeys (getAssignment tableau_n) original_vars
in
if null (violatedBasicVars tableau_n) then
Just assignment
else
Nothing

-----------------------------------------------------------

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)
-- , (M.fromList [(0,1), (1,-1)], LessEquals, -1)
]

example2 = simplex
[ (M.fromList [(0,1), (1,1)], GreaterEquals, 2) ]

4 changes: 2 additions & 2 deletions test/Test.hs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,14 +4,14 @@ import Test.Tasty
import Test.Tasty.HUnit
import Test.Tasty.QuickCheck

-- import qualified TestPropositions as TestProp
import qualified TestNonLinearRealArithmatic as TestNRA
import qualified TestLinearArithmatic as TestLA

main :: IO ()
main = defaultMain tests

tests :: TestTree
tests = testGroup "All tests"
[ TestNRA.tests
-- , TestProp.tests
, TestLA.tests
]
Loading

0 comments on commit 4cd7735

Please sign in to comment.