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

emptyTableau :: Tableau r
emptyTableau :: forall r. Tableau r
emptyTableau = IntMap (Row r)
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 =
    Bool -> Tableau r -> Tableau r
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Tableau r -> Bool
forall r. Tableau r -> Bool
isValidTableau Tableau r
tbl) (Tableau r -> Tableau r) -> Tableau r -> Tableau r
forall a b. (a -> b) -> a -> b
$  -- precondition
    Bool -> Tableau r -> Tableau r
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Tableau r -> Bool
forall r. Tableau r -> Bool
isValidTableau Tableau r
tbl') (Tableau r -> Tableau r) -> Tableau r -> Tableau r
forall a b. (a -> b) -> a -> b
$ -- postcondition
      Tableau r
tbl'
  where
    tbl' :: Tableau r
tbl' = RowIndex -> (IntMap r, r) -> Tableau r -> Tableau r
forall a. RowIndex -> a -> IntMap a -> IntMap a
IM.insert RowIndex
s (IntMap r, r)
row_s (Tableau r -> Tableau r) -> Tableau r -> Tableau r
forall a b. (a -> b) -> a -> b
$ ((IntMap r, r) -> (IntMap r, r)) -> Tableau r -> Tableau r
forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map (IntMap r, r) -> (IntMap r, r)
f (Tableau r -> Tableau r) -> Tableau r -> Tableau r
forall a b. (a -> b) -> a -> b
$ RowIndex -> Tableau r -> Tableau r
forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
r (Tableau r -> Tableau r) -> Tableau r -> Tableau 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 RowIndex -> IntMap r -> Maybe r
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 ->
          ( (r -> Bool) -> IntMap r -> IntMap r
forall a. (a -> Bool) -> IntMap a -> IntMap a
IM.filter (r
0r -> r -> Bool
forall a. Eq a => a -> a -> Bool
/=) (IntMap r -> IntMap r) -> IntMap r -> IntMap r
forall a b. (a -> b) -> a -> b
$ (r -> r -> r) -> IntMap r -> IntMap r -> IntMap r
forall a. (a -> a -> a) -> IntMap a -> IntMap a -> IntMap a
IM.unionWith r -> r -> r
forall a. Num a => a -> a -> a
(+) (RowIndex -> IntMap r -> IntMap r
forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
s IntMap r
row_i) ((r -> r) -> IntMap r -> IntMap r
forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map (r -> r
forall a. Num a => a -> a
negate r
c r -> r -> r
forall a. Num a => a -> a -> a
*) IntMap r
row_r')
          , r
row_i_val r -> r -> r
forall a. Num a => a -> a -> a
- r
cr -> r -> r
forall a. Num a => a -> a -> a
*r
row_r_val'
          )
    (IntMap r
row_r, r
row_r_val) = RowIndex -> Tableau r -> (IntMap r, r)
forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
r Tableau r
tbl
    a_rs :: r
a_rs = IntMap r
row_r IntMap r -> RowIndex -> r
forall a. IntMap a -> RowIndex -> a
IM.! RowIndex
s
    row_r' :: IntMap r
row_r' = (r -> r) -> IntMap r -> IntMap r
forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map (r -> r -> r
forall a. Fractional a => a -> a -> a
/ r
a_rs) (IntMap r -> IntMap r) -> IntMap r -> IntMap r
forall a b. (a -> b) -> a -> b
$ RowIndex -> r -> IntMap r -> IntMap r
forall a. RowIndex -> a -> IntMap a -> IntMap a
IM.insert RowIndex
r r
1 (IntMap r -> IntMap r) -> IntMap r -> IntMap r
forall a b. (a -> b) -> a -> b
$ RowIndex -> IntMap r -> IntMap r
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 r -> r -> r
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 Tableau r -> RowIndex -> Row r
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' = (VarMap r, r) -> (VarMap r, r)
forall {a} {b}. (Eq a, Num a) => (IntMap a, b) -> (IntMap a, b)
g ((VarMap r, r) -> (VarMap r, r)) -> (VarMap r, r) -> (VarMap r, r)
forall a b. (a -> b) -> a -> b
$ ((VarMap r, r) -> (VarMap r, r) -> (VarMap r, r))
-> (VarMap r, r) -> [(VarMap r, r)] -> (VarMap r, r)
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (VarMap r, r) -> (VarMap r, r) -> (VarMap r, r)
forall {a} {b}.
(Num a, Num b) =>
(IntMap a, b) -> (IntMap a, b) -> (IntMap a, b)
f (VarMap r
forall a. IntMap a
IM.empty, r
val0) ([(VarMap r, r)] -> (VarMap r, r))
-> [(VarMap r, r)] -> (VarMap r, r)
forall a b. (a -> b) -> a -> b
$
           [ case RowIndex -> Tableau r -> Maybe (VarMap r, r)
forall a. RowIndex -> IntMap a -> Maybe a
IM.lookup RowIndex
j Tableau r
a of
               Maybe (VarMap r, r)
Nothing -> (RowIndex -> r -> VarMap r
forall a. RowIndex -> a -> IntMap a
IM.singleton RowIndex
j r
x, r
0)
               Just (VarMap r
row,r
val) -> ((r -> r) -> VarMap r -> VarMap r
forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map ((-r
x)r -> r -> r
forall a. Num a => a -> a -> a
*) (RowIndex -> VarMap r -> VarMap r
forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
j VarMap r
row), -r
xr -> r -> r
forall a. Num a => a -> a -> a
*r
val)
           | (RowIndex
j,r
x) <- VarMap r -> [(RowIndex, r)]
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) = ((a -> a -> a) -> IntMap a -> IntMap a -> IntMap a
forall a. (a -> a -> a) -> IntMap a -> IntMap a -> IntMap a
IM.unionWith a -> a -> a
forall a. Num a => a -> a -> a
(+) IntMap a
m1 IntMap a
m2, b
v1b -> b -> b
forall a. Num a => a -> a -> a
+b
v2)
    g :: (IntMap a, b) -> (IntMap a, b)
g (IntMap a
m,b
v) = ((a -> Bool) -> IntMap a -> IntMap a
forall a. (a -> Bool) -> IntMap a -> IntMap a
IM.filter (a
0a -> a -> Bool
forall 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 = Bool -> Tableau r -> Tableau r
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Tableau r -> Bool
forall r. Tableau r -> Bool
isValidTableau Tableau r
tbl) (Tableau r -> Tableau r) -> Tableau r -> Tableau r
forall a b. (a -> b) -> a -> b
$ Bool -> Tableau r -> Tableau r
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Tableau r -> Bool
forall r. Tableau r -> Bool
isValidTableau Tableau r
tbl') (Tableau r -> Tableau r) -> Tableau r -> Tableau r
forall a b. (a -> b) -> a -> b
$ Tableau r
tbl'
  where
    tbl' :: Tableau r
tbl' = RowIndex -> Row r -> Tableau r -> Tableau r
forall a. RowIndex -> a -> IntMap a -> IntMap a
IM.insert RowIndex
i (Tableau r -> Row r -> Row r
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 = Bool -> Tableau r -> Tableau r
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (RowIndex
i RowIndex -> Tableau r -> Bool
forall a. RowIndex -> IntMap a -> Bool
`IM.notMember` Tableau r
tbl) (Tableau r -> Tableau r) -> Tableau r -> Tableau r
forall a b. (a -> b) -> a -> b
$ Tableau r -> RowIndex -> Row r -> Tableau r
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 = Tableau r -> RowIndex -> Row r -> Tableau r
forall r.
(Num r, Eq r) =>
Tableau r -> RowIndex -> Row r -> Tableau r
addRow Tableau r
tbl RowIndex
objRowIndex Row r
row
  where
    row :: Row r
row =
      case RowIndex -> Expr r -> (r, Expr r)
forall r. Num r => RowIndex -> Expr r -> (r, Expr r)
LA.extract RowIndex
LA.unitVar Expr r
e of
        (r
c, Expr r
e') -> (Expr r -> IntMap r
forall r. Expr r -> IntMap r
LA.coeffMap (Expr r -> Expr r
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 RowIndex -> Tableau r -> Maybe (Row r)
forall a. RowIndex -> IntMap a -> Maybe a
IM.lookup RowIndex
objRowIndex Tableau r
from of
    Maybe (Row r)
Nothing -> RowIndex -> Tableau r -> Tableau r
forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
objRowIndex Tableau r
to
    Just Row r
row -> Tableau r -> RowIndex -> Row r -> Tableau r
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 RowIndex -> Tableau r -> Maybe (Row r)
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 = (VarMap r, r) -> r
forall a b. (a, b) -> b
snd ((VarMap r, r) -> r)
-> (Tableau r -> (VarMap r, r)) -> Tableau r -> r
forall b c a. (b -> c) -> (a -> b) -> a -> c
. RowIndex -> Tableau r -> (VarMap r, r)
forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
objRowIndex

isValidTableau :: Tableau r -> Bool
isValidTableau :: forall r. Tableau r -> Bool
isValidTableau Tableau r
tbl =
  [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and [RowIndex
v RowIndex -> IntMap r -> Bool
forall a. RowIndex -> IntMap a -> Bool
`IM.notMember` IntMap r
m | (RowIndex
v, (IntMap r
m,r
_)) <- Tableau r -> [(RowIndex, Row r)]
forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
v RowIndex -> RowIndex -> Bool
forall a. Eq a => a -> a -> Bool
/= RowIndex
objRowIndex] Bool -> Bool -> Bool
&&
  [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and [IntSet -> Bool
IS.null (IntMap r -> IntSet
forall a. IntMap a -> IntSet
IM.keysSet IntMap r
m IntSet -> IntSet -> IntSet
`IS.intersection` IntSet
vs) | (IntMap r
m,r
_) <- Tableau r -> [Row r]
forall a. IntMap a -> [a]
IM.elems Tableau r
tbl']
  where
    tbl' :: Tableau r
tbl' = RowIndex -> Tableau r -> Tableau r
forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
objRowIndex Tableau r
tbl
    vs :: IntSet
vs = Tableau r -> IntSet
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 =
  [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and [r
val r -> r -> Bool
forall a. Ord a => a -> a -> Bool
>= r
0 | (RowIndex
v, (VarMap r
_,r
val)) <- Tableau r -> [(RowIndex, Row r)]
forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
v RowIndex -> RowIndex -> Bool
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 =
  [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and [Bool -> Bool
not (r -> Bool
cmp r
cj) | r
cj <- IntMap r -> [r]
forall a. IntMap a -> [a]
IM.elems ((IntMap r, r) -> IntMap r
forall a b. (a, b) -> a
fst (RowIndex -> Tableau r -> (IntMap r, r)
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
0r -> r -> Bool
forall a. Ord a => a -> a -> Bool
<)
            OptDir
OptMax -> (r
0r -> r -> Bool
forall 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 = Tableau r -> r
forall r. Tableau r -> r
currentObjValue Tableau r
to r -> r -> Bool
forall a. Ord a => a -> a -> Bool
<= Tableau r -> r
forall r. Tableau r -> r
currentObjValue Tableau r
from
isImproving OptDir
OptMax Tableau r
from Tableau r
to = Tableau r -> r
forall r. Tableau r -> r
currentObjValue Tableau r
to r -> r -> Bool
forall a. Ord a => a -> a -> Bool
>= Tableau r -> r
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 = Tableau r -> (Bool, Tableau r)
forall {r}.
(Real r, Fractional r) =>
Tableau r -> (Bool, Tableau r)
go
  where
    go :: Tableau r -> (Bool, Tableau r)
go Tableau r
tbl = Bool -> (Bool, Tableau r) -> (Bool, Tableau r)
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Tableau r -> Bool
forall r. Real r => Tableau r -> Bool
isFeasible Tableau r
tbl) ((Bool, Tableau r) -> (Bool, Tableau r))
-> (Bool, Tableau r) -> (Bool, Tableau r)
forall a b. (a -> b) -> a -> b
$
      case OptDir -> Tableau r -> PivotResult r
forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> PivotResult r
primalPivot OptDir
optdir Tableau r
tbl of
        PivotResult r
PivotFinished  -> Bool -> (Bool, Tableau r) -> (Bool, Tableau r)
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (OptDir -> Tableau r -> Bool
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' -> Bool -> (Bool, Tableau r) -> (Bool, Tableau r)
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (OptDir -> Tableau r -> Tableau r -> Bool
forall r. Real r => OptDir -> Tableau r -> Tableau r -> Bool
isImproving OptDir
optdir Tableau r
tbl Tableau r
tbl') ((Bool, Tableau r) -> (Bool, Tableau r))
-> (Bool, Tableau r) -> (Bool, Tableau r)
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
  | [(RowIndex, r)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, r)]
cs   = PivotResult r
forall r. PivotResult r
PivotFinished
  | [(RowIndex, r)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, r)]
rs   = PivotResult r
forall r. PivotResult r
PivotUnbounded
  | Bool
otherwise = Tableau r -> PivotResult r
forall r. Tableau r -> PivotResult r
PivotSuccess (RowIndex -> RowIndex -> Tableau r -> Tableau r
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
0r -> r -> Bool
forall a. Ord a => a -> a -> Bool
<)
            OptDir
OptMax -> (r
0r -> r -> Bool
forall a. Ord a => a -> a -> Bool
>)
    cs :: [(RowIndex, r)]
cs = [(RowIndex
j,r
cj) | (RowIndex
j,r
cj) <- IntMap r -> [(RowIndex, r)]
forall a. IntMap a -> [(RowIndex, a)]
IM.toList ((IntMap r, r) -> IntMap r
forall a b. (a, b) -> a
fst (RowIndex -> Tableau r -> (IntMap r, r)
forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
objRowIndex Tableau r
tbl)), r -> Bool
cmp r
cj]
    -- smallest subscript rule
    s :: RowIndex
s = (RowIndex, r) -> RowIndex
forall a b. (a, b) -> a
fst ((RowIndex, r) -> RowIndex) -> (RowIndex, r) -> RowIndex
forall a b. (a -> b) -> a -> b
$ [(RowIndex, r)] -> (RowIndex, r)
forall a. (?callStack::CallStack) => [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 r -> r -> r
forall a. Fractional a => a -> a -> a
/ r
y_is)
         | (RowIndex
i, (IntMap r
row_i, r
y_i0)) <- Tableau r -> [(RowIndex, (IntMap r, r))]
forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
i RowIndex -> RowIndex -> Bool
forall a. Eq a => a -> a -> Bool
/= RowIndex
objRowIndex
         , let y_is :: r
y_is = r -> RowIndex -> IntMap r -> r
forall a. a -> RowIndex -> IntMap a -> a
IM.findWithDefault r
0 RowIndex
s IntMap r
row_i, r
y_is r -> r -> Bool
forall a. Ord a => a -> a -> Bool
> r
0
         ]
    r :: RowIndex
r = (RowIndex, r) -> RowIndex
forall a b. (a, b) -> a
fst ((RowIndex, r) -> RowIndex) -> (RowIndex, r) -> RowIndex
forall a b. (a -> b) -> a -> b
$ ((RowIndex, r) -> (RowIndex, r) -> Ordering)
-> [(RowIndex, r)] -> (RowIndex, r)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy (((RowIndex, r) -> r) -> (RowIndex, r) -> (RowIndex, r) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (RowIndex, r) -> r
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 = Tableau r -> (Bool, Tableau r)
forall {r}.
(Real r, Fractional r) =>
Tableau r -> (Bool, Tableau r)
go
  where
    go :: Tableau r -> (Bool, Tableau r)
go Tableau r
tbl = Bool -> (Bool, Tableau r) -> (Bool, Tableau r)
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (OptDir -> Tableau r -> Bool
forall r. Real r => OptDir -> Tableau r -> Bool
isOptimal OptDir
optdir Tableau r
tbl) ((Bool, Tableau r) -> (Bool, Tableau r))
-> (Bool, Tableau r) -> (Bool, Tableau r)
forall a b. (a -> b) -> a -> b
$
      case OptDir -> Tableau r -> PivotResult r
forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> PivotResult r
dualPivot OptDir
optdir Tableau r
tbl of
        PivotResult r
PivotFinished  -> Bool -> (Bool, Tableau r) -> (Bool, Tableau r)
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Tableau r -> Bool
forall r. Real r => Tableau r -> Bool
isFeasible Tableau r
tbl) ((Bool, Tableau r) -> (Bool, Tableau r))
-> (Bool, Tableau r) -> (Bool, Tableau r)
forall a b. (a -> b) -> a -> b
$ (Bool
True, Tableau r
tbl)
        PivotResult r
PivotUnbounded -> (Bool
False, Tableau r
tbl)
        PivotSuccess Tableau r
tbl' -> Bool -> (Bool, Tableau r) -> (Bool, Tableau r)
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (OptDir -> Tableau r -> Tableau r -> Bool
forall r. Real r => OptDir -> Tableau r -> Tableau r -> Bool
isImproving OptDir
optdir Tableau r
tbl' Tableau r
tbl) ((Bool, Tableau r) -> (Bool, Tableau r))
-> (Bool, Tableau r) -> (Bool, Tableau r)
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
  | [(RowIndex, VarMap r)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, VarMap r)]
rs   = PivotResult r
forall r. PivotResult r
PivotFinished
  | [(RowIndex, r)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, r)]
cs   = PivotResult r
forall r. PivotResult r
PivotUnbounded
  | Bool
otherwise = Tableau r -> PivotResult r
forall r. Tableau r -> PivotResult r
PivotSuccess (RowIndex -> RowIndex -> Tableau r -> Tableau r
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)) <- Tableau r -> [(RowIndex, Row r)]
forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
i RowIndex -> RowIndex -> Bool
forall a. Eq a => a -> a -> Bool
/= RowIndex
objRowIndex, r
0 r -> r -> Bool
forall a. Ord a => a -> a -> Bool
> r
y_i0]
    (RowIndex
r, VarMap r
row_r) = [(RowIndex, VarMap r)] -> (RowIndex, VarMap r)
forall a. (?callStack::CallStack) => [a] -> a
head [(RowIndex, VarMap r)]
rs
    cs :: [(RowIndex, r)]
cs = [ (RowIndex
j, if OptDir
optdirOptDir -> OptDir -> Bool
forall a. Eq a => a -> a -> Bool
==OptDir
OptMin then r
y_0j r -> r -> r
forall a. Fractional a => a -> a -> a
/ r
y_rj else - r
y_0j r -> r -> r
forall a. Fractional a => a -> a -> a
/ r
y_rj)
         | (RowIndex
j, r
y_rj) <- VarMap r -> [(RowIndex, r)]
forall a. IntMap a -> [(RowIndex, a)]
IM.toList VarMap r
row_r
         , r
y_rj r -> r -> Bool
forall a. Ord a => a -> a -> Bool
< r
0
         , let y_0j :: r
y_0j = r -> RowIndex -> VarMap r -> r
forall a. a -> RowIndex -> IntMap a -> a
IM.findWithDefault r
0 RowIndex
j VarMap r
obj
         ]
    s :: RowIndex
s = (RowIndex, r) -> RowIndex
forall a b. (a, b) -> a
fst ((RowIndex, r) -> RowIndex) -> (RowIndex, r) -> RowIndex
forall a b. (a -> b) -> a -> b
$ ((RowIndex, r) -> (RowIndex, r) -> Ordering)
-> [(RowIndex, r)] -> (RowIndex, r)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy (((RowIndex, r) -> r) -> (RowIndex, r) -> (RowIndex, r) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (RowIndex, r) -> r
forall a b. (a, b) -> b
snd) [(RowIndex, r)]
cs
    (VarMap r
obj,r
_) = RowIndex -> Tableau r -> Row 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
  | Tableau r -> r
forall r. Tableau r -> r
currentObjValue Tableau r
tbl1' r -> r -> Bool
forall a. Eq a => a -> a -> Bool
/= r
0 = (Bool
False, Tableau r
tbl1')
  | Bool
otherwise = (Bool
True, Tableau r -> Tableau r -> Tableau r
forall r. (Num r, Eq r) => Tableau r -> Tableau r -> Tableau r
copyObjRow Tableau r
tbl (Tableau r -> Tableau r) -> Tableau r -> Tableau r
forall a b. (a -> b) -> a -> b
$ IntSet -> Tableau r -> Tableau r
forall r.
(Real r, Fractional r) =>
IntSet -> Tableau r -> Tableau r
removeArtificialVariables IntSet
avs (Tableau r -> Tableau r) -> Tableau r -> Tableau r
forall a b. (a -> b) -> a -> b
$ Tableau r
tbl1')
  where
    optdir :: OptDir
optdir = OptDir
OptMax
    tbl1 :: Tableau r
tbl1 = Tableau r -> Expr r -> Tableau r
forall r. (Num r, Eq r) => Tableau r -> Expr r -> Tableau r
setObjFun Tableau r
tbl (Expr r -> Tableau r) -> Expr r -> Tableau r
forall a b. (a -> b) -> a -> b
$ Expr r -> Expr r
forall v. AdditiveGroup v => v -> v
negateV (Expr r -> Expr r) -> Expr r -> Expr r
forall a b. (a -> b) -> a -> b
$ [Expr r] -> Expr r
forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV [RowIndex -> Expr r
forall r. Num r => RowIndex -> Expr r
LA.var RowIndex
v | RowIndex
v <- IntSet -> [RowIndex]
IS.toList IntSet
avs]
    tbl1' :: Tableau r
tbl1' = Tableau r -> Tableau r
forall {r}. (Real r, Fractional r) => Tableau r -> Tableau r
go Tableau r
tbl1
    go :: Tableau r -> Tableau r
go Tableau r
tbl2
      | Tableau r -> r
forall r. Tableau r -> r
currentObjValue Tableau r
tbl2 r -> r -> Bool
forall a. Eq a => a -> a -> Bool
== r
0 = Tableau r
tbl2
      | Bool
otherwise =
        case OptDir -> Tableau r -> PivotResult r
forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> PivotResult r
primalPivot OptDir
optdir Tableau r
tbl2 of
          PivotSuccess Tableau r
tbl2' -> Bool -> Tableau r -> Tableau r
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (OptDir -> Tableau r -> Tableau r -> Bool
forall r. Real r => OptDir -> Tableau r -> Tableau r -> Bool
isImproving OptDir
optdir Tableau r
tbl2 Tableau r
tbl2') (Tableau r -> Tableau r) -> Tableau r -> Tableau r
forall a b. (a -> b) -> a -> b
$ Tableau r -> Tableau r
go Tableau r
tbl2'
          PivotResult r
PivotFinished -> Bool -> Tableau r -> Tableau r
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (OptDir -> Tableau r -> Bool
forall r. Real r => OptDir -> Tableau r -> Bool
isOptimal OptDir
optdir Tableau r
tbl2) Tableau r
tbl2
          PivotResult r
PivotUnbounded -> String -> Tableau r
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 = (Tableau r -> RowIndex -> Tableau r)
-> Tableau r -> [RowIndex] -> Tableau r
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Tableau r -> RowIndex -> Tableau r
forall {b}.
(Eq b, Fractional b) =>
IntMap (IntMap b, b) -> RowIndex -> IntMap (IntMap b, b)
f (RowIndex -> Tableau r -> Tableau r
forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
objRowIndex Tableau r
tbl0) (IntSet -> [RowIndex]
IS.toList IntSet
avs)
    tbl2 :: Tableau r
tbl2 = ((IntMap r, r) -> (IntMap r, r)) -> Tableau r -> Tableau r
forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map (\(IntMap r
row,r
val) ->  ((RowIndex -> r -> Bool) -> IntMap r -> IntMap r
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 RowIndex -> IntMap (IntMap b, b) -> Maybe (IntMap b, b)
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) <- IntMap b -> [(RowIndex, b)]
forall a. IntMap a -> [(RowIndex, a)]
IM.toList ((IntMap b, b) -> IntMap b
forall a b. (a, b) -> a
fst (IntMap b, b)
row), b
c b -> b -> Bool
forall a. Eq a => a -> a -> Bool
/= b
0, RowIndex
j RowIndex -> IntSet -> Bool
`IS.notMember` IntSet
avs] of
            [] -> RowIndex -> IntMap (IntMap b, b) -> IntMap (IntMap b, b)
forall a. RowIndex -> IntMap a -> IntMap a
IM.delete RowIndex
i IntMap (IntMap b, b)
tbl
            RowIndex
j:[RowIndex]
_ -> RowIndex
-> RowIndex -> IntMap (IntMap b, b) -> IntMap (IntMap b, b)
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 = Tableau r -> (Bool, Tableau r)
forall {r}.
(Real r, Fractional r) =>
Tableau r -> (Bool, Tableau r)
go
  where
    go :: Tableau r -> (Bool, Tableau r)
go Tableau r
tbl =
      case OptDir -> Tableau r -> Either PDResult (Tableau r)
forall r.
(Real r, Fractional r) =>
OptDir -> Tableau r -> Either PDResult (Tableau r)
pdPivot OptDir
optdir Tableau r
tbl of
        Left PDResult
PDOptimal   -> Bool -> (Bool, Tableau r) -> (Bool, Tableau r)
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Tableau r -> Bool
forall r. Real r => Tableau r -> Bool
isFeasible Tableau r
tbl) ((Bool, Tableau r) -> (Bool, Tableau r))
-> (Bool, Tableau r) -> (Bool, Tableau r)
forall a b. (a -> b) -> a -> b
$ Bool -> (Bool, Tableau r) -> (Bool, Tableau r)
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (OptDir -> Tableau r -> Bool
forall r. Real r => OptDir -> Tableau r -> Bool
isOptimal OptDir
optdir Tableau r
tbl) ((Bool, Tableau r) -> (Bool, Tableau r))
-> (Bool, Tableau r) -> (Bool, Tableau r)
forall a b. (a -> b) -> a -> b
$ (Bool
True, Tableau r
tbl)
        Left PDResult
PDUnsat     -> Bool -> (Bool, Tableau r) -> (Bool, Tableau r)
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Bool -> Bool
not (Tableau r -> Bool
forall r. Real r => Tableau r -> Bool
isFeasible Tableau r
tbl)) ((Bool, Tableau r) -> (Bool, Tableau r))
-> (Bool, Tableau r) -> (Bool, Tableau r)
forall a b. (a -> b) -> a -> b
$ (Bool
False, Tableau r
tbl)
        Left PDResult
PDUnbounded -> Bool -> (Bool, Tableau r) -> (Bool, Tableau r)
forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Bool -> Bool
not (OptDir -> Tableau r -> Bool
forall r. Real r => OptDir -> Tableau r -> Bool
isOptimal OptDir
optdir Tableau r
tbl)) ((Bool, Tableau r) -> (Bool, Tableau r))
-> (Bool, Tableau r) -> (Bool, Tableau r)
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
  | [(Either RowIndex Any, r)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(Either RowIndex Any, r)]
forall {b}. [(Either RowIndex b, r)]
ps Bool -> Bool -> Bool
&& [(Either Any RowIndex, r)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(Either Any RowIndex, r)]
forall {a}. [(Either a RowIndex, r)]
qs = PDResult -> Either PDResult (Tableau r)
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 r -> r -> r
forall a. Num a => a -> a -> a
- r
t) r -> r -> r
forall a. Fractional a => a -> a -> a
/ r
y_ip)
                   | (RowIndex
i, (VarMap r
row_i, r
bi)) <- Tableau r -> [(RowIndex, Row r)]
forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
i RowIndex -> RowIndex -> Bool
forall a. Eq a => a -> a -> Bool
/= RowIndex
objRowIndex
                   , let y_ip :: r
y_ip = r -> RowIndex -> VarMap r -> r
forall a. a -> RowIndex -> IntMap a -> a
IM.findWithDefault r
0 RowIndex
p VarMap r
row_i, r
y_ip r -> r -> Bool
forall a. Ord a => a -> a -> Bool
> r
0
                   ]
              q :: RowIndex
q = (RowIndex, r) -> RowIndex
forall a b. (a, b) -> a
fst ((RowIndex, r) -> RowIndex) -> (RowIndex, r) -> RowIndex
forall a b. (a -> b) -> a -> b
$ ((RowIndex, r) -> (RowIndex, r) -> Ordering)
-> [(RowIndex, r)] -> (RowIndex, r)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy (((RowIndex, r) -> r) -> (RowIndex, r) -> (RowIndex, r) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (RowIndex, r) -> r
forall a b. (a, b) -> b
snd) [(RowIndex, r)]
rs
          in if [(RowIndex, r)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, r)]
rs
             then PDResult -> Either PDResult (Tableau r)
forall a b. a -> Either a b
Left PDResult
PDUnsat
             else Tableau r -> Either PDResult (Tableau r)
forall a b. b -> Either a b
Right (RowIndex -> RowIndex -> Tableau r -> Tableau r
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 Tableau r -> RowIndex -> Row r
forall a. IntMap a -> RowIndex -> a
IM.! RowIndex
q)
              cs :: [(RowIndex, r)]
cs = [ (RowIndex
j, (r
cj'r -> r -> r
forall a. Num a => a -> a -> a
-r
t) r -> r -> r
forall a. Fractional a => a -> a -> a
/ (-r
y_qj))
                   | (RowIndex
j, r
y_qj) <- VarMap r -> [(RowIndex, r)]
forall a. IntMap a -> [(RowIndex, a)]
IM.toList VarMap r
row_q
                   , r
y_qj r -> r -> Bool
forall a. Ord a => a -> a -> Bool
< r
0
                   , let cj :: r
cj  = r -> RowIndex -> VarMap r -> r
forall a. a -> RowIndex -> IntMap a -> a
IM.findWithDefault r
0 RowIndex
j VarMap r
obj
                   , let cj' :: r
cj' = if OptDir
optdirOptDir -> OptDir -> Bool
forall a. Eq a => a -> a -> Bool
==OptDir
OptMax then r
cj else -r
cj
                   ]
              p :: RowIndex
p = (RowIndex, r) -> RowIndex
forall a b. (a, b) -> a
fst ((RowIndex, r) -> RowIndex) -> (RowIndex, r) -> RowIndex
forall a b. (a -> b) -> a -> b
$ ((RowIndex, r) -> (RowIndex, r) -> Ordering)
-> [(RowIndex, r)] -> (RowIndex, r)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy (((RowIndex, r) -> r) -> (RowIndex, r) -> (RowIndex, r) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (RowIndex, r) -> r
forall a b. (a, b) -> b
snd) [(RowIndex, r)]
cs
              (VarMap r
obj,r
_) = RowIndex -> Tableau r -> Row r
forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
objRowIndex Tableau r
tbl
          in if [(RowIndex, r)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(RowIndex, r)]
cs
             then PDResult -> Either PDResult (Tableau r)
forall a b. a -> Either a b
Left PDResult
PDUnbounded -- dual infeasible
             else Tableau r -> Either PDResult (Tableau r)
forall a b. b -> Either a b
Right (RowIndex -> RowIndex -> Tableau r -> Tableau r
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 = [ (RowIndex -> Either a RowIndex
forall a b. b -> Either a b
Right RowIndex
i, r
bi) | (RowIndex
i, (VarMap r
_row_i, r
bi)) <- Tableau r -> [(RowIndex, Row r)]
forall a. IntMap a -> [(RowIndex, a)]
IM.toList Tableau r
tbl, RowIndex
i RowIndex -> RowIndex -> Bool
forall a. Eq a => a -> a -> Bool
/= RowIndex
objRowIndex, r
0 r -> r -> Bool
forall a. Ord a => a -> a -> Bool
> r
bi ]
    ps :: [(Either RowIndex b, r)]
ps = [ (RowIndex -> Either RowIndex b
forall a b. a -> Either a b
Left RowIndex
j, r
cj')
         | (RowIndex
j,r
cj) <- VarMap r -> [(RowIndex, r)]
forall a. IntMap a -> [(RowIndex, a)]
IM.toList (Row r -> VarMap r
forall a b. (a, b) -> a
fst (RowIndex -> Tableau r -> Row r
forall r. RowIndex -> Tableau r -> Row r
lookupRow RowIndex
objRowIndex Tableau r
tbl))
         , let cj' :: r
cj' = if OptDir
optdirOptDir -> OptDir -> Bool
forall a. Eq a => a -> a -> Bool
==OptDir
OptMax then r
cj else -r
cj
         , r
0 r -> r -> Bool
forall a. Ord a => a -> a -> Bool
> r
cj' ]
    (Either RowIndex RowIndex
ret, r
t) = ((Either RowIndex RowIndex, r)
 -> (Either RowIndex RowIndex, r) -> Ordering)
-> [(Either RowIndex RowIndex, r)] -> (Either RowIndex RowIndex, r)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
minimumBy (((Either RowIndex RowIndex, r) -> r)
-> (Either RowIndex RowIndex, r)
-> (Either RowIndex RowIndex, r)
-> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (Either RowIndex RowIndex, r) -> r
forall a b. (a, b) -> b
snd) ([(Either RowIndex RowIndex, r)]
forall {a}. [(Either a RowIndex, r)]
qs [(Either RowIndex RowIndex, r)]
-> [(Either RowIndex RowIndex, r)]
-> [(Either RowIndex RowIndex, r)]
forall a. [a] -> [a] -> [a]
++ [(Either RowIndex RowIndex, r)]
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 ([String] -> String)
-> ([[String]] -> [String]) -> [[String]] -> String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([String] -> String) -> [[String]] -> [String]
forall a b. (a -> b) -> [a] -> [b]
map ([String] -> String
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([String] -> String)
-> ([String] -> [String]) -> [String] -> String
forall b c a. (b -> c) -> (a -> b) -> a -> c
. String -> [String] -> [String]
forall a. a -> [a] -> [a]
intersperse String
",") ([[String]] -> String) -> [[String]] -> String
forall a b. (a -> b) -> a -> b
$ [String]
header [String] -> [[String]] -> [[String]]
forall a. a -> [a] -> [a]
: [[String]]
body
  where
    header :: [String]
    header :: [String]
header = String
"" String -> [String] -> [String]
forall a. a -> [a] -> [a]
: (RowIndex -> String) -> [RowIndex] -> [String]
forall a b. (a -> b) -> [a] -> [b]
map RowIndex -> String
colName [RowIndex]
cols [String] -> [String] -> [String]
forall a. [a] -> [a] -> [a]
++ [String
""]

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

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

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

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

    colName :: ColIndex -> String
    colName :: RowIndex -> String
colName RowIndex
j = String
"x" String -> ShowS
forall a. [a] -> [a] -> [a]
++ RowIndex -> String
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 String -> [String] -> [String]
forall a. a -> [a] -> [a]
: [r -> String
showCell (r -> RowIndex -> IntMap r -> r
forall a. a -> RowIndex -> IntMap a -> a
IM.findWithDefault r
0 RowIndex
j IntMap r
row') | RowIndex
j <- [RowIndex]
cols] [String] -> [String] -> [String]
forall a. [a] -> [a] -> [a]
++ [r -> String
showCell r
row_val]
      where row' :: IntMap r
row' = RowIndex -> r -> IntMap r -> IntMap r
forall a. RowIndex -> a -> IntMap a -> IntMap a
IM.insert RowIndex
i r
1 IntMap r
row

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