{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Arith.Simplex.Textbook
-- 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
  (
  -- * The @Tableau@ type
    Tableau
  , RowIndex
  , ColIndex
  , Row
  , emptyTableau
  , objRowIndex
  , pivot
  , lookupRow
  , addRow
  , setObjFun

  -- * Optimization direction
  , module Data.OptDir

  -- * Reading status
  , currentValue
  , currentObjValue
  , isFeasible
  , isOptimal

  -- * High-level solving functions
  , simplex
  , dualSimplex
  , phaseI
  , primalDualSimplex

  -- * For debugging
  , isValidTableau
  , toCSV
  ) where

import Data.Ord
import Data.List
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import Data.OptDir
import Data.VectorSpace
import Control.Exception

import qualified ToySolver.Data.LA as LA
import ToySolver.Data.IntVar

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

type Tableau r = VarMap (Row r)
{-
tbl ! v == (m, val)
==>
var v .+. m .==. constant val
-}

-- | Basic variables
type RowIndex = Int

-- | Non-basic variables
type ColIndex = Int

type Row r = (VarMap r, r)

data PivotResult r = PivotUnbounded | PivotFinished | PivotSuccess (Tableau r)
  deriving (RowIndex -> PivotResult r -> ShowS
forall r. Show r => RowIndex -> PivotResult r -> ShowS
forall r. Show r => [PivotResult r] -> ShowS
forall r. Show r => PivotResult r -> String
forall a.
(RowIndex -> a -> ShowS)
-> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [PivotResult r] -> ShowS
$cshowList :: forall r. Show r => [PivotResult r] -> ShowS
show :: PivotResult r -> String
$cshow :: forall r. Show r => PivotResult r -> String
showsPrec :: RowIndex -> PivotResult r -> ShowS
$cshowsPrec :: forall r. Show r => RowIndex -> PivotResult r -> ShowS
Show, PivotResult r -> PivotResult r -> Bool
forall r. Eq r => PivotResult r -> PivotResult r -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: PivotResult r -> PivotResult r -> Bool
$c/= :: forall r. Eq r => PivotResult r -> PivotResult r -> Bool
== :: PivotResult r -> PivotResult r -> Bool
$c== :: forall r. Eq r => PivotResult r -> PivotResult r -> Bool
Eq, PivotResult r -> PivotResult r -> Bool
PivotResult r -> PivotResult r -> Ordering
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
forall {r}. Ord r => Eq (PivotResult r)
forall r. Ord r => PivotResult r -> PivotResult r -> Bool
forall r. Ord r => PivotResult r -> PivotResult r -> Ordering
forall r. Ord r => PivotResult r -> PivotResult r -> PivotResult r
min :: PivotResult r -> PivotResult r -> PivotResult r
$cmin :: forall r. Ord r => PivotResult r -> PivotResult r -> PivotResult r
max :: PivotResult r -> PivotResult r -> PivotResult r
$cmax :: forall r. Ord r => PivotResult r -> PivotResult r -> PivotResult r
>= :: PivotResult r -> PivotResult r -> Bool
$c>= :: forall r. Ord r => PivotResult r -> PivotResult r -> Bool
> :: PivotResult r -> PivotResult r -> Bool
$c> :: forall r. Ord r => PivotResult r -> PivotResult r -> Bool
<= :: PivotResult r -> PivotResult r -> Bool
$c<= :: forall r. Ord r => PivotResult r -> PivotResult r -> Bool
< :: PivotResult r -> PivotResult r -> Bool
$c< :: forall r. Ord r => PivotResult r -> PivotResult r -> Bool
compare :: PivotResult r -> PivotResult r -> Ordering
$ccompare :: forall r. Ord r => PivotResult r -> PivotResult r -> Ordering
Ord)

emptyTableau :: Tableau r
emptyTableau :: forall r. Tableau r
emptyTableau = forall a. IntMap a
IM.empty

objRowIndex :: RowIndex
objRowIndex :: RowIndex
objRowIndex = -RowIndex
1

pivot :: (Fractional r, Eq r) => RowIndex -> ColIndex -> Tableau r -> Tableau r
{-# INLINE pivot #-}
{-# SPECIALIZE pivot :: RowIndex -> ColIndex -> Tableau Rational -> Tableau Rational #-}
{-# SPECIALIZE pivot :: RowIndex -> ColIndex -> Tableau Double -> Tableau Double #-}
pivot :: forall r.
(Fractional r, Eq r) =>
RowIndex -> RowIndex -> Tableau r -> Tableau r
pivot RowIndex
r RowIndex
s Tableau r
tbl =
    forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Tableau r -> Bool
isValidTableau Tableau r
tbl) forall a b. (a -> b) -> a -> b
$  -- precondition
    forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Tableau r -> Bool
isValidTableau Tableau r
tbl') forall a b. (a -> b) -> a -> b
$ -- postcondition
      Tableau r
tbl'
  where
    tbl' :: Tableau r
tbl' = forall a. RowIndex -> a -> IntMap a -> IntMap a
IM.insert RowIndex
s (IntMap r, r)
row_s forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map (IntMap r, r) -> (IntMap r, r)
f forall a b. (a -> b) -> a -> b
$ forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
r forall a b. (a -> b) -> a -> b
$ Tableau r
tbl
    f :: (IntMap r, r) -> (IntMap r, r)
f orig :: (IntMap r, r)
orig@(IntMap r
row_i, r
row_i_val) =
      case forall a. RowIndex -> IntMap a -> Maybe a
IM.lookup RowIndex
s IntMap r
row_i of
        Maybe r
Nothing -> (IntMap r, r)
orig
        Just r
c ->
          ( forall a. (a -> Bool) -> IntMap a -> IntMap a
IM.filter (r
0forall a. Eq a => a -> a -> Bool
/=) forall a b. (a -> b) -> a -> b
$ forall a. (a -> a -> a) -> IntMap a -> IntMap a -> IntMap a
IM.unionWith forall a. Num a => a -> a -> a
(+) (forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
s IntMap r
row_i) (forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map (forall a. Num a => a -> a
negate r
c forall a. Num a => a -> a -> a
*) IntMap r
row_r')
          , r
row_i_val forall a. Num a => a -> a -> a
- r
cforall a. Num a => a -> a -> a
*r
row_r_val'
          )
    (IntMap r
row_r, r
row_r_val) = forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
r Tableau r
tbl
    a_rs :: r
a_rs = IntMap r
row_r forall a. IntMap a -> RowIndex -> a
IM.! RowIndex
s
    row_r' :: IntMap r
row_r' = forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map (forall a. Fractional a => a -> a -> a
/ r
a_rs) forall a b. (a -> b) -> a -> b
$ forall a. RowIndex -> a -> IntMap a -> IntMap a
IM.insert RowIndex
r r
1 forall a b. (a -> b) -> a -> b
$ forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
s IntMap r
row_r
    row_r_val' :: r
row_r_val' = r
row_r_val forall a. Fractional a => a -> a -> a
/ r
a_rs
    row_s :: (IntMap r, r)
row_s = (IntMap r
row_r', r
row_r_val')

-- | Lookup a row by basic variable
lookupRow :: RowIndex -> Tableau r -> Row r
lookupRow :: forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
r Tableau r
m = Tableau r
m forall a. IntMap a -> RowIndex -> a
IM.! RowIndex
r

-- 行の基底変数の列が0になるように変形
normalizeRow :: (Num r, Eq r) => Tableau r -> Row r -> Row r
normalizeRow :: forall r. (Num r, Eq r) => Tableau r -> Row r -> Row r
normalizeRow Tableau r
a (VarMap r
row0,r
val0) = (VarMap r, r)
obj'
  where
    obj' :: (VarMap r, r)
obj' = forall {a} {b}. (Eq a, Num a) => (IntMap a, b) -> (IntMap a, b)
g forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' forall {a} {b}.
(Num a, Num b) =>
(IntMap a, b) -> (IntMap a, b) -> (IntMap a, b)
f (forall a. IntMap a
IM.empty, r
val0) forall a b. (a -> b) -> a -> b
$
           [ case forall a. RowIndex -> IntMap a -> Maybe a
IM.lookup RowIndex
j Tableau r
a of
               Maybe (VarMap r, r)
Nothing -> (forall a. RowIndex -> a -> IntMap a
IM.singleton RowIndex
j r
x, r
0)
               Just (VarMap r
row,r
val) -> (forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map ((-r
x)forall a. Num a => a -> a -> a
*) (forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
j VarMap r
row), -r
xforall a. Num a => a -> a -> a
*r
val)
           | (RowIndex
j,r
x) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList VarMap r
row0 ]
    f :: (IntMap a, b) -> (IntMap a, b) -> (IntMap a, b)
f (IntMap a
m1,b
v1) (IntMap a
m2,b
v2) = (forall a. (a -> a -> a) -> IntMap a -> IntMap a -> IntMap a
IM.unionWith forall a. Num a => a -> a -> a
(+) IntMap a
m1 IntMap a
m2, b
v1forall a. Num a => a -> a -> a
+b
v2)
    g :: (IntMap a, b) -> (IntMap a, b)
g (IntMap a
m,b
v) = (forall a. (a -> Bool) -> IntMap a -> IntMap a
IM.filter (a
0forall a. Eq a => a -> a -> Bool
/=) IntMap a
m, b
v)

setRow :: (Num r, Eq r) => Tableau r -> RowIndex -> Row r -> Tableau r
setRow :: forall r.
(Num r, Eq r) =>
Tableau r -> RowIndex -> Row r -> Tableau r
setRow Tableau r
tbl RowIndex
i Row r
row = forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Tableau r -> Bool
isValidTableau Tableau r
tbl) forall a b. (a -> b) -> a -> b
$ forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Tableau r -> Bool
isValidTableau Tableau r
tbl') forall a b. (a -> b) -> a -> b
$ Tableau r
tbl'
  where
    tbl' :: Tableau r
tbl' = forall a. RowIndex -> a -> IntMap a -> IntMap a
IM.insert RowIndex
i (forall r. (Num r, Eq r) => Tableau r -> Row r -> Row r
normalizeRow Tableau r
tbl Row r
row) Tableau r
tbl

addRow :: (Num r, Eq r) => Tableau r -> RowIndex -> Row r -> Tableau r
addRow :: forall r.
(Num r, Eq r) =>
Tableau r -> RowIndex -> Row r -> Tableau r
addRow Tableau r
tbl RowIndex
i Row r
row = forall a. (?callStack::CallStack) => Bool -> a -> a
assert (RowIndex
i forall a. RowIndex -> IntMap a -> Bool
`IM.notMember` Tableau r
tbl) forall a b. (a -> b) -> a -> b
$ forall r.
(Num r, Eq r) =>
Tableau r -> RowIndex -> Row r -> Tableau r
setRow Tableau r
tbl RowIndex
i Row r
row

setObjFun :: (Num r, Eq r) => Tableau r -> LA.Expr r -> Tableau r
setObjFun :: forall r. (Num r, Eq r) => Tableau r -> Expr r -> Tableau r
setObjFun Tableau r
tbl Expr r
e = forall r.
(Num r, Eq r) =>
Tableau r -> RowIndex -> Row r -> Tableau r
addRow Tableau r
tbl RowIndex
objRowIndex (IntMap r, r)
row
  where
    row :: (IntMap r, r)
row =
      case forall r. Num r => RowIndex -> Expr r -> (r, Expr r)
LA.extract RowIndex
LA.unitVar Expr r
e of
        (r
c, Expr r
e') -> (forall r. Expr r -> IntMap r
LA.coeffMap (forall v. AdditiveGroup v => v -> v
negateV Expr r
e'), r
c)

copyObjRow :: (Num r, Eq r) => Tableau r -> Tableau r -> Tableau r
copyObjRow :: forall r. (Num r, Eq r) => Tableau r -> Tableau r -> Tableau r
copyObjRow Tableau r
from Tableau r
to =
  case forall a. RowIndex -> IntMap a -> Maybe a
IM.lookup RowIndex
objRowIndex Tableau r
from of
    Maybe (Row r)
Nothing -> forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
objRowIndex Tableau r
to
    Just Row r
row -> forall r.
(Num r, Eq r) =>
Tableau r -> RowIndex -> Row r -> Tableau r
addRow Tableau r
to RowIndex
objRowIndex Row r
row

currentValue :: Num r => Tableau r -> Var -> r
currentValue :: forall r. Num r => Tableau r -> RowIndex -> r
currentValue Tableau r
tbl RowIndex
v =
  case forall a. RowIndex -> IntMap a -> Maybe a
IM.lookup RowIndex
v Tableau r
tbl of
    Maybe (Row r)
Nothing -> r
0
    Just (VarMap r
_, r
val) -> r
val

currentObjValue :: Tableau r -> r
currentObjValue :: forall r. Tableau r -> r
currentObjValue = forall a b. (a, b) -> b
snd forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
objRowIndex

isValidTableau :: Tableau r -> Bool
isValidTableau :: forall r. Tableau r -> Bool
isValidTableau Tableau r
tbl =
  forall (t :: * -> *). Foldable t => t Bool -> Bool
and [RowIndex
v forall a. RowIndex -> IntMap a -> Bool
`IM.notMember` VarMap r
m | (RowIndex
v, (VarMap r
m,r
_)) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
v forall a. Eq a => a -> a -> Bool
/= RowIndex
objRowIndex] Bool -> Bool -> Bool
&&
  forall (t :: * -> *). Foldable t => t Bool -> Bool
and [IntSet -> Bool
IS.null (forall a. IntMap a -> IntSet
IM.keysSet VarMap r
m IntSet -> IntSet -> IntSet
`IS.intersection` IntSet
vs) | (VarMap r
m,r
_) <- forall a. IntMap a -> [a]
IM.elems Tableau r
tbl']
  where
    tbl' :: Tableau r
tbl' = forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
objRowIndex Tableau r
tbl
    vs :: IntSet
vs = forall a. IntMap a -> IntSet
IM.keysSet Tableau r
tbl'

isFeasible :: Real r => Tableau r -> Bool
isFeasible :: forall r. Real r => Tableau r -> Bool
isFeasible Tableau r
tbl =
  forall (t :: * -> *). Foldable t => t Bool -> Bool
and [r
val forall a. Ord a => a -> a -> Bool
>= r
0 | (RowIndex
v, (VarMap r
_,r
val)) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
v forall a. Eq a => a -> a -> Bool
/= RowIndex
objRowIndex]

isOptimal :: Real r => OptDir -> Tableau r -> Bool
isOptimal :: forall r. Real r => OptDir -> Tableau r -> Bool
isOptimal OptDir
optdir Tableau r
tbl =
  forall (t :: * -> *). Foldable t => t Bool -> Bool
and [Bool -> Bool
not (r -> Bool
cmp r
cj) | r
cj <- forall a. IntMap a -> [a]
IM.elems (forall a b. (a, b) -> a
fst (forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
objRowIndex Tableau r
tbl))]
  where
    cmp :: r -> Bool
cmp = case OptDir
optdir of
            OptDir
OptMin -> (r
0forall a. Ord a => a -> a -> Bool
<)
            OptDir
OptMax -> (r
0forall a. Ord a => a -> a -> Bool
>)

isImproving :: Real r => OptDir -> Tableau r -> Tableau r -> Bool
isImproving :: forall r. Real r => OptDir -> Tableau r -> Tableau r -> Bool
isImproving OptDir
OptMin Tableau r
from Tableau r
to = forall r. Tableau r -> r
currentObjValue Tableau r
to forall a. Ord a => a -> a -> Bool
<= forall r. Tableau r -> r
currentObjValue Tableau r
from
isImproving OptDir
OptMax Tableau r
from Tableau r
to = forall r. Tableau r -> r
currentObjValue Tableau r
to forall a. Ord a => a -> a -> Bool
>= forall r. Tableau r -> r
currentObjValue Tableau r
from

-- ---------------------------------------------------------------------------
-- primal simplex

simplex :: (Real r, Fractional r) => OptDir -> Tableau r -> (Bool, Tableau r)
{-# SPECIALIZE simplex :: OptDir -> Tableau Rational -> (Bool, Tableau Rational) #-}
{-# SPECIALIZE simplex :: OptDir -> Tableau Double -> (Bool, Tableau Double) #-}
simplex :: forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> (Bool, Tableau r)
simplex OptDir
optdir = forall {r}.
(Real r, Fractional r) =>
Tableau r -> (Bool, Tableau r)
go
  where
    go :: Tableau r -> (Bool, Tableau r)
go Tableau r
tbl = forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Real r => Tableau r -> Bool
isFeasible Tableau r
tbl) forall a b. (a -> b) -> a -> b
$
      case forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> PivotResult r
primalPivot OptDir
optdir Tableau r
tbl of
        PivotResult r
PivotFinished  -> forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Real r => OptDir -> Tableau r -> Bool
isOptimal OptDir
optdir Tableau r
tbl) (Bool
True, Tableau r
tbl)
        PivotResult r
PivotUnbounded -> (Bool
False, Tableau r
tbl)
        PivotSuccess Tableau r
tbl' -> forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Real r => OptDir -> Tableau r -> Tableau r -> Bool
isImproving OptDir
optdir Tableau r
tbl Tableau r
tbl') forall a b. (a -> b) -> a -> b
$ Tableau r -> (Bool, Tableau r)
go Tableau r
tbl'

primalPivot :: (Real r, Fractional r) => OptDir -> Tableau r -> PivotResult r
{-# INLINE primalPivot #-}
primalPivot :: forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> PivotResult r
primalPivot OptDir
optdir Tableau r
tbl
  | forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, r)]
cs   = forall r. PivotResult r
PivotFinished
  | forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, r)]
rs   = forall r. PivotResult r
PivotUnbounded
  | Bool
otherwise = forall r. Tableau r -> PivotResult r
PivotSuccess (forall r.
(Fractional r, Eq r) =>
RowIndex -> RowIndex -> Tableau r -> Tableau r
pivot RowIndex
r RowIndex
s Tableau r
tbl)
  where
    cmp :: r -> Bool
cmp = case OptDir
optdir of
            OptDir
OptMin -> (r
0forall a. Ord a => a -> a -> Bool
<)
            OptDir
OptMax -> (r
0forall a. Ord a => a -> a -> Bool
>)
    cs :: [(RowIndex, r)]
cs = [(RowIndex
j,r
cj) | (RowIndex
j,r
cj) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList (forall a b. (a, b) -> a
fst (forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
objRowIndex Tableau r
tbl)), r -> Bool
cmp r
cj]
    -- smallest subscript rule
    s :: RowIndex
s = forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall a. [a] -> a
head [(RowIndex, r)]
cs
    -- classical rule
    --s = fst $ (if optdir==OptMin then maximumBy else minimumBy) (compare `on` snd) cs
    rs :: [(RowIndex, r)]
rs = [ (RowIndex
i, r
y_i0 forall a. Fractional a => a -> a -> a
/ r
y_is)
         | (RowIndex
i, (VarMap r
row_i, r
y_i0)) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
i forall a. Eq a => a -> a -> Bool
/= RowIndex
objRowIndex
         , let y_is :: r
y_is = forall a. a -> RowIndex -> IntMap a -> a
IM.findWithDefault r
0 RowIndex
s VarMap r
row_i, r
y_is forall a. Ord a => a -> a -> Bool
> r
0
         ]
    r :: RowIndex
r = forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy (forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing forall a b. (a, b) -> b
snd) [(RowIndex, r)]
rs

-- ---------------------------------------------------------------------------
-- dual simplex

dualSimplex :: (Real r, Fractional r) => OptDir -> Tableau r -> (Bool, Tableau r)
{-# SPECIALIZE dualSimplex :: OptDir -> Tableau Rational -> (Bool, Tableau Rational) #-}
{-# SPECIALIZE dualSimplex :: OptDir -> Tableau Double -> (Bool, Tableau Double) #-}
dualSimplex :: forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> (Bool, Tableau r)
dualSimplex OptDir
optdir = forall {r}.
(Real r, Fractional r) =>
Tableau r -> (Bool, Tableau r)
go
  where
    go :: Tableau r -> (Bool, Tableau r)
go Tableau r
tbl = forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Real r => OptDir -> Tableau r -> Bool
isOptimal OptDir
optdir Tableau r
tbl) forall a b. (a -> b) -> a -> b
$
      case forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> PivotResult r
dualPivot OptDir
optdir Tableau r
tbl of
        PivotResult r
PivotFinished  -> forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Real r => Tableau r -> Bool
isFeasible Tableau r
tbl) forall a b. (a -> b) -> a -> b
$ (Bool
True, Tableau r
tbl)
        PivotResult r
PivotUnbounded -> (Bool
False, Tableau r
tbl)
        PivotSuccess Tableau r
tbl' -> forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Real r => OptDir -> Tableau r -> Tableau r -> Bool
isImproving OptDir
optdir Tableau r
tbl' Tableau r
tbl) forall a b. (a -> b) -> a -> b
$ Tableau r -> (Bool, Tableau r)
go Tableau r
tbl'

dualPivot :: (Real r, Fractional r) => OptDir -> Tableau r -> PivotResult r
{-# INLINE dualPivot #-}
dualPivot :: forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> PivotResult r
dualPivot OptDir
optdir Tableau r
tbl
  | forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, VarMap r)]
rs   = forall r. PivotResult r
PivotFinished
  | forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, r)]
cs   = forall r. PivotResult r
PivotUnbounded
  | Bool
otherwise = forall r. Tableau r -> PivotResult r
PivotSuccess (forall r.
(Fractional r, Eq r) =>
RowIndex -> RowIndex -> Tableau r -> Tableau r
pivot RowIndex
r RowIndex
s Tableau r
tbl)
  where
    rs :: [(RowIndex, VarMap r)]
rs = [(RowIndex
i, VarMap r
row_i) | (RowIndex
i, (VarMap r
row_i, r
y_i0)) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
i forall a. Eq a => a -> a -> Bool
/= RowIndex
objRowIndex, r
0 forall a. Ord a => a -> a -> Bool
> r
y_i0]
    (RowIndex
r, VarMap r
row_r) = forall a. [a] -> a
head [(RowIndex, VarMap r)]
rs
    cs :: [(RowIndex, r)]
cs = [ (RowIndex
j, if OptDir
optdirforall a. Eq a => a -> a -> Bool
==OptDir
OptMin then r
y_0j forall a. Fractional a => a -> a -> a
/ r
y_rj else - r
y_0j forall a. Fractional a => a -> a -> a
/ r
y_rj)
         | (RowIndex
j, r
y_rj) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList VarMap r
row_r
         , r
y_rj forall a. Ord a => a -> a -> Bool
< r
0
         , let y_0j :: r
y_0j = forall a. a -> RowIndex -> IntMap a -> a
IM.findWithDefault r
0 RowIndex
j VarMap r
obj
         ]
    s :: RowIndex
s = forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy (forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing forall a b. (a, b) -> b
snd) [(RowIndex, r)]
cs
    (VarMap r
obj,r
_) = forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
objRowIndex Tableau r
tbl

-- ---------------------------------------------------------------------------
-- phase I of the two-phased method

phaseI :: (Real r, Fractional r) => Tableau r -> VarSet -> (Bool, Tableau r)
{-# SPECIALIZE phaseI :: Tableau Rational -> VarSet -> (Bool, Tableau Rational) #-}
{-# SPECIALIZE phaseI :: Tableau Double -> VarSet -> (Bool, Tableau Double) #-}
phaseI :: forall r.
(Real r, Fractional r) =>
Tableau r -> IntSet -> (Bool, Tableau r)
phaseI Tableau r
tbl IntSet
avs
  | forall r. Tableau r -> r
currentObjValue Tableau r
tbl1' forall a. Eq a => a -> a -> Bool
/= r
0 = (Bool
False, Tableau r
tbl1')
  | Bool
otherwise = (Bool
True, forall r. (Num r, Eq r) => Tableau r -> Tableau r -> Tableau r
copyObjRow Tableau r
tbl forall a b. (a -> b) -> a -> b
$ forall r.
(Real r, Fractional r) =>
IntSet -> Tableau r -> Tableau r
removeArtificialVariables IntSet
avs forall a b. (a -> b) -> a -> b
$ Tableau r
tbl1')
  where
    optdir :: OptDir
optdir = OptDir
OptMax
    tbl1 :: Tableau r
tbl1 = forall r. (Num r, Eq r) => Tableau r -> Expr r -> Tableau r
setObjFun Tableau r
tbl forall a b. (a -> b) -> a -> b
$ forall v. AdditiveGroup v => v -> v
negateV forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV [forall r. Num r => RowIndex -> Expr r
LA.var RowIndex
v | RowIndex
v <- IntSet -> [RowIndex]
IS.toList IntSet
avs]
    tbl1' :: Tableau r
tbl1' = forall {r}. (Real r, Fractional r) => Tableau r -> Tableau r
go Tableau r
tbl1
    go :: Tableau r -> Tableau r
go Tableau r
tbl2
      | forall r. Tableau r -> r
currentObjValue Tableau r
tbl2 forall a. Eq a => a -> a -> Bool
== r
0 = Tableau r
tbl2
      | Bool
otherwise =
        case forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> PivotResult r
primalPivot OptDir
optdir Tableau r
tbl2 of
          PivotSuccess Tableau r
tbl2' -> forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Real r => OptDir -> Tableau r -> Tableau r -> Bool
isImproving OptDir
optdir Tableau r
tbl2 Tableau r
tbl2') forall a b. (a -> b) -> a -> b
$ Tableau r -> Tableau r
go Tableau r
tbl2'
          PivotResult r
PivotFinished -> forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Real r => OptDir -> Tableau r -> Bool
isOptimal OptDir
optdir Tableau r
tbl2) Tableau r
tbl2
          PivotResult r
PivotUnbounded -> forall a. (?callStack::CallStack) => String -> a
error String
"phaseI: should not happen"

-- post-processing of phaseI
-- pivotを使ってartificial variablesを基底から除いて、削除
removeArtificialVariables :: (Real r, Fractional r) => VarSet -> Tableau r -> Tableau r
removeArtificialVariables :: forall r.
(Real r, Fractional r) =>
IntSet -> Tableau r -> Tableau r
removeArtificialVariables IntSet
avs Tableau r
tbl0 = Tableau r
tbl2
  where
    tbl1 :: Tableau r
tbl1 = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' forall {b}.
(Eq b, Fractional b) =>
IntMap (IntMap b, b) -> RowIndex -> IntMap (IntMap b, b)
f (forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
objRowIndex Tableau r
tbl0) (IntSet -> [RowIndex]
IS.toList IntSet
avs)
    tbl2 :: Tableau r
tbl2 = forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map (\(IntMap r
row,r
val) ->  (forall a. (RowIndex -> a -> Bool) -> IntMap a -> IntMap a
IM.filterWithKey (\RowIndex
j r
_ -> RowIndex
j RowIndex -> IntSet -> Bool
`IS.notMember` IntSet
avs) IntMap r
row, r
val)) Tableau r
tbl1
    f :: IntMap (IntMap b, b) -> RowIndex -> IntMap (IntMap b, b)
f IntMap (IntMap b, b)
tbl RowIndex
i =
      case forall a. RowIndex -> IntMap a -> Maybe a
IM.lookup RowIndex
i IntMap (IntMap b, b)
tbl of
        Maybe (IntMap b, b)
Nothing -> IntMap (IntMap b, b)
tbl
        Just (IntMap b, b)
row ->
          case [RowIndex
j | (RowIndex
j,b
c) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList (forall a b. (a, b) -> a
fst (IntMap b, b)
row), b
c forall a. Eq a => a -> a -> Bool
/= b
0, RowIndex
j RowIndex -> IntSet -> Bool
`IS.notMember` IntSet
avs] of
            [] -> forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
i IntMap (IntMap b, b)
tbl
            RowIndex
j:[RowIndex]
_ -> forall r.
(Fractional r, Eq r) =>
RowIndex -> RowIndex -> Tableau r -> Tableau r
pivot RowIndex
i RowIndex
j IntMap (IntMap b, b)
tbl

-- ---------------------------------------------------------------------------
-- primal-dual simplex

data PDResult = PDUnsat | PDOptimal | PDUnbounded

primalDualSimplex :: (Real r, Fractional r) => OptDir -> Tableau r -> (Bool, Tableau r)
{-# SPECIALIZE primalDualSimplex :: OptDir -> Tableau Rational -> (Bool, Tableau Rational) #-}
{-# SPECIALIZE primalDualSimplex :: OptDir -> Tableau Double -> (Bool, Tableau Double) #-}
primalDualSimplex :: forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> (Bool, Tableau r)
primalDualSimplex OptDir
optdir = forall {r}.
(Real r, Fractional r) =>
Tableau r -> (Bool, Tableau r)
go
  where
    go :: Tableau r -> (Bool, Tableau r)
go Tableau r
tbl =
      case forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> Either PDResult (Tableau r)
pdPivot OptDir
optdir Tableau r
tbl of
        Left PDResult
PDOptimal   -> forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Real r => Tableau r -> Bool
isFeasible Tableau r
tbl) forall a b. (a -> b) -> a -> b
$ forall a. (?callStack::CallStack) => Bool -> a -> a
assert (forall r. Real r => OptDir -> Tableau r -> Bool
isOptimal OptDir
optdir Tableau r
tbl) forall a b. (a -> b) -> a -> b
$ (Bool
True, Tableau r
tbl)
        Left PDResult
PDUnsat     -> forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Bool -> Bool
not (forall r. Real r => Tableau r -> Bool
isFeasible Tableau r
tbl)) forall a b. (a -> b) -> a -> b
$ (Bool
False, Tableau r
tbl)
        Left PDResult
PDUnbounded -> forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Bool -> Bool
not (forall r. Real r => OptDir -> Tableau r -> Bool
isOptimal OptDir
optdir Tableau r
tbl)) forall a b. (a -> b) -> a -> b
$ (Bool
False, Tableau r
tbl)
        Right Tableau r
tbl' -> Tableau r -> (Bool, Tableau r)
go Tableau r
tbl'

pdPivot :: (Real r, Fractional r) => OptDir -> Tableau r -> Either PDResult (Tableau r)
pdPivot :: forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> Either PDResult (Tableau r)
pdPivot OptDir
optdir Tableau r
tbl
  | forall (t :: * -> *) a. Foldable t => t a -> Bool
null forall {b}. [(Either RowIndex b, r)]
ps Bool -> Bool -> Bool
&& forall (t :: * -> *) a. Foldable t => t a -> Bool
null forall {a}. [(Either a RowIndex, r)]
qs = forall a b. a -> Either a b
Left PDResult
PDOptimal -- optimal
  | Bool
otherwise =
      case Either RowIndex RowIndex
ret of
        Left RowIndex
p -> -- primal update
          let rs :: [(RowIndex, r)]
rs = [ (RowIndex
i, (r
bi forall a. Num a => a -> a -> a
- r
t) forall a. Fractional a => a -> a -> a
/ r
y_ip)
                   | (RowIndex
i, (VarMap r
row_i, r
bi)) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
i forall a. Eq a => a -> a -> Bool
/= RowIndex
objRowIndex
                   , let y_ip :: r
y_ip = forall a. a -> RowIndex -> IntMap a -> a
IM.findWithDefault r
0 RowIndex
p VarMap r
row_i, r
y_ip forall a. Ord a => a -> a -> Bool
> r
0
                   ]
              q :: RowIndex
q = forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy (forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing forall a b. (a, b) -> b
snd) [(RowIndex, r)]
rs
          in if forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, r)]
rs
             then forall a b. a -> Either a b
Left PDResult
PDUnsat
             else forall a b. b -> Either a b
Right (forall r.
(Fractional r, Eq r) =>
RowIndex -> RowIndex -> Tableau r -> Tableau r
pivot RowIndex
q RowIndex
p Tableau r
tbl)
        Right RowIndex
q -> -- dual update
          let (VarMap r
row_q, r
_bq) = (Tableau r
tbl forall a. IntMap a -> RowIndex -> a
IM.! RowIndex
q)
              cs :: [(RowIndex, r)]
cs = [ (RowIndex
j, (r
cj'forall a. Num a => a -> a -> a
-r
t) forall a. Fractional a => a -> a -> a
/ (-r
y_qj))
                   | (RowIndex
j, r
y_qj) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList VarMap r
row_q
                   , r
y_qj forall a. Ord a => a -> a -> Bool
< r
0
                   , let cj :: r
cj  = forall a. a -> RowIndex -> IntMap a -> a
IM.findWithDefault r
0 RowIndex
j VarMap r
obj
                   , let cj' :: r
cj' = if OptDir
optdirforall a. Eq a => a -> a -> Bool
==OptDir
OptMax then r
cj else -r
cj
                   ]
              p :: RowIndex
p = forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy (forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing forall a b. (a, b) -> b
snd) [(RowIndex, r)]
cs
              (VarMap r
obj,r
_) = forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
objRowIndex Tableau r
tbl
          in if forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, r)]
cs
             then forall a b. a -> Either a b
Left PDResult
PDUnbounded -- dual infeasible
             else forall a b. b -> Either a b
Right (forall r.
(Fractional r, Eq r) =>
RowIndex -> RowIndex -> Tableau r -> Tableau r
pivot RowIndex
q RowIndex
p Tableau r
tbl)
  where
    qs :: [(Either a RowIndex, r)]
qs = [ (forall a b. b -> Either a b
Right RowIndex
i, r
bi) | (RowIndex
i, (VarMap r
_row_i, r
bi)) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
i forall a. Eq a => a -> a -> Bool
/= RowIndex
objRowIndex, r
0 forall a. Ord a => a -> a -> Bool
> r
bi ]
    ps :: [(Either RowIndex b, r)]
ps = [ (forall a b. a -> Either a b
Left RowIndex
j, r
cj')
         | (RowIndex
j,r
cj) <- forall a. IntMap a -> [(RowIndex, a)]
IM.toList (forall a b. (a, b) -> a
fst (forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
objRowIndex Tableau r
tbl))
         , let cj' :: r
cj' = if OptDir
optdirforall a. Eq a => a -> a -> Bool
==OptDir
OptMax then r
cj else -r
cj
         , r
0 forall a. Ord a => a -> a -> Bool
> r
cj' ]
    (Either RowIndex RowIndex
ret, r
t) = forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy (forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing forall a b. (a, b) -> b
snd) (forall {a}. [(Either a RowIndex, r)]
qs forall a. [a] -> [a] -> [a]
++ forall {b}. [(Either RowIndex b, r)]
ps)

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

toCSV :: (Num r) => (r -> String) -> Tableau r -> String
toCSV :: forall r. Num r => (r -> String) -> Tableau r -> String
toCSV r -> String
showCell Tableau r
tbl = [String] -> String
unlines forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map (forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. a -> [a] -> [a]
intersperse String
",") forall a b. (a -> b) -> a -> b
$ [String]
header forall a. a -> [a] -> [a]
: [[String]]
body
  where
    header :: [String]
    header :: [String]
header = String
"" forall a. a -> [a] -> [a]
: forall a b. (a -> b) -> [a] -> [b]
map RowIndex -> String
colName [RowIndex]
cols forall a. [a] -> [a] -> [a]
++ [String
""]

    body :: [[String]]
    body :: [[String]]
body = [RowIndex -> (IntMap r, r) -> [String]
showRow RowIndex
i (forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
i Tableau r
tbl) | RowIndex
i <- [RowIndex]
rows]

    rows :: [RowIndex]
    rows :: [RowIndex]
rows = forall a. IntMap a -> [RowIndex]
IM.keys (forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
objRowIndex Tableau r
tbl) forall a. [a] -> [a] -> [a]
++ [RowIndex
objRowIndex]

    cols :: [ColIndex]
    cols :: [RowIndex]
cols = [RowIndex
0..RowIndex
colMax]
      where
        colMax :: RowIndex
colMax = forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (-RowIndex
1 forall a. a -> [a] -> [a]
: [RowIndex
c | (IntMap r
row, r
_) <- forall a. IntMap a -> [a]
IM.elems Tableau r
tbl, RowIndex
c <- forall a. IntMap a -> [RowIndex]
IM.keys IntMap r
row])

    rowName :: RowIndex -> String
    rowName :: RowIndex -> String
rowName RowIndex
i = if RowIndex
iforall a. Eq a => a -> a -> Bool
==RowIndex
objRowIndex then String
"obj" else String
"x" forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show RowIndex
i

    colName :: ColIndex -> String
    colName :: RowIndex -> String
colName RowIndex
j = String
"x" forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show RowIndex
j

    showRow :: RowIndex -> (IntMap r, r) -> [String]
showRow RowIndex
i (IntMap r
row, r
row_val) = RowIndex -> String
rowName RowIndex
i forall a. a -> [a] -> [a]
: [r -> String
showCell (forall a. a -> RowIndex -> IntMap a -> a
IM.findWithDefault r
0 RowIndex
j IntMap r
row') | RowIndex
j <- [RowIndex]
cols] forall a. [a] -> [a] -> [a]
++ [r -> String
showCell r
row_val]
      where row' :: IntMap r
row' = forall a. RowIndex -> a -> IntMap a -> IntMap a
IM.insert RowIndex
i r
1 IntMap r
row

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