{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Arith.Simplex.Textbook.LPSolver
-- Copyright   :  (c) Masahiro Sakai 2011
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- Naïve implementation of Simplex method
--
-- Reference:
--
-- * <http://www.math.cuhk.edu.hk/~wei/lpch3.pdf>
--
-----------------------------------------------------------------------------

module ToySolver.Arith.Simplex.Textbook.LPSolver
  (
  -- * Solver type
    Solver
  , emptySolver

  -- * LP monad
  , LP
  , getTableau
  , putTableau

  -- * Problem specification
  , newVar
  , addConstraint
  , addConstraintWithArtificialVariable
  , tableau
  , define

  -- * Solving
  , phaseI
  , simplex
  , dualSimplex
  , OptResult (..)
  , twoPhaseSimplex
  , primalDualSimplex

  -- * Extract results
  , getModel

  -- * Utilities
  , collectNonnegVars
  ) where

import Control.Exception (assert)
import Control.Monad
import Control.Monad.State
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import Data.OptDir
import Data.VectorSpace

import Data.Interval ((<=!), (>=!), (==!), (<!), (>!), (/=!))
import qualified Data.Interval as Interval

import ToySolver.Data.OrdRel
import qualified ToySolver.Data.LA as LA
import ToySolver.Data.IntVar
import qualified ToySolver.Arith.Simplex.Textbook as Simplex
import qualified ToySolver.Arith.BoundsInference as BI

-- ---------------------------------------------------------------------------
-- Solver type

type Solver r = (Var, Simplex.Tableau r, VarSet, VarMap (LA.Expr r))

emptySolver :: VarSet -> Solver r
emptySolver :: forall r. VarSet -> Solver r
emptySolver VarSet
vs = (Var
1 Var -> Var -> Var
forall a. Num a => a -> a -> a
+ [Var] -> Var
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ((-Var
1) Var -> [Var] -> [Var]
forall a. a -> [a] -> [a]
: VarSet -> [Var]
IS.toList VarSet
vs), Tableau r
forall r. Tableau r
Simplex.emptyTableau, VarSet
IS.empty, IntMap (Expr r)
forall a. IntMap a
IM.empty)

-- ---------------------------------------------------------------------------
-- LP Monad

type LP r = State (Solver r)

-- | Allocate a new /non-negative/ variable.
newVar :: LP r Var
newVar :: forall r. LP r Var
newVar = do
  (Var
x,Tableau r
tbl,VarSet
avs,VarMap (Expr r)
defs) <- StateT (Solver r) Identity (Solver r)
forall s (m :: * -> *). MonadState s m => m s
get
  Solver r -> StateT (Solver r) Identity ()
forall s (m :: * -> *). MonadState s m => s -> m ()
put (Var
xVar -> Var -> Var
forall a. Num a => a -> a -> a
+Var
1,Tableau r
tbl,VarSet
avs,VarMap (Expr r)
defs)
  Var -> LP r Var
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return Var
x

getTableau :: LP r (Simplex.Tableau r)
getTableau :: forall r. LP r (Tableau r)
getTableau = do
  (Var
_,Tableau r
tbl,VarSet
_,VarMap (Expr r)
_) <- StateT (Solver r) Identity (Solver r)
forall s (m :: * -> *). MonadState s m => m s
get
  Tableau r -> LP r (Tableau r)
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return Tableau r
tbl

putTableau :: Simplex.Tableau r -> LP r ()
putTableau :: forall r. Tableau r -> LP r ()
putTableau Tableau r
tbl = do
  (Var
x,Tableau r
_,VarSet
avs,VarMap (Expr r)
defs) <- StateT (Solver r) Identity (Solver r)
forall s (m :: * -> *). MonadState s m => m s
get
  Solver r -> LP r ()
forall s (m :: * -> *). MonadState s m => s -> m ()
put (Var
x,Tableau r
tbl,VarSet
avs,VarMap (Expr r)
defs)

addArtificialVariable :: Var -> LP r ()
addArtificialVariable :: forall r. Var -> LP r ()
addArtificialVariable Var
v = do
  (Var
x,Tableau r
tbl,VarSet
avs,VarMap (Expr r)
defs) <- StateT (Solver r) Identity (Solver r)
forall s (m :: * -> *). MonadState s m => m s
get
  Solver r -> LP r ()
forall s (m :: * -> *). MonadState s m => s -> m ()
put (Var
x, Tableau r
tbl, Var -> VarSet -> VarSet
IS.insert Var
v VarSet
avs, VarMap (Expr r)
defs)

getArtificialVariables :: LP r VarSet
getArtificialVariables :: forall r. LP r VarSet
getArtificialVariables = do
  (Var
_,Tableau r
_,VarSet
avs,VarMap (Expr r)
_) <- StateT (Solver r) Identity (Solver r)
forall s (m :: * -> *). MonadState s m => m s
get
  VarSet -> LP r VarSet
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return VarSet
avs

clearArtificialVariables :: LP r ()
clearArtificialVariables :: forall r. LP r ()
clearArtificialVariables = do
  (Var
x,Tableau r
tbl,VarSet
_,VarMap (Expr r)
defs) <- StateT (Solver r) Identity (Solver r)
forall s (m :: * -> *). MonadState s m => m s
get
  Solver r -> LP r ()
forall s (m :: * -> *). MonadState s m => s -> m ()
put (Var
x, Tableau r
tbl, VarSet
IS.empty, VarMap (Expr r)
defs)

define :: Var -> LA.Expr r -> LP r ()
define :: forall r. Var -> Expr r -> LP r ()
define Var
v Expr r
e = do
  (Var
x,Tableau r
tbl,VarSet
avs,IntMap (Expr r)
defs) <- StateT (Solver r) Identity (Solver r)
forall s (m :: * -> *). MonadState s m => m s
get
  Solver r -> LP r ()
forall s (m :: * -> *). MonadState s m => s -> m ()
put (Var
x,Tableau r
tbl,VarSet
avs, Var -> Expr r -> IntMap (Expr r) -> IntMap (Expr r)
forall a. Var -> a -> IntMap a -> IntMap a
IM.insert Var
v Expr r
e IntMap (Expr r)
defs)

getDefs :: LP r (VarMap (LA.Expr r))
getDefs :: forall r. LP r (VarMap (Expr r))
getDefs = do
  (Var
_,Tableau r
_,VarSet
_,VarMap (Expr r)
defs) <- StateT (Solver r) Identity (Solver r)
forall s (m :: * -> *). MonadState s m => m s
get
  VarMap (Expr r) -> LP r (VarMap (Expr r))
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return VarMap (Expr r)
defs

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

-- | Add a contraint and maintain feasibility condition by introducing artificial variable (if necessary).
--
-- * Disequality is not supported.
--
-- * Unlike 'addConstraint', an equality contstraint becomes one row with an artificial variable.
--
addConstraintWithArtificialVariable :: Real r => LA.Atom r -> LP r ()
addConstraintWithArtificialVariable :: forall r. Real r => Atom r -> LP r ()
addConstraintWithArtificialVariable Atom r
c = do
  Atom r
c2 <- Atom r -> LP r (Atom r)
forall r. (Num r, Eq r) => Atom r -> LP r (Atom r)
expandDefs' Atom r
c
  let (Expr r
e, RelOp
rop, r
b) = Atom r -> (Expr r, RelOp, r)
forall r. Real r => Atom r -> (Expr r, RelOp, r)
normalizeConstraint Atom r
c2
  Bool -> LP r () -> LP r ()
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (r
b r -> r -> Bool
forall a. Ord a => a -> a -> Bool
>= r
0) (LP r () -> LP r ()) -> LP r () -> LP r ()
forall a b. (a -> b) -> a -> b
$ () -> LP r ()
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return ()
  Tableau r
tbl <- LP r (Tableau r)
forall r. LP r (Tableau r)
getTableau
  case RelOp
rop of
    -- x≥0 is trivially true, since x is non-negative.
    RelOp
Ge | r
br -> r -> Bool
forall a. Eq a => a -> a -> Bool
==r
0 Bool -> Bool -> Bool
&& Expr r -> Bool
forall r. Real r => Expr r -> Bool
isSingleVar Expr r
e -> () -> LP r ()
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return ()
    -- -x≤0 is trivially true, since x is non-negative.
    RelOp
Le | r
br -> r -> Bool
forall a. Eq a => a -> a -> Bool
==r
0 Bool -> Bool -> Bool
&& Expr r -> Bool
forall r. Real r => Expr r -> Bool
isSingleNegatedVar Expr r
e -> () -> LP r ()
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return ()

    RelOp
Le -> do
      Var
v <- LP r Var
forall r. LP r Var
newVar -- slack variable
      Tableau r -> LP r ()
forall r. Tableau r -> LP r ()
putTableau (Tableau r -> LP r ()) -> Tableau r -> LP r ()
forall a b. (a -> b) -> a -> b
$ Tableau r -> Var -> Row r -> Tableau r
forall r. (Num r, Eq r) => Tableau r -> Var -> Row r -> Tableau r
Simplex.addRow Tableau r
tbl Var
v (Expr r -> IntMap r
forall r. Expr r -> IntMap r
LA.coeffMap Expr r
e, r
b)
    RelOp
Ge -> do
      Var
v1 <- LP r Var
forall r. LP r Var
newVar -- surplus variable
      Var
v2 <- LP r Var
forall r. LP r Var
newVar -- artificial variable
      Tableau r -> LP r ()
forall r. Tableau r -> LP r ()
putTableau (Tableau r -> LP r ()) -> Tableau r -> LP r ()
forall a b. (a -> b) -> a -> b
$ Tableau r -> Var -> Row r -> Tableau r
forall r. (Num r, Eq r) => Tableau r -> Var -> Row r -> Tableau r
Simplex.addRow Tableau r
tbl Var
v2 (Expr r -> IntMap r
forall r. Expr r -> IntMap r
LA.coeffMap (Expr r
e Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^-^ Var -> Expr r
forall r. Num r => Var -> Expr r
LA.var Var
v1), r
b)
      Var -> LP r ()
forall r. Var -> LP r ()
addArtificialVariable Var
v2
    RelOp
Eql -> do
      Var
v <- LP r Var
forall r. LP r Var
newVar -- artificial variable
      Tableau r -> LP r ()
forall r. Tableau r -> LP r ()
putTableau (Tableau r -> LP r ()) -> Tableau r -> LP r ()
forall a b. (a -> b) -> a -> b
$ Tableau r -> Var -> Row r -> Tableau r
forall r. (Num r, Eq r) => Tableau r -> Var -> Row r -> Tableau r
Simplex.addRow Tableau r
tbl Var
v (Expr r -> IntMap r
forall r. Expr r -> IntMap r
LA.coeffMap Expr r
e, r
b)
      Var -> LP r ()
forall r. Var -> LP r ()
addArtificialVariable Var
v
    RelOp
_ -> [Char] -> LP r ()
forall a. (?callStack::CallStack) => [Char] -> a
error ([Char] -> LP r ()) -> [Char] -> LP r ()
forall a b. (a -> b) -> a -> b
$ [Char]
"ToySolver.LPSolver.addConstraintWithArtificialVariable does not support " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ RelOp -> [Char]
forall a. Show a => a -> [Char]
show RelOp
rop

-- | Add a contraint, without maintaining feasibilty condition of tableaus.
--
-- * Disequality is not supported.
--
-- * Unlike 'addConstraintWithArtificialVariable', an equality constraint becomes two rows.
--
addConstraint :: Real r => LA.Atom r -> LP r ()
addConstraint :: forall r. Real r => Atom r -> LP r ()
addConstraint Atom r
c = do
  OrdRel Expr r
lhs RelOp
rop Expr r
rhs <- Atom r -> LP r (Atom r)
forall r. (Num r, Eq r) => Atom r -> LP r (Atom r)
expandDefs' Atom r
c
  let
    (r
b', Expr r
e) = Var -> Expr r -> (r, Expr r)
forall r. Num r => Var -> Expr r -> (r, Expr r)
LA.extract Var
LA.unitVar (Expr r
lhs Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^-^ Expr r
rhs)
    b :: r
b = - r
b'
  case RelOp
rop of
    RelOp
Le -> Expr r -> r -> LP r ()
forall {r}. Real r => Expr r -> r -> StateT (Solver r) Identity ()
f Expr r
e r
b
    RelOp
Ge -> Expr r -> r -> LP r ()
forall {r}. Real r => Expr r -> r -> StateT (Solver r) Identity ()
f (Expr r -> Expr r
forall v. AdditiveGroup v => v -> v
negateV Expr r
e) (r -> r
forall a. Num a => a -> a
negate r
b)
    RelOp
Eql -> do
      -- Unlike addConstraintWithArtificialVariable, an equality constraint becomes two rows.
      Expr r -> r -> LP r ()
forall {r}. Real r => Expr r -> r -> StateT (Solver r) Identity ()
f Expr r
e r
b
      Expr r -> r -> LP r ()
forall {r}. Real r => Expr r -> r -> StateT (Solver r) Identity ()
f (Expr r -> Expr r
forall v. AdditiveGroup v => v -> v
negateV Expr r
e) (r -> r
forall a. Num a => a -> a
negate r
b)
    RelOp
_ -> [Char] -> LP r ()
forall a. (?callStack::CallStack) => [Char] -> a
error ([Char] -> LP r ()) -> [Char] -> LP r ()
forall a b. (a -> b) -> a -> b
$ [Char]
"ToySolver.LPSolver.addConstraint does not support " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ RelOp -> [Char]
forall a. Show a => a -> [Char]
show RelOp
rop
  where
    -- -x≤b with b≥0 is trivially true.
    f :: Expr r -> r -> StateT (Solver r) Identity ()
f Expr r
e r
b | Expr r -> Bool
forall r. Real r => Expr r -> Bool
isSingleNegatedVar Expr r
e Bool -> Bool -> Bool
&& r
0 r -> r -> Bool
forall a. Ord a => a -> a -> Bool
<= r
b = () -> StateT (Solver r) Identity ()
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return ()
    f Expr r
e r
b = do
      Tableau r
tbl <- LP r (Tableau r)
forall r. LP r (Tableau r)
getTableau
      Var
v <- LP r Var
forall r. LP r Var
newVar -- slack variable
      Tableau r -> StateT (Solver r) Identity ()
forall r. Tableau r -> LP r ()
putTableau (Tableau r -> StateT (Solver r) Identity ())
-> Tableau r -> StateT (Solver r) Identity ()
forall a b. (a -> b) -> a -> b
$ Tableau r -> Var -> Row r -> Tableau r
forall r. (Num r, Eq r) => Tableau r -> Var -> Row r -> Tableau r
Simplex.addRow Tableau r
tbl Var
v (Expr r -> IntMap r
forall r. Expr r -> IntMap r
LA.coeffMap Expr r
e, r
b)

isSingleVar :: Real r => LA.Expr r -> Bool
isSingleVar :: forall r. Real r => Expr r -> Bool
isSingleVar Expr r
e =
  case Expr r -> [(r, Var)]
forall r. Expr r -> [(r, Var)]
LA.terms Expr r
e of
    [(r
1,Var
_)] -> Bool
True
    [(r, Var)]
_ -> Bool
False

isSingleNegatedVar :: Real r => LA.Expr r -> Bool
isSingleNegatedVar :: forall r. Real r => Expr r -> Bool
isSingleNegatedVar Expr r
e =
  case Expr r -> [(r, Var)]
forall r. Expr r -> [(r, Var)]
LA.terms Expr r
e of
    [(-1,Var
_)] -> Bool
True
    [(r, Var)]
_ -> Bool
False

expandDefs :: (Num r, Eq r) => LA.Expr r -> LP r (LA.Expr r)
expandDefs :: forall r. (Num r, Eq r) => Expr r -> LP r (Expr r)
expandDefs Expr r
e = do
  VarMap (Expr r)
defs <- LP r (VarMap (Expr r))
forall r. LP r (VarMap (Expr r))
getDefs
  Expr r -> LP r (Expr r)
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (Expr r -> LP r (Expr r)) -> Expr r -> LP r (Expr r)
forall a b. (a -> b) -> a -> b
$ VarMap (Expr r) -> Expr r -> Expr r
forall r. (Num r, Eq r) => VarMap (Expr r) -> Expr r -> Expr r
LA.applySubst VarMap (Expr r)
defs Expr r
e

expandDefs' :: (Num r, Eq r) => LA.Atom r -> LP r (LA.Atom r)
expandDefs' :: forall r. (Num r, Eq r) => Atom r -> LP r (Atom r)
expandDefs' (OrdRel Expr r
lhs RelOp
op Expr r
rhs) = do
  Expr r
lhs' <- Expr r -> LP r (Expr r)
forall r. (Num r, Eq r) => Expr r -> LP r (Expr r)
expandDefs Expr r
lhs
  Expr r
rhs' <- Expr r -> LP r (Expr r)
forall r. (Num r, Eq r) => Expr r -> LP r (Expr r)
expandDefs Expr r
rhs
  OrdRel (Expr r) -> LP r (OrdRel (Expr r))
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (OrdRel (Expr r) -> LP r (OrdRel (Expr r)))
-> OrdRel (Expr r) -> LP r (OrdRel (Expr r))
forall a b. (a -> b) -> a -> b
$ Expr r -> RelOp -> Expr r -> OrdRel (Expr r)
forall e. e -> RelOp -> e -> OrdRel e
OrdRel Expr r
lhs' RelOp
op Expr r
rhs'

tableau :: (RealFrac r) => [LA.Atom r] -> LP r ()
tableau :: forall r. RealFrac r => [Atom r] -> LP r ()
tableau [Atom r]
cs = do
  let (VarSet
nonnegVars, [Atom r]
cs') = [Atom r] -> VarSet -> (VarSet, [Atom r])
forall r. RealFrac r => [Atom r] -> VarSet -> (VarSet, [Atom r])
collectNonnegVars [Atom r]
cs VarSet
IS.empty
      fvs :: VarSet
fvs = [Atom r] -> VarSet
forall a. Variables a => a -> VarSet
vars [Atom r]
cs VarSet -> VarSet -> VarSet
`IS.difference` VarSet
nonnegVars
  [Var] -> (Var -> LP r ()) -> LP r ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ (VarSet -> [Var]
IS.toList VarSet
fvs) ((Var -> LP r ()) -> LP r ()) -> (Var -> LP r ()) -> LP r ()
forall a b. (a -> b) -> a -> b
$ \Var
v -> do
    Var
v1 <- LP r Var
forall r. LP r Var
newVar
    Var
v2 <- LP r Var
forall r. LP r Var
newVar
    Var -> Expr r -> LP r ()
forall r. Var -> Expr r -> LP r ()
define Var
v (Var -> Expr r
forall r. Num r => Var -> Expr r
LA.var Var
v1 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^-^ Var -> Expr r
forall r. Num r => Var -> Expr r
LA.var Var
v2)
  (Atom r -> LP r ()) -> [Atom r] -> LP r ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Atom r -> LP r ()
forall r. Real r => Atom r -> LP r ()
addConstraint [Atom r]
cs'

getModel :: Fractional r => VarSet -> LP r (Model r)
getModel :: forall r. Fractional r => VarSet -> LP r (Model r)
getModel VarSet
vs = do
  Tableau r
tbl <- LP r (Tableau r)
forall r. LP r (Tableau r)
getTableau
  VarMap (Expr r)
defs <- LP r (VarMap (Expr r))
forall r. LP r (VarMap (Expr r))
getDefs
  let vs' :: VarSet
vs' = (VarSet
vs VarSet -> VarSet -> VarSet
`IS.difference` VarMap (Expr r) -> VarSet
forall a. IntMap a -> VarSet
IM.keysSet VarMap (Expr r)
defs) VarSet -> VarSet -> VarSet
`IS.union` [VarSet] -> VarSet
forall (f :: * -> *). Foldable f => f VarSet -> VarSet
IS.unions [Expr r -> VarSet
forall a. Variables a => a -> VarSet
vars Expr r
e | Expr r
e <- VarMap (Expr r) -> [Expr r]
forall a. IntMap a -> [a]
IM.elems VarMap (Expr r)
defs]
      m0 :: IntMap r
m0 = [(Var, r)] -> IntMap r
forall a. [(Var, a)] -> IntMap a
IM.fromAscList [(Var
v, Tableau r -> Var -> r
forall r. Num r => Tableau r -> Var -> r
Simplex.currentValue Tableau r
tbl Var
v) | Var
v <- VarSet -> [Var]
IS.toAscList VarSet
vs']
  IntMap r -> LP r (IntMap r)
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (IntMap r -> LP r (IntMap r)) -> IntMap r -> LP r (IntMap r)
forall a b. (a -> b) -> a -> b
$ (Var -> r -> Bool) -> IntMap r -> IntMap r
forall a. (Var -> a -> Bool) -> IntMap a -> IntMap a
IM.filterWithKey (\Var
k r
_ -> Var
k Var -> VarSet -> Bool
`IS.member` VarSet
vs) (IntMap r -> IntMap r) -> IntMap r -> IntMap r
forall a b. (a -> b) -> a -> b
$ (Expr r -> r) -> VarMap (Expr r) -> IntMap r
forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map (IntMap r -> Expr r -> r
forall m e v. Eval m e v => m -> e -> v
LA.eval IntMap r
m0) VarMap (Expr r)
defs IntMap r -> IntMap r -> IntMap r
forall a. IntMap a -> IntMap a -> IntMap a
`IM.union` IntMap r
m0

phaseI :: (Fractional r, Real r) => LP r Bool
phaseI :: forall r. (Fractional r, Real r) => LP r Bool
phaseI = do
  LP r ()
forall r. Real r => LP r ()
introduceArtificialVariables
  Tableau r
tbl <- LP r (Tableau r)
forall r. LP r (Tableau r)
getTableau
  VarSet
avs <- LP r VarSet
forall r. LP r VarSet
getArtificialVariables
  let (Bool
ret, Tableau r
tbl') = Tableau r -> VarSet -> (Bool, Tableau r)
forall r.
(Real r, Fractional r) =>
Tableau r -> VarSet -> (Bool, Tableau r)
Simplex.phaseI Tableau r
tbl VarSet
avs
  Tableau r -> LP r ()
forall r. Tableau r -> LP r ()
putTableau Tableau r
tbl'
  Bool -> LP r () -> LP r ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
ret LP r ()
forall r. LP r ()
clearArtificialVariables
  Bool -> LP r Bool
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return Bool
ret

introduceArtificialVariables :: (Real r) => LP r ()
introduceArtificialVariables :: forall r. Real r => LP r ()
introduceArtificialVariables = do
  Tableau r
tbl <- LP r (Tableau r)
forall r. LP r (Tableau r)
getTableau
  Tableau r
tbl' <- ([(Var, (VarMap r, r))] -> Tableau r)
-> StateT (Solver r) Identity [(Var, (VarMap r, r))]
-> LP r (Tableau r)
forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
liftM [(Var, (VarMap r, r))] -> Tableau r
forall a. [(Var, a)] -> IntMap a
IM.fromList (StateT (Solver r) Identity [(Var, (VarMap r, r))]
 -> LP r (Tableau r))
-> StateT (Solver r) Identity [(Var, (VarMap r, r))]
-> LP r (Tableau r)
forall a b. (a -> b) -> a -> b
$ [(Var, (VarMap r, r))]
-> ((Var, (VarMap r, r))
    -> StateT (Solver r) Identity (Var, (VarMap r, r)))
-> StateT (Solver r) Identity [(Var, (VarMap r, r))]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM (Tableau r -> [(Var, (VarMap r, r))]
forall a. IntMap a -> [(Var, a)]
IM.toList Tableau r
tbl) (((Var, (VarMap r, r))
  -> StateT (Solver r) Identity (Var, (VarMap r, r)))
 -> StateT (Solver r) Identity [(Var, (VarMap r, r))])
-> ((Var, (VarMap r, r))
    -> StateT (Solver r) Identity (Var, (VarMap r, r)))
-> StateT (Solver r) Identity [(Var, (VarMap r, r))]
forall a b. (a -> b) -> a -> b
$ \(Var
v,(VarMap r
e,r
rhs)) -> do
    if r
rhs r -> r -> Bool
forall a. Ord a => a -> a -> Bool
>= r
0 then do
      (Var, (VarMap r, r))
-> StateT (Solver r) Identity (Var, (VarMap r, r))
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (Var
v,(VarMap r
e,r
rhs)) -- v + e == rhs
    else do
      Var
a <- LP r Var
forall r. LP r Var
newVar
      Var -> LP r ()
forall r. Var -> LP r ()
addArtificialVariable Var
a
      (Var, (VarMap r, r))
-> StateT (Solver r) Identity (Var, (VarMap r, r))
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (Var
a, (Var -> r -> VarMap r -> VarMap r
forall a. Var -> a -> IntMap a -> IntMap a
IM.insert Var
v (-r
1) ((r -> r) -> VarMap r -> VarMap r
forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map r -> r
forall a. Num a => a -> a
negate VarMap r
e), -r
rhs)) -- a - (v + e) == -rhs
  Tableau r -> LP r ()
forall r. Tableau r -> LP r ()
putTableau Tableau r
tbl'

simplex :: (Fractional r, Real r) => OptDir -> LA.Expr r -> LP r Bool
simplex :: forall r. (Fractional r, Real r) => OptDir -> Expr r -> LP r Bool
simplex OptDir
optdir Expr r
obj = do
  Tableau r
tbl <- LP r (Tableau r)
forall r. LP r (Tableau r)
getTableau
  VarMap (Expr r)
defs <- LP r (VarMap (Expr r))
forall r. LP r (VarMap (Expr r))
getDefs
  let (Bool
ret, Tableau r
tbl') = OptDir -> Tableau r -> (Bool, Tableau r)
forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> (Bool, Tableau r)
Simplex.simplex OptDir
optdir (Tableau r -> Expr r -> Tableau r
forall r. (Num r, Eq r) => Tableau r -> Expr r -> Tableau r
Simplex.setObjFun Tableau r
tbl (VarMap (Expr r) -> Expr r -> Expr r
forall r. (Num r, Eq r) => VarMap (Expr r) -> Expr r -> Expr r
LA.applySubst VarMap (Expr r)
defs Expr r
obj))
  Tableau r -> LP r ()
forall r. Tableau r -> LP r ()
putTableau Tableau r
tbl'
  Bool -> LP r Bool
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return Bool
ret

dualSimplex :: (Fractional r, Real r) => OptDir -> LA.Expr r -> LP r Bool
dualSimplex :: forall r. (Fractional r, Real r) => OptDir -> Expr r -> LP r Bool
dualSimplex OptDir
optdir Expr r
obj = do
  Tableau r
tbl <- LP r (Tableau r)
forall r. LP r (Tableau r)
getTableau
  VarMap (Expr r)
defs <- LP r (VarMap (Expr r))
forall r. LP r (VarMap (Expr r))
getDefs
  let (Bool
ret, Tableau r
tbl') = OptDir -> Tableau r -> (Bool, Tableau r)
forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> (Bool, Tableau r)
Simplex.dualSimplex OptDir
optdir (Tableau r -> Expr r -> Tableau r
forall r. (Num r, Eq r) => Tableau r -> Expr r -> Tableau r
Simplex.setObjFun Tableau r
tbl (VarMap (Expr r) -> Expr r -> Expr r
forall r. (Num r, Eq r) => VarMap (Expr r) -> Expr r -> Expr r
LA.applySubst VarMap (Expr r)
defs Expr r
obj))
  Tableau r -> LP r ()
forall r. Tableau r -> LP r ()
putTableau Tableau r
tbl'
  Bool -> LP r Bool
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return Bool
ret

-- | results of optimization
data OptResult = Optimum | Unsat | Unbounded
  deriving (Var -> OptResult -> [Char] -> [Char]
[OptResult] -> [Char] -> [Char]
OptResult -> [Char]
(Var -> OptResult -> [Char] -> [Char])
-> (OptResult -> [Char])
-> ([OptResult] -> [Char] -> [Char])
-> Show OptResult
forall a.
(Var -> a -> [Char] -> [Char])
-> (a -> [Char]) -> ([a] -> [Char] -> [Char]) -> Show a
$cshowsPrec :: Var -> OptResult -> [Char] -> [Char]
showsPrec :: Var -> OptResult -> [Char] -> [Char]
$cshow :: OptResult -> [Char]
show :: OptResult -> [Char]
$cshowList :: [OptResult] -> [Char] -> [Char]
showList :: [OptResult] -> [Char] -> [Char]
Show, OptResult -> OptResult -> Bool
(OptResult -> OptResult -> Bool)
-> (OptResult -> OptResult -> Bool) -> Eq OptResult
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: OptResult -> OptResult -> Bool
== :: OptResult -> OptResult -> Bool
$c/= :: OptResult -> OptResult -> Bool
/= :: OptResult -> OptResult -> Bool
Eq, Eq OptResult
Eq OptResult =>
(OptResult -> OptResult -> Ordering)
-> (OptResult -> OptResult -> Bool)
-> (OptResult -> OptResult -> Bool)
-> (OptResult -> OptResult -> Bool)
-> (OptResult -> OptResult -> Bool)
-> (OptResult -> OptResult -> OptResult)
-> (OptResult -> OptResult -> OptResult)
-> Ord OptResult
OptResult -> OptResult -> Bool
OptResult -> OptResult -> Ordering
OptResult -> OptResult -> OptResult
forall a.
Eq a =>
(a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: OptResult -> OptResult -> Ordering
compare :: OptResult -> OptResult -> Ordering
$c< :: OptResult -> OptResult -> Bool
< :: OptResult -> OptResult -> Bool
$c<= :: OptResult -> OptResult -> Bool
<= :: OptResult -> OptResult -> Bool
$c> :: OptResult -> OptResult -> Bool
> :: OptResult -> OptResult -> Bool
$c>= :: OptResult -> OptResult -> Bool
>= :: OptResult -> OptResult -> Bool
$cmax :: OptResult -> OptResult -> OptResult
max :: OptResult -> OptResult -> OptResult
$cmin :: OptResult -> OptResult -> OptResult
min :: OptResult -> OptResult -> OptResult
Ord)

twoPhaseSimplex :: (Fractional r, Real r) => OptDir -> LA.Expr r -> LP r OptResult
twoPhaseSimplex :: forall r.
(Fractional r, Real r) =>
OptDir -> Expr r -> LP r OptResult
twoPhaseSimplex OptDir
optdir Expr r
obj = do
  Bool
ret <- LP r Bool
forall r. (Fractional r, Real r) => LP r Bool
phaseI
  if Bool -> Bool
not Bool
ret then
    OptResult -> LP r OptResult
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Unsat
  else do
    Bool
ret <- OptDir -> Expr r -> LP r Bool
forall r. (Fractional r, Real r) => OptDir -> Expr r -> LP r Bool
simplex OptDir
optdir Expr r
obj
    if Bool
ret then
      OptResult -> LP r OptResult
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Optimum
    else
      OptResult -> LP r OptResult
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Unbounded

primalDualSimplex :: (Fractional r, Real r) => OptDir -> LA.Expr r -> LP r OptResult
primalDualSimplex :: forall r.
(Fractional r, Real r) =>
OptDir -> Expr r -> LP r OptResult
primalDualSimplex OptDir
optdir Expr r
obj = do
  Tableau r
tbl <- LP r (Tableau r)
forall r. LP r (Tableau r)
getTableau
  VarMap (Expr r)
defs <- LP r (VarMap (Expr r))
forall r. LP r (VarMap (Expr r))
getDefs
  let (Bool
ret, Tableau r
tbl') = OptDir -> Tableau r -> (Bool, Tableau r)
forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> (Bool, Tableau r)
Simplex.primalDualSimplex OptDir
optdir (Tableau r -> Expr r -> Tableau r
forall r. (Num r, Eq r) => Tableau r -> Expr r -> Tableau r
Simplex.setObjFun Tableau r
tbl (VarMap (Expr r) -> Expr r -> Expr r
forall r. (Num r, Eq r) => VarMap (Expr r) -> Expr r -> Expr r
LA.applySubst VarMap (Expr r)
defs Expr r
obj))
  Tableau r -> LP r ()
forall r. Tableau r -> LP r ()
putTableau Tableau r
tbl'
  if Bool
ret then
    OptResult -> LP r OptResult
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Optimum
  else if Bool -> Bool
not (Tableau r -> Bool
forall r. Real r => Tableau r -> Bool
Simplex.isFeasible Tableau r
tbl') then
    OptResult -> LP r OptResult
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Unsat
  else
    OptResult -> LP r OptResult
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Unbounded

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

-- convert right hand side to be non-negative
normalizeConstraint :: forall r. Real r => LA.Atom r -> (LA.Expr r, RelOp, r)
normalizeConstraint :: forall r. Real r => Atom r -> (Expr r, RelOp, r)
normalizeConstraint (OrdRel Expr r
a RelOp
op Expr r
b)
  | r
rhs r -> r -> Bool
forall a. Ord a => a -> a -> Bool
< r
0   = (Expr r -> Expr r
forall v. AdditiveGroup v => v -> v
negateV Expr r
lhs, RelOp -> RelOp
flipOp RelOp
op, r -> r
forall a. Num a => a -> a
negate r
rhs)
  | Bool
otherwise = (Expr r
lhs, RelOp
op, r
rhs)
  where
    (r
c, Expr r
lhs) = Var -> Expr r -> (r, Expr r)
forall r. Num r => Var -> Expr r -> (r, Expr r)
LA.extract Var
LA.unitVar (Expr r
a Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^-^ Expr r
b)
    rhs :: r
rhs = - r
c

collectNonnegVars :: forall r. (RealFrac r) => [LA.Atom r] -> VarSet -> (VarSet, [LA.Atom r])
collectNonnegVars :: forall r. RealFrac r => [Atom r] -> VarSet -> (VarSet, [Atom r])
collectNonnegVars [Atom r]
cs VarSet
ivs = (VarSet
nonnegVars, [Atom r]
cs)
  where
    vs :: VarSet
vs = [Atom r] -> VarSet
forall a. Variables a => a -> VarSet
vars [Atom r]
cs
    bounds :: BoundsEnv r
bounds = BoundsEnv r -> [Atom r] -> VarSet -> Var -> BoundsEnv r
forall r.
RealFrac r =>
BoundsEnv r -> [Atom r] -> VarSet -> Var -> BoundsEnv r
BI.inferBounds BoundsEnv r
initialBounds [Atom r]
cs VarSet
ivs Var
1000
      where
        initialBounds :: BoundsEnv r
initialBounds = [(Var, Interval r)] -> BoundsEnv r
forall a. [(Var, a)] -> IntMap a
IM.fromAscList [(Var
v, Interval r
forall r. Ord r => Interval r
Interval.whole) | Var
v <- VarSet -> [Var]
IS.toAscList VarSet
vs]
    nonnegVars :: VarSet
nonnegVars = (Var -> Bool) -> VarSet -> VarSet
IS.filter (\Var
v -> Interval r
0 Interval r -> Interval r -> Bool
forall r. Ord r => Interval r -> Interval r -> Bool
<=! (BoundsEnv r
bounds BoundsEnv r -> Var -> Interval r
forall a. IntMap a -> Var -> a
IM.! Var
v)) VarSet
vs

    isTriviallyTrue :: LA.Atom r -> Bool
    isTriviallyTrue :: Atom r -> Bool
isTriviallyTrue (OrdRel Expr r
a RelOp
op Expr r
b) =
      case RelOp
op of
        RelOp
Le -> Interval r
i Interval r -> Interval r -> Bool
forall r. Ord r => Interval r -> Interval r -> Bool
<=! Interval r
0
        RelOp
Ge -> Interval r
i Interval r -> Interval r -> Bool
forall r. Ord r => Interval r -> Interval r -> Bool
>=! Interval r
0
        RelOp
Lt -> Interval r
i Interval r -> Interval r -> Bool
forall r. Ord r => Interval r -> Interval r -> Bool
<! Interval r
0
        RelOp
Gt -> Interval r
i Interval r -> Interval r -> Bool
forall r. Ord r => Interval r -> Interval r -> Bool
>! Interval r
0
        RelOp
Eql -> Interval r
i Interval r -> Interval r -> Bool
forall r. Ord r => Interval r -> Interval r -> Bool
==! Interval r
0
        RelOp
NEq -> Interval r
i Interval r -> Interval r -> Bool
forall r. Ord r => Interval r -> Interval r -> Bool
/=! Interval r
0
      where
        i :: Interval r
i = BoundsEnv r -> Expr r -> Interval r
forall r.
(Real r, Fractional r) =>
BoundsEnv r -> Expr r -> Interval r
LA.computeInterval BoundsEnv r
bounds (Expr r
a Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^-^ Expr r
b)

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