{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE TupleSections #-}

{-|
Module      : Linear.Simplex.Simplex
Description : Implements the twoPhaseSimplex method
Copyright   : (c) Junaid Rasheed, 2020-2022
License     : BSD-3
Maintainer  : jrasheed178@gmail.com
Stability   : experimental

Module implementing the two-phase simplex method.
'findFeasibleSolution' performs phase one of the two-phase simplex method.
'optimizeFeasibleSystem' performs phase two of the two-phase simplex method.
'twoPhaseSimplex' performs both phases of the two-phase simplex method. 
-}
module Linear.Simplex.Simplex (findFeasibleSolution, optimizeFeasibleSystem, twoPhaseSimplex) where
import Linear.Simplex.Types
import Linear.Simplex.Util
import Prelude hiding (EQ);
import Data.List
import Data.Bifunctor
import Data.Maybe (fromMaybe, mapMaybe)
import Data.Ratio (numerator, denominator, (%))
-- import Debug.Trace (trace)

trace :: p -> p -> p
trace p
s p
a = p
a

-- |Find a feasible solution for the given system of 'PolyConstraint's by performing the first phase of the two-phase simplex method
-- All 'Integer' variables in the 'PolyConstraint' must be positive.
-- If the system is infeasible, return 'Nothing'
-- Otherwise, return the feasible system in 'DictionaryForm' as well as a list of slack variables, a list artificial variables, and the objective variable.
findFeasibleSolution :: [PolyConstraint] -> Maybe (DictionaryForm, [Integer], [Integer], Integer)
findFeasibleSolution :: [PolyConstraint]
-> Maybe (DictionaryForm, [Integer], [Integer], Integer)
findFeasibleSolution [PolyConstraint]
unsimplifiedSystem = 
  if [Integer] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Integer]
artificialVars -- No artificial vars, we have a feasible system
    then (DictionaryForm, [Integer], [Integer], Integer)
-> Maybe (DictionaryForm, [Integer], [Integer], Integer)
forall a. a -> Maybe a
Just (DictionaryForm
systemWithBasicVarsAsDictionary, [Integer]
slackVars, [Integer]
artificialVars, Integer
objectiveVar)
    else 
      case DictionaryForm -> Maybe DictionaryForm
simplexPivot (ObjectiveFunction -> Integer -> (Integer, VarConstMap)
createObjectiveDict ObjectiveFunction
artificialObjective Integer
objectiveVar (Integer, VarConstMap) -> DictionaryForm -> DictionaryForm
forall a. a -> [a] -> [a]
: DictionaryForm
systemWithBasicVarsAsDictionary) of
        Just DictionaryForm
phase1Dict ->
          let
            eliminateArtificialVarsFromPhase1Tableau :: DictionaryForm
eliminateArtificialVarsFromPhase1Tableau = ((Integer, VarConstMap) -> (Integer, VarConstMap))
-> DictionaryForm -> DictionaryForm
forall a b. (a -> b) -> [a] -> [b]
map ((VarConstMap -> VarConstMap)
-> (Integer, VarConstMap) -> (Integer, VarConstMap)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second (((Integer, Rational) -> Bool) -> VarConstMap -> VarConstMap
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Integer
v, Rational
_) -> Integer
v Integer -> [Integer] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Integer]
artificialVars))) DictionaryForm
phase1Dict
          in
            case Integer -> DictionaryForm -> Maybe VarConstMap
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Integer
objectiveVar DictionaryForm
eliminateArtificialVarsFromPhase1Tableau of
              Maybe VarConstMap
Nothing -> String
-> Maybe (DictionaryForm, [Integer], [Integer], Integer)
-> Maybe (DictionaryForm, [Integer], [Integer], Integer)
forall {p} {p}. p -> p -> p
trace String
"objective row not found in phase 1 tableau" Maybe (DictionaryForm, [Integer], [Integer], Integer)
forall a. Maybe a
Nothing -- Should this be an error?
              Just VarConstMap
row ->
                if Rational -> Maybe Rational -> Rational
forall a. a -> Maybe a -> a
fromMaybe Rational
0 (Integer -> VarConstMap -> Maybe Rational
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup (-Integer
1) VarConstMap
row) Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0
                  then (DictionaryForm, [Integer], [Integer], Integer)
-> Maybe (DictionaryForm, [Integer], [Integer], Integer)
forall a. a -> Maybe a
Just (DictionaryForm
eliminateArtificialVarsFromPhase1Tableau, [Integer]
slackVars, [Integer]
artificialVars, Integer
objectiveVar)
                  else String
-> Maybe (DictionaryForm, [Integer], [Integer], Integer)
-> Maybe (DictionaryForm, [Integer], [Integer], Integer)
forall {p} {p}. p -> p -> p
trace String
"rhs not zero after phase 1, thus original tableau is infeasible" Maybe (DictionaryForm, [Integer], [Integer], Integer)
forall a. Maybe a
Nothing 
        Maybe DictionaryForm
Nothing -> Maybe (DictionaryForm, [Integer], [Integer], Integer)
forall a. Maybe a
Nothing
  where
    system :: [PolyConstraint]
system = [PolyConstraint] -> [PolyConstraint]
simplifySystem [PolyConstraint]
unsimplifiedSystem

    maxVar :: Integer
maxVar =
      [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum ([Integer] -> Integer) -> [Integer] -> Integer
forall a b. (a -> b) -> a -> b
$ (PolyConstraint -> Integer) -> [PolyConstraint] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map 
      (\case
          LEQ VarConstMap
vcm Rational
_ -> [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (((Integer, Rational) -> Integer) -> VarConstMap -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer, Rational) -> Integer
forall a b. (a, b) -> a
fst VarConstMap
vcm)
          GEQ VarConstMap
vcm Rational
_ -> [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (((Integer, Rational) -> Integer) -> VarConstMap -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer, Rational) -> Integer
forall a b. (a, b) -> a
fst VarConstMap
vcm)
          EQ VarConstMap
vcm Rational
_  -> [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (((Integer, Rational) -> Integer) -> VarConstMap -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer, Rational) -> Integer
forall a b. (a, b) -> a
fst VarConstMap
vcm)
      ) 
      [PolyConstraint]
system

    ([(Maybe Integer, PolyConstraint)]
systemWithSlackVars, [Integer]
slackVars) = [PolyConstraint]
-> Integer
-> [Integer]
-> ([(Maybe Integer, PolyConstraint)], [Integer])
systemInStandardForm [PolyConstraint]
system Integer
maxVar []

    maxVarWithSlackVars :: Integer
maxVarWithSlackVars = if [Integer] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Integer]
slackVars then Integer
maxVar else [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Integer]
slackVars

    (Tableau
systemWithBasicVars, [Integer]
artificialVars) = [(Maybe Integer, PolyConstraint)]
-> Integer -> (Tableau, [Integer])
systemWithArtificialVars [(Maybe Integer, PolyConstraint)]
systemWithSlackVars Integer
maxVarWithSlackVars 

    finalMaxVar :: Integer
finalMaxVar        = if [Integer] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Integer]
artificialVars then Integer
maxVarWithSlackVars else [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Integer]
artificialVars

    systemWithBasicVarsAsDictionary :: DictionaryForm
systemWithBasicVarsAsDictionary = Tableau -> DictionaryForm
tableauInDictionaryForm Tableau
systemWithBasicVars
    
    artificialObjective :: ObjectiveFunction
artificialObjective = DictionaryForm -> [Integer] -> ObjectiveFunction
createArtificialObjective DictionaryForm
systemWithBasicVarsAsDictionary [Integer]
artificialVars
    
    objectiveVar :: Integer
objectiveVar  = Integer
finalMaxVar Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1

    -- |Convert a system of 'PolyConstraint's to standard form; a system of only equations ('EQ').
    -- Add slack vars where necessary.
    -- This may give you an infeasible system if slack vars are negative when original variables are zero.
    -- If a constraint is already EQ, set the basic var to Nothing.
    -- Final system is a list of equalities for the given system. 
    -- To be feasible, all vars must be >= 0.
    systemInStandardForm :: [PolyConstraint] -> Integer -> [Integer] -> ([(Maybe Integer, PolyConstraint)], [Integer])
    systemInStandardForm :: [PolyConstraint]
-> Integer
-> [Integer]
-> ([(Maybe Integer, PolyConstraint)], [Integer])
systemInStandardForm []  Integer
_       [Integer]
sVars = ([], [Integer]
sVars)
    systemInStandardForm (EQ VarConstMap
v Rational
r : [PolyConstraint]
xs) Integer
maxVar [Integer]
sVars = ((Maybe Integer
forall a. Maybe a
Nothing, VarConstMap -> Rational -> PolyConstraint
EQ VarConstMap
v Rational
r) (Maybe Integer, PolyConstraint)
-> [(Maybe Integer, PolyConstraint)]
-> [(Maybe Integer, PolyConstraint)]
forall a. a -> [a] -> [a]
: [(Maybe Integer, PolyConstraint)]
newSystem, [Integer]
newSlackVars) 
      where
        ([(Maybe Integer, PolyConstraint)]
newSystem, [Integer]
newSlackVars) = [PolyConstraint]
-> Integer
-> [Integer]
-> ([(Maybe Integer, PolyConstraint)], [Integer])
systemInStandardForm [PolyConstraint]
xs Integer
maxVar [Integer]
sVars
    systemInStandardForm (LEQ VarConstMap
v Rational
r : [PolyConstraint]
xs) Integer
maxVar  [Integer]
sVars = ((Integer -> Maybe Integer
forall a. a -> Maybe a
Just Integer
newSlackVar, VarConstMap -> Rational -> PolyConstraint
EQ (VarConstMap
v VarConstMap -> VarConstMap -> VarConstMap
forall a. [a] -> [a] -> [a]
++ [(Integer
newSlackVar, Rational
1)]) Rational
r) (Maybe Integer, PolyConstraint)
-> [(Maybe Integer, PolyConstraint)]
-> [(Maybe Integer, PolyConstraint)]
forall a. a -> [a] -> [a]
: [(Maybe Integer, PolyConstraint)]
newSystem, [Integer]
newSlackVars)
      where
        newSlackVar :: Integer
newSlackVar = Integer
maxVar Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1
        ([(Maybe Integer, PolyConstraint)]
newSystem, [Integer]
newSlackVars) = [PolyConstraint]
-> Integer
-> [Integer]
-> ([(Maybe Integer, PolyConstraint)], [Integer])
systemInStandardForm [PolyConstraint]
xs Integer
newSlackVar (Integer
newSlackVar Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: [Integer]
sVars)
    systemInStandardForm (GEQ VarConstMap
v Rational
r : [PolyConstraint]
xs) Integer
maxVar  [Integer]
sVars = ((Integer -> Maybe Integer
forall a. a -> Maybe a
Just Integer
newSlackVar, VarConstMap -> Rational -> PolyConstraint
EQ (VarConstMap
v VarConstMap -> VarConstMap -> VarConstMap
forall a. [a] -> [a] -> [a]
++ [(Integer
newSlackVar, -Rational
1)]) Rational
r) (Maybe Integer, PolyConstraint)
-> [(Maybe Integer, PolyConstraint)]
-> [(Maybe Integer, PolyConstraint)]
forall a. a -> [a] -> [a]
: [(Maybe Integer, PolyConstraint)]
newSystem, [Integer]
newSlackVars)
      where
        newSlackVar :: Integer
newSlackVar = Integer
maxVar Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1
        ([(Maybe Integer, PolyConstraint)]
newSystem, [Integer]
newSlackVars) = [PolyConstraint]
-> Integer
-> [Integer]
-> ([(Maybe Integer, PolyConstraint)], [Integer])
systemInStandardForm [PolyConstraint]
xs Integer
newSlackVar (Integer
newSlackVar Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: [Integer]
sVars)

    -- |Add artificial vars to a system of 'PolyConstraint's.
    -- Artificial vars are added when:
    --  Basic var is Nothing (When the original constraint was already an EQ).
    --  Slack var is equal to a negative value (this is infeasible, all vars need to be >= 0).
    --  Final system will be a feasible artificial system.
    -- We keep track of artificial vars in the second item of the returned pair so they can be eliminated once phase 1 is complete.
    -- If an artificial var would normally be negative, we negate the row so we can keep artificial variables equal to 1
    systemWithArtificialVars :: [(Maybe Integer, PolyConstraint)] -> Integer -> (Tableau, [Integer])
    systemWithArtificialVars :: [(Maybe Integer, PolyConstraint)]
-> Integer -> (Tableau, [Integer])
systemWithArtificialVars [] Integer
_                                = ([],[])
    systemWithArtificialVars ((Maybe Integer
mVar, EQ VarConstMap
v Rational
r) : [(Maybe Integer, PolyConstraint)]
pcs) Integer
maxVar  =
      case Maybe Integer
mVar of
        Maybe Integer
Nothing ->
          if Rational
r Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
0 
            then 
              ((Integer
newArtificialVar, (VarConstMap
v VarConstMap -> VarConstMap -> VarConstMap
forall a. [a] -> [a] -> [a]
++ [(Integer
newArtificialVar, Rational
1)], Rational
r)) (Integer, (VarConstMap, Rational)) -> Tableau -> Tableau
forall a. a -> [a] -> [a]
: Tableau
newSystemWithNewMaxVar, Integer
newArtificialVar Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: [Integer]
artificialVarsWithNewMaxVar)
            else 
              ((Integer
newArtificialVar, (VarConstMap
v VarConstMap -> VarConstMap -> VarConstMap
forall a. [a] -> [a] -> [a]
++ [(Integer
newArtificialVar, -Rational
1)], Rational
r)) (Integer, (VarConstMap, Rational)) -> Tableau -> Tableau
forall a. a -> [a] -> [a]
: Tableau
newSystemWithNewMaxVar, Integer
newArtificialVar Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: [Integer]
artificialVarsWithNewMaxVar)
        Just Integer
basicVar ->
          case Integer -> VarConstMap -> Maybe Rational
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Integer
basicVar VarConstMap
v of
            Just Rational
basicVarCoeff ->
              if Rational
r Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Rational
0
                then ((Integer
basicVar, (VarConstMap
v, Rational
r)) (Integer, (VarConstMap, Rational)) -> Tableau -> Tableau
forall a. a -> [a] -> [a]
: Tableau
newSystemWithoutNewMaxVar, [Integer]
artificialVarsWithoutNewMaxVar)
                else
                  if Rational
r Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
> Rational
0
                    then 
                      if Rational
basicVarCoeff Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
0 -- Should only be 1 in the standard call path
                        then ((Integer
basicVar, (VarConstMap
v, Rational
r)) (Integer, (VarConstMap, Rational)) -> Tableau -> Tableau
forall a. a -> [a] -> [a]
: Tableau
newSystemWithoutNewMaxVar, [Integer]
artificialVarsWithoutNewMaxVar)
                        else ((Integer
newArtificialVar, (VarConstMap
v VarConstMap -> VarConstMap -> VarConstMap
forall a. [a] -> [a] -> [a]
++ [(Integer
newArtificialVar, Rational
1)], Rational
r)) (Integer, (VarConstMap, Rational)) -> Tableau -> Tableau
forall a. a -> [a] -> [a]
: Tableau
newSystemWithNewMaxVar, Integer
newArtificialVar Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: [Integer]
artificialVarsWithNewMaxVar) -- Slack var is negative, r is positive (when original constraint was GEQ)
                    else -- r < 0
                      if Rational
basicVarCoeff Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
0 -- Should only be -1 in the standard call path
                        then ((Integer
basicVar, (VarConstMap
v, Rational
r)) (Integer, (VarConstMap, Rational)) -> Tableau -> Tableau
forall a. a -> [a] -> [a]
: Tableau
newSystemWithoutNewMaxVar, [Integer]
artificialVarsWithoutNewMaxVar)
                        else ((Integer
newArtificialVar, (VarConstMap
v VarConstMap -> VarConstMap -> VarConstMap
forall a. [a] -> [a] -> [a]
++ [(Integer
newArtificialVar, -Rational
1)], Rational
r)) (Integer, (VarConstMap, Rational)) -> Tableau -> Tableau
forall a. a -> [a] -> [a]
: Tableau
newSystemWithNewMaxVar, Integer
newArtificialVar Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: [Integer]
artificialVarsWithNewMaxVar) -- Slack var is negative, r is negative (when original constraint was LEQ)
      where
        newArtificialVar :: Integer
newArtificialVar = Integer
maxVar Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1

        (Tableau
newSystemWithNewMaxVar, [Integer]
artificialVarsWithNewMaxVar) = [(Maybe Integer, PolyConstraint)]
-> Integer -> (Tableau, [Integer])
systemWithArtificialVars [(Maybe Integer, PolyConstraint)]
pcs Integer
newArtificialVar

        (Tableau
newSystemWithoutNewMaxVar, [Integer]
artificialVarsWithoutNewMaxVar) = [(Maybe Integer, PolyConstraint)]
-> Integer -> (Tableau, [Integer])
systemWithArtificialVars [(Maybe Integer, PolyConstraint)]
pcs Integer
maxVar

    -- |Create an artificial objective using the given 'Integer' list of artificialVars and the given 'DictionaryForm'.
    -- The artificial 'ObjectiveFunction' is the negated sum of all artificial vars.
    createArtificialObjective :: DictionaryForm -> [Integer] -> ObjectiveFunction
    createArtificialObjective :: DictionaryForm -> [Integer] -> ObjectiveFunction
createArtificialObjective DictionaryForm
rows [Integer]
artificialVars = VarConstMap -> ObjectiveFunction
Max VarConstMap
negatedSumWithoutArtificialVars
      where
        rowsToAdd :: DictionaryForm
rowsToAdd = ((Integer, VarConstMap) -> Bool)
-> DictionaryForm -> DictionaryForm
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Integer
i, VarConstMap
_) -> Integer
i Integer -> [Integer] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Integer]
artificialVars) DictionaryForm
rows
        negatedRows :: [VarConstMap]
negatedRows = ((Integer, VarConstMap) -> VarConstMap)
-> DictionaryForm -> [VarConstMap]
forall a b. (a -> b) -> [a] -> [b]
map (\(Integer
_, VarConstMap
vcm) -> ((Integer, Rational) -> (Integer, Rational))
-> VarConstMap -> VarConstMap
forall a b. (a -> b) -> [a] -> [b]
map ((Rational -> Rational)
-> (Integer, Rational) -> (Integer, Rational)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second Rational -> Rational
forall a. Num a => a -> a
negate) VarConstMap
vcm) DictionaryForm
rowsToAdd
        negatedSum :: VarConstMap
negatedSum = VarConstMap -> VarConstMap
foldSumVarConstMap ((VarConstMap -> VarConstMap
forall a. Ord a => [a] -> [a]
sort (VarConstMap -> VarConstMap)
-> ([VarConstMap] -> VarConstMap) -> [VarConstMap] -> VarConstMap
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [VarConstMap] -> VarConstMap
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat) [VarConstMap]
negatedRows) 
        negatedSumWithoutArtificialVars :: VarConstMap
negatedSumWithoutArtificialVars = ((Integer, Rational) -> Bool) -> VarConstMap -> VarConstMap
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Integer
v, Rational
_) -> Integer
v Integer -> [Integer] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Integer]
artificialVars) VarConstMap
negatedSum


-- |Optimize a feasible system by performing the second phase of the two-phase simplex method.
-- We first pass an 'ObjectiveFunction'.
-- Then, the feasible system in 'DictionaryForm' as well as a list of slack variables, a list artificial variables, and the objective variable.
-- Returns a pair with the first item being the 'Integer' variable equal to the 'ObjectiveFunction'
-- and the second item being a map of the values of all 'Integer' variables appearing in the system, including the 'ObjectiveFunction'.
optimizeFeasibleSystem :: ObjectiveFunction -> DictionaryForm -> [Integer] -> [Integer] -> Integer -> Maybe (Integer, [(Integer, Rational)])
optimizeFeasibleSystem :: ObjectiveFunction
-> DictionaryForm
-> [Integer]
-> [Integer]
-> Integer
-> Maybe (Integer, VarConstMap)
optimizeFeasibleSystem ObjectiveFunction
unsimplifiedObjFunction DictionaryForm
phase1Dict [Integer]
slackVars [Integer]
artificialVars Integer
objectiveVar =
  if [Integer] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Integer]
artificialVars
    then Tableau -> (Integer, VarConstMap)
displayResults (Tableau -> (Integer, VarConstMap))
-> (DictionaryForm -> Tableau)
-> DictionaryForm
-> (Integer, VarConstMap)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. DictionaryForm -> Tableau
dictionaryFormToTableau (DictionaryForm -> (Integer, VarConstMap))
-> Maybe DictionaryForm -> Maybe (Integer, VarConstMap)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> DictionaryForm -> Maybe DictionaryForm
simplexPivot (ObjectiveFunction -> Integer -> (Integer, VarConstMap)
createObjectiveDict ObjectiveFunction
objFunction Integer
objectiveVar (Integer, VarConstMap) -> DictionaryForm -> DictionaryForm
forall a. a -> [a] -> [a]
: DictionaryForm
phase1Dict)
    else Tableau -> (Integer, VarConstMap)
displayResults (Tableau -> (Integer, VarConstMap))
-> (DictionaryForm -> Tableau)
-> DictionaryForm
-> (Integer, VarConstMap)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. DictionaryForm -> Tableau
dictionaryFormToTableau (DictionaryForm -> (Integer, VarConstMap))
-> Maybe DictionaryForm -> Maybe (Integer, VarConstMap)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> DictionaryForm -> Maybe DictionaryForm
simplexPivot (ObjectiveFunction -> Integer -> (Integer, VarConstMap)
createObjectiveDict ObjectiveFunction
phase2ObjFunction Integer
objectiveVar (Integer, VarConstMap) -> DictionaryForm -> DictionaryForm
forall a. a -> [a] -> [a]
: DictionaryForm -> DictionaryForm
forall a. [a] -> [a]
tail DictionaryForm
phase1Dict)
  where
    objFunction :: ObjectiveFunction
objFunction = ObjectiveFunction -> ObjectiveFunction
simplifyObjectiveFunction ObjectiveFunction
unsimplifiedObjFunction

    displayResults :: Tableau -> (Integer, [(Integer, Rational)])
    displayResults :: Tableau -> (Integer, VarConstMap)
displayResults Tableau
tableau =
      (
        Integer
objectiveVar,
        case ObjectiveFunction
objFunction of
          Max VarConstMap
_ -> 
            ((Integer, (VarConstMap, Rational)) -> (Integer, Rational))
-> Tableau -> VarConstMap
forall a b. (a -> b) -> [a] -> [b]
map 
            (((VarConstMap, Rational) -> Rational)
-> (Integer, (VarConstMap, Rational)) -> (Integer, Rational)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second (VarConstMap, Rational) -> Rational
forall a b. (a, b) -> b
snd) 
            (Tableau -> VarConstMap) -> Tableau -> VarConstMap
forall a b. (a -> b) -> a -> b
$ ((Integer, (VarConstMap, Rational)) -> Bool) -> Tableau -> Tableau
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Integer
basicVar,(VarConstMap, Rational)
_) -> Integer
basicVar Integer -> [Integer] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Integer]
slackVars [Integer] -> [Integer] -> [Integer]
forall a. [a] -> [a] -> [a]
++ [Integer]
artificialVars) Tableau
tableau
          Min VarConstMap
_ -> 
            ((Integer, (VarConstMap, Rational)) -> (Integer, Rational))
-> Tableau -> VarConstMap
forall a b. (a -> b) -> [a] -> [b]
map -- We maximized -objVar, so we negate the objVar to get the final value
            (\(Integer
basicVar, (VarConstMap, Rational)
row) -> if Integer
basicVar Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
objectiveVar then (Integer
basicVar, Rational -> Rational
forall a. Num a => a -> a
negate ((VarConstMap, Rational) -> Rational
forall a b. (a, b) -> b
snd (VarConstMap, Rational)
row)) else (Integer
basicVar, (VarConstMap, Rational) -> Rational
forall a b. (a, b) -> b
snd (VarConstMap, Rational)
row))
            (Tableau -> VarConstMap) -> Tableau -> VarConstMap
forall a b. (a -> b) -> a -> b
$ ((Integer, (VarConstMap, Rational)) -> Bool) -> Tableau -> Tableau
forall a. (a -> Bool) -> [a] -> [a]
filter (\(Integer
basicVar,(VarConstMap, Rational)
_) -> Integer
basicVar Integer -> [Integer] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [Integer]
slackVars [Integer] -> [Integer] -> [Integer]
forall a. [a] -> [a] -> [a]
++ [Integer]
artificialVars) Tableau
tableau
      )

    phase2Objective :: VarConstMap
phase2Objective = 
      (VarConstMap -> VarConstMap
foldSumVarConstMap (VarConstMap -> VarConstMap)
-> (VarConstMap -> VarConstMap) -> VarConstMap -> VarConstMap
forall b c a. (b -> c) -> (a -> b) -> a -> c
. VarConstMap -> VarConstMap
forall a. Ord a => [a] -> [a]
sort) (VarConstMap -> VarConstMap) -> VarConstMap -> VarConstMap
forall a b. (a -> b) -> a -> b
$
        ((Integer, Rational) -> VarConstMap) -> VarConstMap -> VarConstMap
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap
        (\(Integer
var, Rational
coeff) ->
          case Integer -> DictionaryForm -> Maybe VarConstMap
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Integer
var DictionaryForm
phase1Dict of
            Maybe VarConstMap
Nothing -> [(Integer
var, Rational
coeff)]
            Just VarConstMap
row -> ((Integer, Rational) -> (Integer, Rational))
-> VarConstMap -> VarConstMap
forall a b. (a -> b) -> [a] -> [b]
map ((Rational -> Rational)
-> (Integer, Rational) -> (Integer, Rational)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second (Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
*Rational
coeff)) VarConstMap
row
        )  
        (ObjectiveFunction -> VarConstMap
getObjective ObjectiveFunction
objFunction)

    phase2ObjFunction :: ObjectiveFunction
phase2ObjFunction = if ObjectiveFunction -> Bool
isMax ObjectiveFunction
objFunction then VarConstMap -> ObjectiveFunction
Max VarConstMap
phase2Objective else VarConstMap -> ObjectiveFunction
Min VarConstMap
phase2Objective

-- |Perform the two phase simplex method with a given 'ObjectiveFunction' a system of 'PolyConstraint's.
-- Assumes the 'ObjectiveFunction' and 'PolyConstraint' is not empty. 
-- Returns a pair with the first item being the 'Integer' variable equal to the 'ObjectiveFunction'
-- and the second item being a map of the values of all 'Integer' variables appearing in the system, including the 'ObjectiveFunction'.
twoPhaseSimplex :: ObjectiveFunction -> [PolyConstraint] -> Maybe (Integer, [(Integer, Rational)])
twoPhaseSimplex :: ObjectiveFunction
-> [PolyConstraint] -> Maybe (Integer, VarConstMap)
twoPhaseSimplex ObjectiveFunction
objFunction [PolyConstraint]
unsimplifiedSystem = 
  case [PolyConstraint]
-> Maybe (DictionaryForm, [Integer], [Integer], Integer)
findFeasibleSolution [PolyConstraint]
unsimplifiedSystem of
    Just r :: (DictionaryForm, [Integer], [Integer], Integer)
r@(DictionaryForm
phase1Dict, [Integer]
slackVars, [Integer]
artificialVars, Integer
objectiveVar) -> ObjectiveFunction
-> DictionaryForm
-> [Integer]
-> [Integer]
-> Integer
-> Maybe (Integer, VarConstMap)
optimizeFeasibleSystem ObjectiveFunction
objFunction DictionaryForm
phase1Dict [Integer]
slackVars [Integer]
artificialVars Integer
objectiveVar
    Maybe (DictionaryForm, [Integer], [Integer], Integer)
Nothing -> Maybe (Integer, VarConstMap)
forall a. Maybe a
Nothing

-- |Perform the simplex pivot algorithm on a system with basic vars, assume that the first row is the 'ObjectiveFunction'.
simplexPivot :: DictionaryForm -> Maybe DictionaryForm
simplexPivot :: DictionaryForm -> Maybe DictionaryForm
simplexPivot DictionaryForm
dictionary = 
  String -> Maybe DictionaryForm -> Maybe DictionaryForm
forall {p} {p}. p -> p -> p
trace (DictionaryForm -> String
forall a. Show a => a -> String
show DictionaryForm
dictionary) (Maybe DictionaryForm -> Maybe DictionaryForm)
-> Maybe DictionaryForm -> Maybe DictionaryForm
forall a b. (a -> b) -> a -> b
$
  case (Integer, VarConstMap) -> Maybe Integer
mostPositive (DictionaryForm -> (Integer, VarConstMap)
forall a. [a] -> a
head DictionaryForm
dictionary) of
    Maybe Integer
Nothing -> 
      String
-> (String
    -> (DictionaryForm -> Maybe DictionaryForm)
    -> DictionaryForm
    -> Maybe DictionaryForm)
-> String
-> (DictionaryForm -> Maybe DictionaryForm)
-> DictionaryForm
-> Maybe DictionaryForm
forall {p} {p}. p -> p -> p
trace String
"all neg \n"
      String
-> (DictionaryForm -> Maybe DictionaryForm)
-> DictionaryForm
-> Maybe DictionaryForm
forall {p} {p}. p -> p -> p
trace (DictionaryForm -> String
forall a. Show a => a -> String
show DictionaryForm
dictionary)
      DictionaryForm -> Maybe DictionaryForm
forall a. a -> Maybe a
Just DictionaryForm
dictionary
    Just Integer
pivotNonBasicVar -> 
      let
        mPivotBasicVar :: Maybe Integer
mPivotBasicVar = DictionaryForm
-> Integer -> Maybe Integer -> Maybe Rational -> Maybe Integer
ratioTest (DictionaryForm -> DictionaryForm
forall a. [a] -> [a]
tail DictionaryForm
dictionary) Integer
pivotNonBasicVar Maybe Integer
forall a. Maybe a
Nothing Maybe Rational
forall a. Maybe a
Nothing
      in
        case Maybe Integer
mPivotBasicVar of
          Maybe Integer
Nothing -> String -> Maybe DictionaryForm -> Maybe DictionaryForm
forall {p} {p}. p -> p -> p
trace (String
"Ratio test failed on non-basic var: " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show Integer
pivotNonBasicVar String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
"\n" String -> String -> String
forall a. [a] -> [a] -> [a]
++ DictionaryForm -> String
forall a. Show a => a -> String
show DictionaryForm
dictionary) Maybe DictionaryForm
forall a. Maybe a
Nothing
          Just Integer
pivotBasicVar -> 
            String
-> (String
    -> (DictionaryForm -> Maybe DictionaryForm)
    -> DictionaryForm
    -> Maybe DictionaryForm)
-> String
-> (DictionaryForm -> Maybe DictionaryForm)
-> DictionaryForm
-> Maybe DictionaryForm
forall {p} {p}. p -> p -> p
trace String
"one pos \n"
            String
-> (DictionaryForm -> Maybe DictionaryForm)
-> DictionaryForm
-> Maybe DictionaryForm
forall {p} {p}. p -> p -> p
trace (DictionaryForm -> String
forall a. Show a => a -> String
show DictionaryForm
dictionary)
            DictionaryForm -> Maybe DictionaryForm
simplexPivot (Integer -> Integer -> DictionaryForm -> DictionaryForm
pivot Integer
pivotBasicVar Integer
pivotNonBasicVar DictionaryForm
dictionary )
  where
    ratioTest :: DictionaryForm -> Integer -> Maybe Integer -> Maybe Rational -> Maybe Integer
    ratioTest :: DictionaryForm
-> Integer -> Maybe Integer -> Maybe Rational -> Maybe Integer
ratioTest []                    Integer
_               Maybe Integer
mCurrentMinBasicVar Maybe Rational
_           = Maybe Integer
mCurrentMinBasicVar
    ratioTest ((Integer
basicVar, VarConstMap
lp) : DictionaryForm
xs) Integer
mostNegativeVar Maybe Integer
mCurrentMinBasicVar Maybe Rational
mCurrentMin =
      case Integer -> VarConstMap -> Maybe Rational
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Integer
mostNegativeVar VarConstMap
lp of
        Maybe Rational
Nothing                         -> DictionaryForm
-> Integer -> Maybe Integer -> Maybe Rational -> Maybe Integer
ratioTest DictionaryForm
xs Integer
mostNegativeVar Maybe Integer
mCurrentMinBasicVar Maybe Rational
mCurrentMin
        Just Rational
currentCoeff ->
          let 
            rhs :: Rational
rhs = Rational -> Maybe Rational -> Rational
forall a. a -> Maybe a -> a
fromMaybe Rational
0 (Integer -> VarConstMap -> Maybe Rational
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup (-Integer
1) VarConstMap
lp)
          in
            if Rational
currentCoeff Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
0 Bool -> Bool -> Bool
|| Rational
rhs Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
< Rational
0
              then 
                -- trace (show currentCoeff)
                DictionaryForm
-> Integer -> Maybe Integer -> Maybe Rational -> Maybe Integer
ratioTest DictionaryForm
xs Integer
mostNegativeVar Maybe Integer
mCurrentMinBasicVar Maybe Rational
mCurrentMin -- rhs was already in right side in original tableau, so should be above zero
                                                                              -- Coeff needs to be negative since it has been moved to the RHS
              else
                case Maybe Rational
mCurrentMin of
                  Maybe Rational
Nothing         -> DictionaryForm
-> Integer -> Maybe Integer -> Maybe Rational -> Maybe Integer
ratioTest DictionaryForm
xs Integer
mostNegativeVar (Integer -> Maybe Integer
forall a. a -> Maybe a
Just Integer
basicVar) (Rational -> Maybe Rational
forall a. a -> Maybe a
Just (Rational
rhs Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
currentCoeff))
                  Just Rational
currentMin ->
                    if (Rational
rhs Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
currentCoeff) Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
currentMin
                      then DictionaryForm
-> Integer -> Maybe Integer -> Maybe Rational -> Maybe Integer
ratioTest DictionaryForm
xs Integer
mostNegativeVar (Integer -> Maybe Integer
forall a. a -> Maybe a
Just Integer
basicVar) (Rational -> Maybe Rational
forall a. a -> Maybe a
Just (Rational
rhs Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational
currentCoeff))
                      else DictionaryForm
-> Integer -> Maybe Integer -> Maybe Rational -> Maybe Integer
ratioTest DictionaryForm
xs Integer
mostNegativeVar Maybe Integer
mCurrentMinBasicVar Maybe Rational
mCurrentMin

    mostPositive :: (Integer, VarConstMap) -> Maybe Integer
    mostPositive :: (Integer, VarConstMap) -> Maybe Integer
mostPositive (Integer
_, VarConstMap
lp) = 
      case VarConstMap
-> Maybe (Integer, Rational) -> Maybe (Integer, Rational)
findLargestCoeff VarConstMap
lp Maybe (Integer, Rational)
forall a. Maybe a
Nothing of
        Just (Integer
largestVar, Rational
largestCoeff) ->
          if Rational
largestCoeff Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
0 
            then Maybe Integer
forall a. Maybe a
Nothing
            else Integer -> Maybe Integer
forall a. a -> Maybe a
Just Integer
largestVar
        Maybe (Integer, Rational)
Nothing -> String -> Maybe Integer -> Maybe Integer
forall {p} {p}. p -> p -> p
trace String
"No variables in first row when looking for most positive" Maybe Integer
forall a. Maybe a
Nothing

      where
        findLargestCoeff :: VarConstMap -> Maybe (Integer, Rational) -> Maybe (Integer, Rational)
        findLargestCoeff :: VarConstMap
-> Maybe (Integer, Rational) -> Maybe (Integer, Rational)
findLargestCoeff [] Maybe (Integer, Rational)
mCurrentMax                  = Maybe (Integer, Rational)
mCurrentMax
        findLargestCoeff ((Integer
var, Rational
coeff) : VarConstMap
xs) Maybe (Integer, Rational)
mCurrentMax = 
          if Integer
var Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== (-Integer
1) 
            then VarConstMap
-> Maybe (Integer, Rational) -> Maybe (Integer, Rational)
findLargestCoeff VarConstMap
xs Maybe (Integer, Rational)
mCurrentMax
            else 
              case Maybe (Integer, Rational)
mCurrentMax of
                Maybe (Integer, Rational)
Nothing         -> VarConstMap
-> Maybe (Integer, Rational) -> Maybe (Integer, Rational)
findLargestCoeff VarConstMap
xs ((Integer, Rational) -> Maybe (Integer, Rational)
forall a. a -> Maybe a
Just (Integer
var, Rational
coeff))
                Just (Integer, Rational)
currentMax ->
                  if (Integer, Rational) -> Rational
forall a b. (a, b) -> b
snd (Integer, Rational)
currentMax Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
>= Rational
coeff 
                    then VarConstMap
-> Maybe (Integer, Rational) -> Maybe (Integer, Rational)
findLargestCoeff VarConstMap
xs Maybe (Integer, Rational)
mCurrentMax
                    else VarConstMap
-> Maybe (Integer, Rational) -> Maybe (Integer, Rational)
findLargestCoeff VarConstMap
xs ((Integer, Rational) -> Maybe (Integer, Rational)
forall a. a -> Maybe a
Just (Integer
var, Rational
coeff))

    -- |Pivot a dictionary using the two given variables.
    -- The first variable is the leaving (non-basic) variable.
    -- The second variable is the entering (basic) variable.
    -- Expects the entering variable to be present in the row containing the leaving variable.
    -- Expects each row to have a unique basic variable.
    -- Expects each basic variable to not appear on the RHS of any equation.
    pivot :: Integer -> Integer -> DictionaryForm -> DictionaryForm
    pivot :: Integer -> Integer -> DictionaryForm -> DictionaryForm
pivot Integer
leavingVariable Integer
enteringVariable DictionaryForm
rows =
      case Integer -> VarConstMap -> Maybe Rational
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Integer
enteringVariable VarConstMap
basicRow of
        Just Rational
nonBasicCoeff ->
          DictionaryForm
updatedRows
          where
            -- Move entering variable to basis, update other variables in row appropriately
            pivotEquation :: (Integer, VarConstMap)
pivotEquation = (Integer
enteringVariable, ((Integer, Rational) -> (Integer, Rational))
-> VarConstMap -> VarConstMap
forall a b. (a -> b) -> [a] -> [b]
map ((Rational -> Rational)
-> (Integer, Rational) -> (Integer, Rational)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second (Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/ Rational -> Rational
forall a. Num a => a -> a
negate Rational
nonBasicCoeff)) ((Integer
leavingVariable, -Rational
1) (Integer, Rational) -> VarConstMap -> VarConstMap
forall a. a -> [a] -> [a]
: ((Integer, Rational) -> Bool) -> VarConstMap -> VarConstMap
forall a. (a -> Bool) -> [a] -> [a]
filter ((Integer
enteringVariable Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/=) (Integer -> Bool)
-> ((Integer, Rational) -> Integer) -> (Integer, Rational) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer, Rational) -> Integer
forall a b. (a, b) -> a
fst) VarConstMap
basicRow))
            -- Substitute pivot equation into other rows
            updatedRows :: DictionaryForm
updatedRows =
              ((Integer, VarConstMap) -> (Integer, VarConstMap))
-> DictionaryForm -> DictionaryForm
forall a b. (a -> b) -> [a] -> [b]
map
              (\(Integer
basicVar, VarConstMap
vMap) ->
                if Integer
leavingVariable Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
basicVar
                  then (Integer, VarConstMap)
pivotEquation
                  else
                    case Integer -> VarConstMap -> Maybe Rational
forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Integer
enteringVariable VarConstMap
vMap of
                      Just Rational
subsCoeff -> (Integer
basicVar, (VarConstMap -> VarConstMap
foldSumVarConstMap (VarConstMap -> VarConstMap)
-> (VarConstMap -> VarConstMap) -> VarConstMap -> VarConstMap
forall b c a. (b -> c) -> (a -> b) -> a -> c
. VarConstMap -> VarConstMap
forall a. Ord a => [a] -> [a]
sort) (((Integer, Rational) -> (Integer, Rational))
-> VarConstMap -> VarConstMap
forall a b. (a -> b) -> [a] -> [b]
map ((Rational -> Rational)
-> (Integer, Rational) -> (Integer, Rational)
forall (p :: * -> * -> *) b c a.
Bifunctor p =>
(b -> c) -> p a b -> p a c
second (Rational
subsCoeff Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
*)) ((Integer, VarConstMap) -> VarConstMap
forall a b. (a, b) -> b
snd (Integer, VarConstMap)
pivotEquation) VarConstMap -> VarConstMap -> VarConstMap
forall a. [a] -> [a] -> [a]
++ ((Integer, Rational) -> Bool) -> VarConstMap -> VarConstMap
forall a. (a -> Bool) -> [a] -> [a]
filter ((Integer
enteringVariable Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/=) (Integer -> Bool)
-> ((Integer, Rational) -> Integer) -> (Integer, Rational) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer, Rational) -> Integer
forall a b. (a, b) -> a
fst) VarConstMap
vMap))
                      Maybe Rational
Nothing -> (Integer
basicVar, VarConstMap
vMap)
              )
              DictionaryForm
rows
        Maybe Rational
Nothing -> String -> DictionaryForm -> DictionaryForm
forall {p} {p}. p -> p -> p
trace String
"non basic variable not found in basic row" DictionaryForm
forall a. HasCallStack => a
undefined
      where
        (Integer
_, VarConstMap
basicRow) = DictionaryForm -> (Integer, VarConstMap)
forall a. [a] -> a
head (DictionaryForm -> (Integer, VarConstMap))
-> DictionaryForm -> (Integer, VarConstMap)
forall a b. (a -> b) -> a -> b
$ ((Integer, VarConstMap) -> Bool)
-> DictionaryForm -> DictionaryForm
forall a. (a -> Bool) -> [a] -> [a]
filter ((Integer
leavingVariable Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==) (Integer -> Bool)
-> ((Integer, VarConstMap) -> Integer)
-> (Integer, VarConstMap)
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer, VarConstMap) -> Integer
forall a b. (a, b) -> a
fst) DictionaryForm
rows