diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,3 +1,20 @@ # Changelog for simplex-haskell ## Unreleased changes + +## [v0.2.0.0](https://github.com/rasheedja/LPPaver/tree/v0.2.0.0) + +- Setup CI +- Use fourmolu formatter +- Add better types +- Use lens +- Use RecordDot syntax +- Add logging +- Improve Docs +- More Tests +- Bump Stackage LTS +- Rename Linear.Simplex.Simplex -> Linear.Simplex.TwoPhase.Simplex + +## [v0.1.0.0](https://github.com/rasheedja/LPPaver/tree/v0.1.0.0) + +- Initial release diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright Junaid Rasheed (c) 2020-2022 +Copyright Junaid Rasheed (c) 2020-2023 All rights reserved. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -4,14 +4,14 @@ ## Quick Overview -The `Linear.Simplex.Simplex` module contain both phases of the simplex method. +The `Linear.Simplex.Solver.TwoPhase` module contain both phases of the two-phase simplex method. ### Phase One Phase one is implemented by `findFeasibleSolution`: ```haskell -findFeasibleSolution :: [PolyConstraint] -> Maybe (DictionaryForm, [Integer], [Integer], Integer) +findFeasibleSolution :: (MonadIO m, MonadLogger m) => [PolyConstraint] -> m (Maybe FeasibleSystem) ``` `findFeasibleSolution` takes a list of `PolyConstraint`s. @@ -19,85 +19,97 @@ `PolyConstraint` is defined as: ```haskell -data PolyConstraint = - LEQ VarConstMap Rational | - GEQ VarConstMap Rational | - EQ VarConstMap Rational deriving (Show, Eq); +data PolyConstraint + = LEQ {lhs :: VarLitMapSum, rhs :: SimplexNum} + | GEQ {lhs :: VarLitMapSum, rhs :: SimplexNum} + | EQ {lhs :: VarLitMapSum, rhs :: SimplexNum} + deriving (Show, Read, Eq, Generic) ``` -And `VarConstMap` is defined as: +`SimplexNum` is an alias for `Rational`, and `VarLitMapSum` is an alias for `VarLitMap`, which is an alias for `Map Var SimplexNum`. +`Var` is an alias of `Int`. -```haskell -type VarConstMap = [(Integer, Rational)] -``` +A `VarLitMapSum` is read as `Integer` variables mapped to their `Rational` coefficients, with an implicit `+` between each entry. +For example: `Map.fromList [(1, 2), (2, (-3)), (1, 3)]` is equivalent to `(2x1 + (-3x2) + 3x1)`. -A `VarConstMap` is treated as a list of `Integer` variables mapped to their `Rational` coefficients, with an implicit `+` between each element in the list. -For example: `[(1, 2), (2, (-3)), (1, 3)]` is equivalent to `(2x1 + (-3x2) + 3x1)`. +And a `PolyConstraint` is an inequality/equality where the LHS is a `VarLitMapSum` and the RHS is a `Rational`. +For example: `LEQ (Map.fromList [(1, 2), (2, (-3)), (1, 3)] 60)` is equivalent to `(2x1 + (-3x2) + 3x1) <= 60`. -And a `PolyConstraint` is an inequality/equality where the LHS is a `VarConstMap` and the RHS is a `Rational`. -For example: `LEQ [(1, 2), (2, (-3)), (1, 3)] 60` is equivalent to `(2x1 + (-3x2) + 3x1) <= 60`. +Passing a `[PolyConstraint]` to `findFeasibleSolution` will return a `FeasibleSystem` if a feasible solution exists: -Passing a `[PolyConstraint]` to `findFeasibleSolution` will return a feasible solution if it exists as well as a list of slack variables, artificial variables, and a variable that can be safely used to represent the objective for phase two. -`Nothing` is returned if the given `[PolyConstraint]` is infeasible. -The feasible system is returned as the type `DictionaryForm`: +```haskell +data FeasibleSystem = FeasibleSystem + { dict :: Dict + , slackVars :: [Var] + , artificialVars :: [Var] + , objectiveVar :: Var + } + deriving (Show, Read, Eq, Generic) +``` ```haskell -type DictionaryForm = [(Integer, VarConstMap)] +type Dict = M.Map Var DictValue + +data DictValue = DictValue + { varMapSum :: VarLitMapSum + , constant :: SimplexNum + } + deriving (Show, Read, Eq, Generic) ``` -`DictionaryForm` can be thought of as a list of equations, where the `Integer` represents a basic variable on the LHS that is equal to the RHS represented as a `VarConstMap`. In this `VarConstMap`, the `Integer` -1 is used internally to represent a `Rational` number. +`Dict` can be thought of as a set of equations, where the key represents a basic variable on the LHS of the equation +that is equal to the RHS represented as a `DictValue` value. ### Phase Two `optimizeFeasibleSystem` performs phase two of the simplex method, and has the type: ```haskell -data ObjectiveFunction = Max VarConstMap | Min VarConstMap deriving (Show, Eq) -optimizeFeasibleSystem :: ObjectiveFunction -> DictionaryForm -> [Integer] -> [Integer] -> Integer -> Maybe (Integer, [(Integer, Rational)]) +optimizeFeasibleSystem :: (MonadIO m, MonadLogger m) => ObjectiveFunction -> FeasibleSystem -> m (Maybe Result) + +data ObjectiveFunction = Max {objective :: VarLitMapSum} | Min {objective :: VarLitMapSum} + +data Result = Result + { objectiveVar :: Var + , varValMap :: VarLitMap + } + deriving (Show, Read, Eq, Generic) ``` -We first pass an `ObjectiveFunction`. -Then we give a feasible system in `DictionaryForm`, a list of slack variables, a list of artificial variables, and a variable to represent the objective. -`optimizeFeasibleSystem` Maximizes/Minimizes the linear equation represented as a `VarConstMap` in the given `ObjectiveFunction`. -The first item of the returned pair is the `Integer` variable representing the objective. -The second item is a list of `Integer` variables mapped to their optimized values. -If a variable is not in this list, the variable is equal to 0. +We give `optimizeFeasibleSystem` an `ObjectiveFunction` along with a `FeasibleSystem`. ### Two-Phase Simplex + `twoPhaseSimplex` performs both phases of the simplex method. It has the type: + ```haskell -twoPhaseSimplex :: ObjectiveFunction -> [PolyConstraint] -> Maybe (Integer, [(Integer, Rational)]) +twoPhaseSimplex :: (MonadIO m, MonadLogger m) => ObjectiveFunction -> [PolyConstraint] -> m (Maybe Result) ``` -The return type is the same as that of `optimizeFeasibleSystem` ### Extracting Results -The result of the objective function is present in the return type of both `twoPhaseSimplex` and `optimizeFeasibleSystem`, but this can be difficult to grok in systems with many variables, so the following function will extract the value of the objective function for you. +The result of the objective function is present in the returned `Result` type of both `twoPhaseSimplex` and `optimizeFeasibleSystem`, but this can be difficult to grok in systems with many variables, so the following function will extract the value of the objective function for you. + ```haskell -extractObjectiveValue :: Maybe (Integer, [(Integer, Rational)]) -> Maybe Rational +dictionaryFormToTableau :: Dict -> Tableau ``` There are similar functions for `DictionaryForm` as well as other custom types in the module `Linear.Simplex.Util`. -## Usage notes - -You must only use positive `Integer` variables in a `VarConstMap`. -This implementation assumes that the user only provides positive `Integer` variables; the `Integer` -1, for example, is sometimes used to represent a `Rational` number. - ## Example ```haskell exampleFunction :: (ObjectiveFunction, [PolyConstraint]) exampleFunction = ( - Max [(1, 3), (2, 5)], -- 3x1 + 5x2 + Max {objective = Map.fromList [(1, 3), (2, 5)]}, -- 3x1 + 5x2 [ - LEQ [(1, 3), (2, 1)] 15, -- 3x1 + x2 <= 15 - LEQ [(1, 1), (2, 1)] 7, -- x1 + x2 <= 7 - LEQ [(2, 1)] 4, -- x2 <= 4 - LEQ [(1, -1), (2, 2)] 6 -- -x1 + 2x2 <= 6 + LEQ {lhs = Map.fromList [(1, 3), (2, 1)], rhs = 15}, -- 3x1 + x2 <= 15 + LEQ {lhs = Map.fromList [(1, 1), (2, 1)], rhs = 7}, -- x1 + x2 <= 7 + LEQ {lhs = Map.fromList [(2, 1)], rhs = 4}, -- x2 <= 4 + LEQ {lhs = Map.fromList [(1, -1), (2, 2)], rhs = 6} -- -x1 + 2x2 <= 6 ] ) @@ -105,14 +117,17 @@ ``` The result of the call above is: + ```haskell -Just - (7, -- Integer representing objective function - [ - (7,29 % 1), -- Value for variable 7, so max(3x1 + 5x2) = 29. - (1,3 % 1), -- Value for variable 1, so x1 = 3 - (2,4 % 1) -- Value for variable 2, so x2 = 4 - ] +Just + (Result + { objectiveVar = 7 -- Integer representing objective function + , varValMap = Map.fromList + [ (7, 29) -- Value for variable 7, so max(3x1 + 5x2) = 29. + , (1, 3) -- Value for variable 1, so x1 = 3 + , (2, 4) -- Value for variable 2, so x2 = 4 + ] + } ) ``` diff --git a/Setup.hs b/Setup.hs deleted file mode 100644 --- a/Setup.hs +++ /dev/null @@ -1,2 +0,0 @@ -import Distribution.Simple -main = defaultMain diff --git a/simplex-method.cabal b/simplex-method.cabal --- a/simplex-method.cabal +++ b/simplex-method.cabal @@ -1,11 +1,11 @@ cabal-version: 1.12 --- This file has been generated from package.yaml by hpack version 0.34.4. +-- This file has been generated from package.yaml by hpack version 0.36.0. -- -- see: https://github.com/sol/hpack name: simplex-method -version: 0.1.0.0 +version: 0.2.0.0 synopsis: Implementation of the two-phase simplex method in exact rational arithmetic description: Please see the README on GitHub at category: Math, Maths, Mathematics, Optimisation, Optimization, Linear Programming @@ -28,15 +28,23 @@ library exposed-modules: Linear.Simplex.Prettify - Linear.Simplex.Simplex + Linear.Simplex.Solver.TwoPhase Linear.Simplex.Types Linear.Simplex.Util other-modules: Paths_simplex_method hs-source-dirs: src + default-extensions: + DataKinds DeriveFunctor DeriveGeneric DisambiguateRecordFields DuplicateRecordFields FlexibleContexts LambdaCase OverloadedLabels OverloadedRecordDot OverloadedStrings RecordWildCards TemplateHaskell TupleSections TypeApplications NamedFieldPuns build-depends: - base >=4.7 && <5 + base >=4.14 && <5 + , containers >=0.6.5.1 && <0.7 + , generic-lens >=2.2.0 && <2.3 + , lens >=5.2.2 && <5.3 + , monad-logger >=0.3.40 && <0.4 + , text >=2.0.2 && <2.1 + , time >=1.12.2 && <1.13 default-language: Haskell2010 test-suite simplex-haskell-test @@ -47,7 +55,15 @@ Paths_simplex_method hs-source-dirs: test + default-extensions: + DataKinds DeriveFunctor DeriveGeneric DisambiguateRecordFields DuplicateRecordFields FlexibleContexts LambdaCase OverloadedLabels OverloadedRecordDot OverloadedStrings RecordWildCards TemplateHaskell TupleSections TypeApplications NamedFieldPuns build-depends: - base >=4.7 && <5 + base >=4.14 && <5 + , containers >=0.6.5.1 && <0.7 + , generic-lens >=2.2.0 && <2.3 + , lens >=5.2.2 && <5.3 + , monad-logger >=0.3.40 && <0.4 , simplex-method + , text >=2.0.2 && <2.1 + , time >=1.12.2 && <1.13 default-language: Haskell2010 diff --git a/src/Linear/Simplex/Prettify.hs b/src/Linear/Simplex/Prettify.hs --- a/src/Linear/Simplex/Prettify.hs +++ b/src/Linear/Simplex/Prettify.hs @@ -1,39 +1,44 @@ -{-| -Module : Linear.Simplex.Prettify -Description : Prettifier for "Linear.Simplex.Types" types -Copyright : (c) Junaid Rasheed, 2020-2022 -License : BSD-3 -Maintainer : jrasheed178@gmail.com -Stability : experimental +{-# LANGUAGE ImportQualifiedPost #-} +{-# LANGUAGE RankNTypes #-} -Converts "Linear.Simplex.Types" types into human-readable 'String's --} +-- | +-- Module : Linear.Simplex.Prettify +-- Description : Prettifier for "Linear.Simplex.Types" types +-- Copyright : (c) Junaid Rasheed, 2020-2023 +-- License : BSD-3 +-- Maintainer : jrasheed178@gmail.com +-- Stability : experimental +-- +-- Converts "Linear.Simplex.Types" types into human-readable 'String's module Linear.Simplex.Prettify where -import Linear.Simplex.Types as T +import Control.Lens +import Data.Generics.Labels () +import Data.Map qualified as M import Data.Ratio +import Linear.Simplex.Types --- |Convert a 'VarConstMap' into a human-readable 'String' -prettyShowVarConstMap :: VarConstMap -> String -prettyShowVarConstMap [] = "" -prettyShowVarConstMap [(v, c)] = prettyShowRational c ++ " * x" ++ show v ++ "" +-- | Convert a 'VarConstMap' into a human-readable 'String' +prettyShowVarConstMap :: VarLitMapSum -> String +prettyShowVarConstMap = aux . M.toList where - prettyShowRational r = - if r < 0 - then "(" ++ r' ++ ")" - else r' + aux [] = "" + aux ((vName, vCoeff) : vs) = prettyShowRational vCoeff ++ " * " ++ show vName ++ " + " ++ aux vs where - r' = if denominator r == 1 then show (numerator r) else show (numerator r) ++ " / " ++ show (numerator r) - -prettyShowVarConstMap ((v, c) : vcs) = prettyShowVarConstMap [(v, c)] ++ " + " ++ prettyShowVarConstMap vcs + prettyShowRational r = + if r < 0 + then "(" ++ r' ++ ")" + else r' + where + r' = if denominator r == 1 then show (numerator r) else show (numerator r) ++ " / " ++ show (numerator r) --- |Convert a 'PolyConstraint' into a human-readable 'String' +-- | Convert a 'PolyConstraint' into a human-readable 'String' prettyShowPolyConstraint :: PolyConstraint -> String prettyShowPolyConstraint (LEQ vcm r) = prettyShowVarConstMap vcm ++ " <= " ++ show r prettyShowPolyConstraint (GEQ vcm r) = prettyShowVarConstMap vcm ++ " >= " ++ show r -prettyShowPolyConstraint (T.EQ vcm r) = prettyShowVarConstMap vcm ++ " == " ++ show r +prettyShowPolyConstraint (Linear.Simplex.Types.EQ vcm r) = prettyShowVarConstMap vcm ++ " == " ++ show r --- |Convert an 'ObjectiveFunction' into a human-readable 'String' +-- | Convert an 'ObjectiveFunction' into a human-readable 'String' prettyShowObjectiveFunction :: ObjectiveFunction -> String prettyShowObjectiveFunction (Min vcm) = "min: " ++ prettyShowVarConstMap vcm prettyShowObjectiveFunction (Max vcm) = "max: " ++ prettyShowVarConstMap vcm diff --git a/src/Linear/Simplex/Simplex.hs b/src/Linear/Simplex/Simplex.hs deleted file mode 100644 --- a/src/Linear/Simplex/Simplex.hs +++ /dev/null @@ -1,289 +0,0 @@ -{-# 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 s a = 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 unsimplifiedSystem = - if null artificialVars -- No artificial vars, we have a feasible system - then Just (systemWithBasicVarsAsDictionary, slackVars, artificialVars, objectiveVar) - else - case simplexPivot (createObjectiveDict artificialObjective objectiveVar : systemWithBasicVarsAsDictionary) of - Just phase1Dict -> - let - eliminateArtificialVarsFromPhase1Tableau = map (second (filter (\(v, _) -> v `notElem` artificialVars))) phase1Dict - in - case lookup objectiveVar eliminateArtificialVarsFromPhase1Tableau of - Nothing -> trace "objective row not found in phase 1 tableau" Nothing -- Should this be an error? - Just row -> - if fromMaybe 0 (lookup (-1) row) == 0 - then Just (eliminateArtificialVarsFromPhase1Tableau, slackVars, artificialVars, objectiveVar) - else trace "rhs not zero after phase 1, thus original tableau is infeasible" Nothing - Nothing -> Nothing - where - system = simplifySystem unsimplifiedSystem - - maxVar = - maximum $ map - (\case - LEQ vcm _ -> maximum (map fst vcm) - GEQ vcm _ -> maximum (map fst vcm) - EQ vcm _ -> maximum (map fst vcm) - ) - system - - (systemWithSlackVars, slackVars) = systemInStandardForm system maxVar [] - - maxVarWithSlackVars = if null slackVars then maxVar else maximum slackVars - - (systemWithBasicVars, artificialVars) = systemWithArtificialVars systemWithSlackVars maxVarWithSlackVars - - finalMaxVar = if null artificialVars then maxVarWithSlackVars else maximum artificialVars - - systemWithBasicVarsAsDictionary = tableauInDictionaryForm systemWithBasicVars - - artificialObjective = createArtificialObjective systemWithBasicVarsAsDictionary artificialVars - - objectiveVar = finalMaxVar + 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 [] _ sVars = ([], sVars) - systemInStandardForm (EQ v r : xs) maxVar sVars = ((Nothing, EQ v r) : newSystem, newSlackVars) - where - (newSystem, newSlackVars) = systemInStandardForm xs maxVar sVars - systemInStandardForm (LEQ v r : xs) maxVar sVars = ((Just newSlackVar, EQ (v ++ [(newSlackVar, 1)]) r) : newSystem, newSlackVars) - where - newSlackVar = maxVar + 1 - (newSystem, newSlackVars) = systemInStandardForm xs newSlackVar (newSlackVar : sVars) - systemInStandardForm (GEQ v r : xs) maxVar sVars = ((Just newSlackVar, EQ (v ++ [(newSlackVar, -1)]) r) : newSystem, newSlackVars) - where - newSlackVar = maxVar + 1 - (newSystem, newSlackVars) = systemInStandardForm xs newSlackVar (newSlackVar : 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 [] _ = ([],[]) - systemWithArtificialVars ((mVar, EQ v r) : pcs) maxVar = - case mVar of - Nothing -> - if r >= 0 - then - ((newArtificialVar, (v ++ [(newArtificialVar, 1)], r)) : newSystemWithNewMaxVar, newArtificialVar : artificialVarsWithNewMaxVar) - else - ((newArtificialVar, (v ++ [(newArtificialVar, -1)], r)) : newSystemWithNewMaxVar, newArtificialVar : artificialVarsWithNewMaxVar) - Just basicVar -> - case lookup basicVar v of - Just basicVarCoeff -> - if r == 0 - then ((basicVar, (v, r)) : newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) - else - if r > 0 - then - if basicVarCoeff >= 0 -- Should only be 1 in the standard call path - then ((basicVar, (v, r)) : newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) - else ((newArtificialVar, (v ++ [(newArtificialVar, 1)], r)) : newSystemWithNewMaxVar, newArtificialVar : artificialVarsWithNewMaxVar) -- Slack var is negative, r is positive (when original constraint was GEQ) - else -- r < 0 - if basicVarCoeff <= 0 -- Should only be -1 in the standard call path - then ((basicVar, (v, r)) : newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) - else ((newArtificialVar, (v ++ [(newArtificialVar, -1)], r)) : newSystemWithNewMaxVar, newArtificialVar : artificialVarsWithNewMaxVar) -- Slack var is negative, r is negative (when original constraint was LEQ) - where - newArtificialVar = maxVar + 1 - - (newSystemWithNewMaxVar, artificialVarsWithNewMaxVar) = systemWithArtificialVars pcs newArtificialVar - - (newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) = systemWithArtificialVars pcs 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 rows artificialVars = Max negatedSumWithoutArtificialVars - where - rowsToAdd = filter (\(i, _) -> i `elem` artificialVars) rows - negatedRows = map (\(_, vcm) -> map (second negate) vcm) rowsToAdd - negatedSum = foldSumVarConstMap ((sort . concat) negatedRows) - negatedSumWithoutArtificialVars = filter (\(v, _) -> v `notElem` artificialVars) 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 unsimplifiedObjFunction phase1Dict slackVars artificialVars objectiveVar = - if null artificialVars - then displayResults . dictionaryFormToTableau <$> simplexPivot (createObjectiveDict objFunction objectiveVar : phase1Dict) - else displayResults . dictionaryFormToTableau <$> simplexPivot (createObjectiveDict phase2ObjFunction objectiveVar : tail phase1Dict) - where - objFunction = simplifyObjectiveFunction unsimplifiedObjFunction - - displayResults :: Tableau -> (Integer, [(Integer, Rational)]) - displayResults tableau = - ( - objectiveVar, - case objFunction of - Max _ -> - map - (second snd) - $ filter (\(basicVar,_) -> basicVar `notElem` slackVars ++ artificialVars) tableau - Min _ -> - map -- We maximized -objVar, so we negate the objVar to get the final value - (\(basicVar, row) -> if basicVar == objectiveVar then (basicVar, negate (snd row)) else (basicVar, snd row)) - $ filter (\(basicVar,_) -> basicVar `notElem` slackVars ++ artificialVars) tableau - ) - - phase2Objective = - (foldSumVarConstMap . sort) $ - concatMap - (\(var, coeff) -> - case lookup var phase1Dict of - Nothing -> [(var, coeff)] - Just row -> map (second (*coeff)) row - ) - (getObjective objFunction) - - phase2ObjFunction = if isMax objFunction then Max phase2Objective else Min 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 objFunction unsimplifiedSystem = - case findFeasibleSolution unsimplifiedSystem of - Just r@(phase1Dict, slackVars, artificialVars, objectiveVar) -> optimizeFeasibleSystem objFunction phase1Dict slackVars artificialVars objectiveVar - Nothing -> 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 dictionary = - trace (show dictionary) $ - case mostPositive (head dictionary) of - Nothing -> - trace "all neg \n" - trace (show dictionary) - Just dictionary - Just pivotNonBasicVar -> - let - mPivotBasicVar = ratioTest (tail dictionary) pivotNonBasicVar Nothing Nothing - in - case mPivotBasicVar of - Nothing -> trace ("Ratio test failed on non-basic var: " ++ show pivotNonBasicVar ++ "\n" ++ show dictionary) Nothing - Just pivotBasicVar -> - trace "one pos \n" - trace (show dictionary) - simplexPivot (pivot pivotBasicVar pivotNonBasicVar dictionary ) - where - ratioTest :: DictionaryForm -> Integer -> Maybe Integer -> Maybe Rational -> Maybe Integer - ratioTest [] _ mCurrentMinBasicVar _ = mCurrentMinBasicVar - ratioTest ((basicVar, lp) : xs) mostNegativeVar mCurrentMinBasicVar mCurrentMin = - case lookup mostNegativeVar lp of - Nothing -> ratioTest xs mostNegativeVar mCurrentMinBasicVar mCurrentMin - Just currentCoeff -> - let - rhs = fromMaybe 0 (lookup (-1) lp) - in - if currentCoeff >= 0 || rhs < 0 - then - -- trace (show currentCoeff) - ratioTest xs mostNegativeVar mCurrentMinBasicVar 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 mCurrentMin of - Nothing -> ratioTest xs mostNegativeVar (Just basicVar) (Just (rhs / currentCoeff)) - Just currentMin -> - if (rhs / currentCoeff) >= currentMin - then ratioTest xs mostNegativeVar (Just basicVar) (Just (rhs / currentCoeff)) - else ratioTest xs mostNegativeVar mCurrentMinBasicVar mCurrentMin - - mostPositive :: (Integer, VarConstMap) -> Maybe Integer - mostPositive (_, lp) = - case findLargestCoeff lp Nothing of - Just (largestVar, largestCoeff) -> - if largestCoeff <= 0 - then Nothing - else Just largestVar - Nothing -> trace "No variables in first row when looking for most positive" Nothing - - where - findLargestCoeff :: VarConstMap -> Maybe (Integer, Rational) -> Maybe (Integer, Rational) - findLargestCoeff [] mCurrentMax = mCurrentMax - findLargestCoeff ((var, coeff) : xs) mCurrentMax = - if var == (-1) - then findLargestCoeff xs mCurrentMax - else - case mCurrentMax of - Nothing -> findLargestCoeff xs (Just (var, coeff)) - Just currentMax -> - if snd currentMax >= coeff - then findLargestCoeff xs mCurrentMax - else findLargestCoeff xs (Just (var, 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 leavingVariable enteringVariable rows = - case lookup enteringVariable basicRow of - Just nonBasicCoeff -> - updatedRows - where - -- Move entering variable to basis, update other variables in row appropriately - pivotEquation = (enteringVariable, map (second (/ negate nonBasicCoeff)) ((leavingVariable, -1) : filter ((enteringVariable /=) . fst) basicRow)) - -- Substitute pivot equation into other rows - updatedRows = - map - (\(basicVar, vMap) -> - if leavingVariable == basicVar - then pivotEquation - else - case lookup enteringVariable vMap of - Just subsCoeff -> (basicVar, (foldSumVarConstMap . sort) (map (second (subsCoeff *)) (snd pivotEquation) ++ filter ((enteringVariable /=) . fst) vMap)) - Nothing -> (basicVar, vMap) - ) - rows - Nothing -> trace "non basic variable not found in basic row" undefined - where - (_, basicRow) = head $ filter ((leavingVariable ==) . fst) rows diff --git a/src/Linear/Simplex/Solver/TwoPhase.hs b/src/Linear/Simplex/Solver/TwoPhase.hs new file mode 100644 --- /dev/null +++ b/src/Linear/Simplex/Solver/TwoPhase.hs @@ -0,0 +1,570 @@ +-- | +-- Module : Linear.Simplex.Simplex.TwoPhase +-- Description : Implements the twoPhaseSimplex method +-- Copyright : (c) Junaid Rasheed, 2020-2023 +-- 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.Solver.TwoPhase (findFeasibleSolution, optimizeFeasibleSystem, twoPhaseSimplex) where + +import Prelude hiding (EQ) + +import Control.Lens +import Control.Monad (unless) +import Control.Monad.IO.Class (MonadIO) +import Control.Monad.Logger +import Data.Bifunctor +import Data.List +import qualified Data.Map as M +import Data.Maybe (fromJust, fromMaybe, mapMaybe) +import Data.Ratio (denominator, numerator, (%)) +import qualified Data.Text as Text +import GHC.Real (Ratio) +import Linear.Simplex.Types +import Linear.Simplex.Util + +-- | Find a feasible solution for the given system of 'PolyConstraint's by performing the first phase of the two-phase simplex method +-- All variables in the 'PolyConstraint' must be positive. +-- If the system is infeasible, return 'Nothing' +-- Otherwise, return the feasible system in 'Dict' as well as a list of slack variables, a list artificial variables, and the objective variable. +findFeasibleSolution :: (MonadIO m, MonadLogger m) => [PolyConstraint] -> m (Maybe FeasibleSystem) +findFeasibleSolution unsimplifiedSystem = do + logMsg LevelInfo $ "findFeasibleSolution: Looking for solution for " <> showT unsimplifiedSystem + if null artificialVars -- No artificial vars, we have a feasible system + then do + logMsg LevelInfo "findFeasibleSolution: Feasible solution found with no artificial vars" + pure . Just $ FeasibleSystem systemWithBasicVarsAsDictionary slackVars artificialVars objectiveVar + else do + logMsg LevelInfo $ + "findFeasibleSolution: Needed to create artificial vars. System with artificial vars (in Tableau form) " + <> showT systemWithBasicVars + mPhase1Dict <- simplexPivot artificialPivotObjective systemWithBasicVarsAsDictionary + case mPhase1Dict of + Just phase1Dict -> do + logMsg LevelInfo $ + "findFeasibleSolution: System after pivoting with objective" + <> showT artificialPivotObjective + <> ": " + <> showT phase1Dict + let eliminateArtificialVarsFromPhase1Tableau = + M.map + ( \DictValue {..} -> + DictValue + { varMapSum = M.filterWithKey (\k _ -> k `notElem` artificialVars) varMapSum + , .. + } + ) + phase1Dict + case M.lookup objectiveVar eliminateArtificialVarsFromPhase1Tableau of + Nothing -> do + logMsg LevelWarn $ + "findFeasibleSolution: Objective row not found after eliminatiing artificial vars. This is unexpected. System without artificial vars (in Dict form) " + <> showT eliminateArtificialVarsFromPhase1Tableau + -- If the objecitve row is not found, the system is feasible iff + -- the artificial vars sum to zero. The value of an artificial + -- variable is 0 if non-basic, and the RHS of the row if basic + let artificialVarsVals = map (\v -> maybe 0 (.constant) (M.lookup v eliminateArtificialVarsFromPhase1Tableau)) artificialVars + let artificialVarsValsSum = sum artificialVarsVals + if artificialVarsValsSum == 0 + then do + logMsg LevelInfo $ + "findFeasibleSolution: Artifical variables sum up to 0, thus original tableau is feasible. System without artificial vars (in Dict form) " + <> showT eliminateArtificialVarsFromPhase1Tableau + pure . Just $ + FeasibleSystem + { dict = eliminateArtificialVarsFromPhase1Tableau + , slackVars = slackVars + , artificialVars = artificialVars + , objectiveVar = objectiveVar + } + else do + logMsg LevelInfo $ + "findFeasibleSolution: Artifical variables sum up to " + <> showT artificialVarsValsSum + <> ", thus original tableau is infeasible. System without artificial vars (in Dict form) " + <> showT eliminateArtificialVarsFromPhase1Tableau + pure Nothing + Just row -> + if row.constant == 0 + then do + logMsg LevelInfo $ + "findFeasibleSolution: Objective RHS is zero after pivoting, thus original tableau is feasible. feasible system (in Dict form) " + <> showT eliminateArtificialVarsFromPhase1Tableau + pure . Just $ + FeasibleSystem + { dict = eliminateArtificialVarsFromPhase1Tableau + , slackVars = slackVars + , artificialVars = artificialVars + , objectiveVar = objectiveVar + } + else do + unless (row.constant < 0) $ do + let errMsg = + "findFeasibleSolution: Objective RHS is negative after pivoting. This should be impossible. System without artificial vars (in Dict form) " + <> show eliminateArtificialVarsFromPhase1Tableau + logMsg LevelError $ Text.pack errMsg + error errMsg + logMsg LevelInfo $ + "findFeasibleSolution: Objective RHS not zero after phase 1, thus original tableau is infeasible. System without artificial vars (in Dict form) " + <> showT eliminateArtificialVarsFromPhase1Tableau + pure Nothing + Nothing -> do + logMsg LevelInfo $ + "findFeasibleSolution: Infeasible solution found, could not pivot with objective " + <> showT artificialPivotObjective + <> " over system (in Dict form) " + <> showT systemWithBasicVarsAsDictionary + pure Nothing + where + system = simplifySystem unsimplifiedSystem + + maxVar = + maximum $ + map + ( \case + LEQ vcm _ -> maximum (map fst $ M.toList vcm) + GEQ vcm _ -> maximum (map fst $ M.toList vcm) + EQ vcm _ -> maximum (map fst $ M.toList vcm) + ) + system + + (systemWithSlackVars, slackVars) = systemInStandardForm system maxVar [] + + maxVarWithSlackVars = if null slackVars then maxVar else maximum slackVars + + (systemWithBasicVars, artificialVars) = systemWithArtificialVars systemWithSlackVars maxVarWithSlackVars + + finalMaxVar = if null artificialVars then maxVarWithSlackVars else maximum artificialVars + + systemWithBasicVarsAsDictionary = tableauInDictionaryForm systemWithBasicVars + + artificialPivotObjective = createArtificialPivotObjective systemWithBasicVarsAsDictionary artificialVars + + objectiveVar = finalMaxVar + 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] -> Var -> [Var] -> ([(Maybe Var, PolyConstraint)], [Var]) + systemInStandardForm [] _ sVars = ([], sVars) + systemInStandardForm (EQ v r : xs) maxVar sVars = ((Nothing, EQ v r) : newSystem, newSlackVars) + where + (newSystem, newSlackVars) = systemInStandardForm xs maxVar sVars + systemInStandardForm (LEQ v r : xs) maxVar sVars = ((Just newSlackVar, EQ (M.insert newSlackVar 1 v) r) : newSystem, newSlackVars) + where + newSlackVar = maxVar + 1 + (newSystem, newSlackVars) = systemInStandardForm xs newSlackVar (newSlackVar : sVars) + systemInStandardForm (GEQ v r : xs) maxVar sVars = ((Just newSlackVar, EQ (M.insert newSlackVar (-1) v) r) : newSystem, newSlackVars) + where + newSlackVar = maxVar + 1 + (newSystem, newSlackVars) = systemInStandardForm xs newSlackVar (newSlackVar : 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 Var, PolyConstraint)] -> Var -> (Tableau, [Var]) + systemWithArtificialVars [] _ = (M.empty, []) + systemWithArtificialVars ((mVar, EQ v r) : pcs) maxVar = + case mVar of + Nothing -> + if r >= 0 + then + ( M.insert newArtificialVar (TableauRow {lhs = M.insert newArtificialVar 1 v, rhs = r}) newSystemWithNewMaxVar + , newArtificialVar : artificialVarsWithNewMaxVar + ) + else + ( M.insert newArtificialVar (TableauRow {lhs = M.insert newArtificialVar (-1) v, rhs = r}) newSystemWithNewMaxVar + , newArtificialVar : artificialVarsWithNewMaxVar + ) + Just basicVar -> + case M.lookup basicVar v of + Just basicVarCoeff -> + if r == 0 + then (M.insert basicVar (TableauRow {lhs = v, rhs = r}) newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) + else + if r > 0 + then + if basicVarCoeff >= 0 -- Should only be 1 in the standard call path + then (M.insert basicVar (TableauRow {lhs = v, rhs = r}) newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) + else + ( M.insert newArtificialVar (TableauRow {lhs = M.insert newArtificialVar 1 v, rhs = r}) newSystemWithNewMaxVar + , newArtificialVar : artificialVarsWithNewMaxVar -- Slack var is negative, r is positive (when original constraint was GEQ) + ) + else -- r < 0 + + if basicVarCoeff <= 0 -- Should only be -1 in the standard call path + then (M.insert basicVar (TableauRow {lhs = v, rhs = r}) newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) + else + ( M.insert newArtificialVar (TableauRow {lhs = M.insert newArtificialVar (-1) v, rhs = r}) newSystemWithNewMaxVar + , newArtificialVar : artificialVarsWithNewMaxVar -- Slack var is negative, r is negative (when original constraint was LEQ) + ) + Nothing -> error "1" -- undefined + where + newArtificialVar = maxVar + 1 + + (newSystemWithNewMaxVar, artificialVarsWithNewMaxVar) = systemWithArtificialVars pcs newArtificialVar + + (newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar) = systemWithArtificialVars pcs maxVar + systemWithArtificialVars _ _ = error "systemWithArtificialVars: given system includes non-EQ constraints" + + -- \| Takes a 'Dict' and a '[Var]' as input and returns a 'PivotObjective'. + -- The 'Dict' represents the tableau of a linear program with artificial + -- variables, and '[Var]' represents the artificial variables. + + -- The function first filters out the rows of the tableau that correspond + -- to the artificial variables, and negates them. It then computes the sum + -- of the negated rows, which represents the 'PivotObjective'. + createArtificialPivotObjective :: Dict -> [Var] -> PivotObjective + createArtificialPivotObjective rows artificialVars = + PivotObjective + { variable = objectiveVar + , function = foldVarLitMap $ map (.varMapSum) negatedRowsWithoutArtificialVars + , constant = sum $ map (.constant) negatedRowsWithoutArtificialVars + } + where + -- Filter out non-artificial entries + rowsToAdd = M.filterWithKey (\k _ -> k `elem` artificialVars) rows + negatedRows = M.map (\(DictValue rowVarMapSum rowConstant) -> DictValue (M.map negate rowVarMapSum) (negate rowConstant)) rowsToAdd + -- Negate rows, discard keys and artificial vars since the pivot objective does not care about them + negatedRowsWithoutArtificialVars = + map + ( \(_, DictValue {..}) -> + DictValue + { varMapSum = M.map negate $ M.filterWithKey (\k _ -> k `notElem` artificialVars) varMapSum + , constant = negate constant + } + ) + $ M.toList rowsToAdd + +-- | 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 :: (MonadIO m, MonadLogger m) => ObjectiveFunction -> FeasibleSystem -> m (Maybe Result) +optimizeFeasibleSystem objFunction fsys@(FeasibleSystem {dict = phase1Dict, ..}) = do + logMsg LevelInfo $ + "optimizeFeasibleSystem: Optimizing feasible system " <> showT fsys <> " with objective " <> showT objFunction + if null artificialVars + then do + logMsg LevelInfo $ + "optimizeFeasibleSystem: No artificial vars, system is feasible. Pivoting system (in dict form) " + <> showT phase1Dict + <> " with objective " + <> showT normalObjective + fmap (displayResults . dictionaryFormToTableau) <$> simplexPivot normalObjective phase1Dict + else do + logMsg LevelInfo $ + "optimizeFeasibleSystem: Artificial vars present. Pivoting system (in dict form) " + <> showT phase1Dict + <> " with objective " + <> showT adjustedObjective + fmap (displayResults . dictionaryFormToTableau) <$> simplexPivot adjustedObjective phase1Dict + where + -- \| displayResults takes a 'Tableau' and returns a 'Result'. The 'Tableau' + -- represents the final tableau of a linear program after the simplex + -- algorithm has been applied. The 'Result' contains the value of the + -- objective variable and a map of the values of all variables appearing + -- in the system, including the objective variable. + -- + -- The function first filters out the rows of the tableau that correspond + -- to the slack and artificial variables. It then extracts the values of + -- the remaining variables and stores them in a map. If the objective + -- function is a maximization problem, the map contains the values of the + -- variables as they appear in the final tableau. If the objective function + -- is a minimization problem, the map contains the values of the variables + -- as they appear in the final tableau, except for the objective variable, + -- which is negated. + displayResults :: Tableau -> Result + displayResults tableau = + Result + { objectiveVar = objectiveVar + , varValMap = extractVarVals + } + where + extractVarVals = + let tableauWithOriginalVars = + M.filterWithKey + ( \basicVarName _ -> + basicVarName `notElem` slackVars ++ artificialVars + ) + tableau + in case objFunction of + Max _ -> + M.map + ( \tableauRow -> + tableauRow.rhs + ) + tableauWithOriginalVars + Min _ -> + M.mapWithKey -- We maximized -objVar, so we negate the objVar to get the final value + ( \basicVarName tableauRow -> + if basicVarName == objectiveVar + then negate $ tableauRow.rhs + else tableauRow.rhs + ) + tableauWithOriginalVars + + -- \| Objective to use when optimising the linear program if no artificial + -- variables were necessary in the first phase. It is essentially the original + -- objective function, with a potential change of sign based on the type of + -- problem (Maximization or Minimization). + normalObjective :: PivotObjective + normalObjective = + PivotObjective + { variable = objectiveVar + , function = if isMax objFunction then objFunction.objective else M.map negate objFunction.objective + , constant = 0 + } + + -- \| Objective to use when optimising the linear program if artificial + -- variables were necessary in the first phase. It is an adjustment to the + -- original objective function, where the linear coefficients are modified + -- by back-substitution of the values of the artificial variables. + adjustedObjective :: PivotObjective + adjustedObjective = + PivotObjective + { variable = objectiveVar + , function = calcVarMap + , constant = calcConstants + } + where + -- \| Compute the adjustment to the constant term of the objective + -- function. It adds up the products of the original coefficients and + -- the corresponding constant term (rhs) of each artificial variable + -- in the phase 1 'Dict'. + calcConstants :: SimplexNum + calcConstants = + sum + $ map + ( \(var, coeff) -> + let multiplyWith = if isMax objFunction then coeff else -coeff + in case M.lookup var phase1Dict of + Nothing -> 0 + Just row -> row.constant * multiplyWith + ) + $ M.toList objFunction.objective + + -- \| Compute the adjustment to the coefficients of the original + -- variables in the objective function. It performs back-substitution + -- of the variables in the original objective function using the + -- current value of each artificial variable in the phase 1 'Dict'. + calcVarMap :: VarLitMapSum + calcVarMap = + foldVarLitMap $ + map + ( M.fromList + . ( \(var, coeff) -> + let multiplyWith = if isMax objFunction then coeff else -coeff + in case M.lookup var phase1Dict of + Nothing -> + [(var, multiplyWith)] + Just row -> map (second (* multiplyWith)) (M.toList $ row.varMapSum) + ) + ) + (M.toList objFunction.objective) + +-- | 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 :: (MonadIO m, MonadLogger m) => ObjectiveFunction -> [PolyConstraint] -> m (Maybe Result) +twoPhaseSimplex objFunction unsimplifiedSystem = do + logMsg LevelInfo $ + "twoPhaseSimplex: Solving system " <> showT unsimplifiedSystem <> " with objective " <> showT objFunction + phase1Result <- findFeasibleSolution unsimplifiedSystem + case phase1Result of + Just feasibleSystem -> do + logMsg LevelInfo $ + "twoPhaseSimplex: Feasible system found for " + <> showT unsimplifiedSystem + <> "; Feasible system: " + <> showT feasibleSystem + optimizedSystem <- optimizeFeasibleSystem objFunction feasibleSystem + logMsg LevelInfo $ + "twoPhaseSimplex: Optimized system found for " + <> showT unsimplifiedSystem + <> "; Optimized system: " + <> showT optimizedSystem + pure optimizedSystem + Nothing -> do + logMsg LevelInfo $ "twoPhaseSimplex: Phase 1 gives infeasible result for " <> showT unsimplifiedSystem + pure Nothing + +-- | Perform the simplex pivot algorithm on a system with basic vars, assume that the first row is the 'ObjectiveFunction'. +simplexPivot :: (MonadIO m, MonadLogger m) => PivotObjective -> Dict -> m (Maybe Dict) +simplexPivot objective@(PivotObjective {variable = objectiveVar, function = objectiveFunc, constant = objectiveConstant}) dictionary = do + logMsg LevelInfo $ + "simplexPivot: Pivoting with objective " <> showT objective <> " over system (in Dict form) " <> showT dictionary + case mostPositive objectiveFunc of + Nothing -> do + logMsg LevelInfo $ + "simplexPivot: Pivoting complete as no positive variables found in objective " + <> showT objective + <> " over system (in Dict form) " + <> showT dictionary + pure $ Just (insertPivotObjectiveToDict objective dictionary) + Just pivotNonBasicVar -> do + logMsg LevelInfo $ + "simplexPivot: Non-basic pivoting variable in objective, determined by largest coefficient = " <> showT pivotNonBasicVar + let mPivotBasicVar = ratioTest dictionary pivotNonBasicVar Nothing Nothing + case mPivotBasicVar of + Nothing -> do + logMsg LevelInfo $ + "simplexPivot: Ratio test failed with non-basic variable " + <> showT pivotNonBasicVar + <> " over system (in Dict form) " + <> showT dictionary + pure Nothing + Just pivotBasicVar -> do + logMsg LevelInfo $ "simplexPivot: Basic pivoting variable determined by ratio test " <> showT pivotBasicVar + logMsg LevelInfo $ + "simplexPivot: Pivoting with basic var " + <> showT pivotBasicVar + <> ", non-basic var " + <> showT pivotNonBasicVar + <> ", objective " + <> showT objective + <> " over system (in Dict form) " + <> showT dictionary + let pivotResult = pivot pivotBasicVar pivotNonBasicVar (insertPivotObjectiveToDict objective dictionary) + pivotedObj = + let pivotedObjEntry = fromMaybe (error "simplexPivot: Can't find objective after pivoting") $ M.lookup objectiveVar pivotResult + in objective & #function .~ pivotedObjEntry.varMapSum & #constant .~ pivotedObjEntry.constant + pivotedDict = M.delete objectiveVar pivotResult + logMsg LevelInfo $ + "simplexPivot: Pivoted, Recursing with new pivoting objective " + <> showT pivotedObj + <> " for new pivoted system (in Dict form) " + <> showT pivotedDict + simplexPivot + pivotedObj + pivotedDict + where + ratioTest :: Dict -> Var -> Maybe Var -> Maybe Rational -> Maybe Var + ratioTest dict = aux (M.toList dict) + where + aux :: [(Var, DictValue)] -> Var -> Maybe Var -> Maybe Rational -> Maybe Var + aux [] _ mCurrentMinBasicVar _ = mCurrentMinBasicVar + aux (x@(basicVar, dictEquation) : xs) mostNegativeVar mCurrentMinBasicVar mCurrentMin = + case M.lookup mostNegativeVar dictEquation.varMapSum of + Nothing -> aux xs mostNegativeVar mCurrentMinBasicVar mCurrentMin + Just currentCoeff -> + let dictEquationConstant = dictEquation.constant + in if currentCoeff >= 0 || dictEquationConstant < 0 + then aux xs mostNegativeVar mCurrentMinBasicVar mCurrentMin + else case mCurrentMin of + Nothing -> aux xs mostNegativeVar (Just basicVar) (Just (dictEquationConstant / currentCoeff)) + Just currentMin -> + if (dictEquationConstant / currentCoeff) >= currentMin + then aux xs mostNegativeVar (Just basicVar) (Just (dictEquationConstant / currentCoeff)) + else aux xs mostNegativeVar mCurrentMinBasicVar mCurrentMin + + mostPositive :: VarLitMapSum -> Maybe Var + mostPositive varLitMap = + case findLargestCoeff (M.toList varLitMap) Nothing of + Just (largestVarName, largestVarCoeff) -> + if largestVarCoeff <= 0 + then Nothing + else Just largestVarName + Nothing -> Nothing + where + findLargestCoeff :: [(Var, SimplexNum)] -> Maybe (Var, SimplexNum) -> Maybe (Var, SimplexNum) + findLargestCoeff [] mCurrentMax = mCurrentMax + findLargestCoeff (v@(vName, vCoeff) : vs) mCurrentMax = + case mCurrentMax of + Nothing -> findLargestCoeff vs (Just v) + Just (_, currentMaxCoeff) -> + if currentMaxCoeff >= vCoeff + then findLargestCoeff vs mCurrentMax + else findLargestCoeff vs (Just v) + + -- 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 :: Var -> Var -> Dict -> Dict + pivot leavingVariable enteringVariable dict = + case M.lookup enteringVariable (dictEntertingRow.varMapSum) of + Just enteringVariableCoeff -> + updatedRows + where + -- Move entering variable to basis, update other variables in row appropriately + pivotEnteringRow :: DictValue + pivotEnteringRow = + dictEntertingRow + & #varMapSum + %~ ( \basicEquation -> + -- uncurry + M.insert + leavingVariable + (-1) + (filterOutEnteringVarTerm basicEquation) + & traverse + %~ divideByNegatedEnteringVariableCoeff + ) + & #constant + %~ divideByNegatedEnteringVariableCoeff + where + newEnteringVarTerm = (leavingVariable, -1) + divideByNegatedEnteringVariableCoeff = (/ negate enteringVariableCoeff) + + -- Substitute pivot equation into other rows + updatedRows :: Dict + updatedRows = + M.fromList $ map (uncurry f2) $ M.toList dict + where + f entryVar entryVal = + if leavingVariable == entryVar + then pivotEnteringRow + else case M.lookup enteringVariable (entryVal.varMapSum) of + Just subsCoeff -> + entryVal + & #varMapSum + .~ combineVarLitMapSums + (pivotEnteringRow.varMapSum <&> (subsCoeff *)) + (filterOutEnteringVarTerm (entryVal.varMapSum)) + & #constant + .~ ((subsCoeff * (pivotEnteringRow.constant)) + entryVal.constant) + Nothing -> entryVal + + f2 :: Var -> DictValue -> (Var, DictValue) + f2 entryVar entryVal = + if leavingVariable == entryVar + then (enteringVariable, pivotEnteringRow) + else case M.lookup enteringVariable (entryVal.varMapSum) of + Just subsCoeff -> + ( entryVar + , entryVal + & #varMapSum + .~ combineVarLitMapSums + (pivotEnteringRow.varMapSum <&> (subsCoeff *)) + (filterOutEnteringVarTerm (entryVal.varMapSum)) + & #constant + .~ ((subsCoeff * (pivotEnteringRow.constant)) + entryVal.constant) + ) + Nothing -> (entryVar, entryVal) + Nothing -> error "pivot: non basic variable not found in basic row" + where + -- \| The entering row, i.e., the row in the dict which is the value of + -- leavingVariable. + dictEntertingRow = + fromMaybe + (error "pivot: Basic variable not found in Dict") + $ M.lookup leavingVariable dict + + filterOutEnteringVarTerm = M.filterWithKey (\vName _ -> vName /= enteringVariable) diff --git a/src/Linear/Simplex/Types.hs b/src/Linear/Simplex/Types.hs --- a/src/Linear/Simplex/Types.hs +++ b/src/Linear/Simplex/Types.hs @@ -1,46 +1,123 @@ -{-| -Module : Linear.Simplex.Types -Description : Custom types -Copyright : (c) Junaid Rasheed, 2020-2022 -License : BSD-3 -Maintainer : jrasheed178@gmail.com -Stability : experimental --} +-- | +-- Module : Linear.Simplex.Types +-- Description : Custom types +-- Copyright : (c) Junaid Rasheed, 2020-2023 +-- License : BSD-3 +-- Maintainer : jrasheed178@gmail.com +-- Stability : experimental module Linear.Simplex.Types where --- |List of 'Integer' variables with their 'Rational' coefficients. --- There is an implicit addition between elements in this list. --- Users must only provide positive integer variables. --- --- Example: [(2, 3), (6, (-1), (2, 1))] is equivalent to 3x2 + (-x6) + x2. -type VarConstMap = [(Integer, Rational)] +import Control.Lens +import Data.Generics.Labels () +import Data.List (sort) +import qualified Data.Map as M +import GHC.Generics (Generic) --- |For specifying constraints in a system. --- The LHS is a 'VarConstMap', and the RHS, is a 'Rational' number. --- LEQ [(1, 2), (2, 1)] 3.5 is equivalent to 2x1 + x2 <= 3.5. --- Users must only provide positive integer variables. --- --- Example: LEQ [(2, 3), (6, (-1), (2, 1))] 12.3 is equivalent to 3x2 + (-x6) + x2 <= 12.3. -data PolyConstraint = - LEQ VarConstMap Rational | - GEQ VarConstMap Rational | - EQ VarConstMap Rational deriving (Show, Eq); +type Var = Int --- |Create an objective function. --- We can either 'Max'imize or 'Min'imize a 'VarConstMap'. -data ObjectiveFunction = Max VarConstMap | Min VarConstMap deriving (Show, Eq) +type SimplexNum = Rational --- |A 'Tableau' of equations. --- Each pair in the list is a row. --- The first item in the pair specifies which 'Integer' variable is basic in the equation. --- The second item in the pair is an equation. --- The 'VarConstMap' in the second equation is a list of variables with their coefficients. --- The RHS of the equation is a 'Rational' constant. -type Tableau = [(Integer, (VarConstMap, Rational))] +type SystemRow = PolyConstraint --- |Type representing equations. --- Each pair in the list is one equation. --- The first item of the pair is the basic variable, and is on the LHS of the equation with a coefficient of one. --- The RHS is represented using a `VarConstMap`. --- The integer variable -1 is used to represent a 'Rational' on the RHS -type DictionaryForm = [(Integer, VarConstMap)] +type System = [SystemRow] + +-- A 'Tableau' where the basic variable may be empty. +-- All non-empty basic vars are slack vars +data SystemWithSlackVarRow = SystemInStandardFormRow + { mSlackVar :: Maybe Var + -- ^ This is Nothing iff the row does not have a slack variable + , row :: TableauRow + } + +type SystemWithSlackVars = [SystemWithSlackVarRow] + +data FeasibleSystem = FeasibleSystem + { dict :: Dict + , slackVars :: [Var] + , artificialVars :: [Var] + , objectiveVar :: Var + } + deriving (Show, Read, Eq, Generic) + +data Result = Result + { objectiveVar :: Var + , varValMap :: VarLitMap + -- TODO: + -- Maybe VarLitMap + -- , feasible :: Bool + -- , optimisable :: Bool + } + deriving (Show, Read, Eq, Generic) + +data SimplexMeta = SimplexMeta + { objective :: ObjectiveFunction + , feasibleSystem :: Maybe FeasibleSystem + , optimisedResult :: Maybe Result + } + +type VarLitMap = M.Map Var SimplexNum + +-- | List of variables with their 'SimplexNum' coefficients. +-- There is an implicit addition between elements in this list. +-- +-- Example: [Var "x" 3, Var "y" -1, Var "z" 1] is equivalent to 3x + (-y) + z. +type VarLitMapSum = VarLitMap + +-- | For specifying constraints in a system. +-- The LHS is a 'Vars', and the RHS, is a 'SimplexNum' number. +-- LEQ [(1, 2), (2, 1)] 3.5 is equivalent to 2x1 + x2 <= 3.5. +-- Users must only provide positive integer variables. +-- +-- Example: LEQ [Var "x" 3, Var "y" -1, Var "x" 1] 12.3 is equivalent to 3x + (-y) + x <= 12.3. +data PolyConstraint + = LEQ {lhs :: VarLitMapSum, rhs :: SimplexNum} + | GEQ {lhs :: VarLitMapSum, rhs :: SimplexNum} + | EQ {lhs :: VarLitMapSum, rhs :: SimplexNum} + deriving (Show, Read, Eq, Generic) + +-- | Create an objective function. +-- We can either 'Max'imize or 'Min'imize a 'VarTermSum'. +data ObjectiveFunction = Max {objective :: VarLitMapSum} | Min {objective :: VarLitMapSum} + deriving (Show, Read, Eq, Generic) + +-- | TODO: Maybe we want this type +-- TODO: A better/alternative name +data Equation = Equation + { lhs :: VarLitMapSum + , rhs :: SimplexNum + } + +-- | Value for 'Tableau'. lhs = rhs. +data TableauRow = TableauRow + { lhs :: VarLitMapSum + , rhs :: SimplexNum + } + deriving (Show, Read, Eq, Generic) + +-- | A simplex 'Tableu' of equations. +-- Each entry in the map is a row. +type Tableau = M.Map Var TableauRow + +-- | Values for a 'Dict'. +data DictValue = DictValue + { varMapSum :: VarLitMapSum + , constant :: SimplexNum + } + deriving (Show, Read, Eq, Generic) + +-- | A simplex 'Dict' +-- One quation represents the objective function. +-- Each pair in the list is one equation in the system we're working with. +-- data Dict = Dict +-- { objective :: DictObjective +-- , entries :: DictEntries +-- } +-- deriving (Show, Read, Eq, Generic) +type Dict = M.Map Var DictValue + +data PivotObjective = PivotObjective + { variable :: Var + , function :: VarLitMapSum + , constant :: SimplexNum + } + deriving (Show, Read, Eq, Generic) diff --git a/src/Linear/Simplex/Util.hs b/src/Linear/Simplex/Util.hs --- a/src/Linear/Simplex/Util.hs +++ b/src/Linear/Simplex/Util.hs @@ -1,153 +1,184 @@ -{-# LANGUAGE LambdaCase #-} - -{-| -Module : Linear.Simplex.Util -Description : Helper functions -Copyright : (c) Junaid Rasheed, 2020-2022 -License : BSD-3 -Maintainer : jrasheed178@gmail.com -Stability : experimental - -Helper functions for performing the two-phase simplex method. --} +-- | +-- Module : Linear.Simplex.Util +-- Description : Helper functions +-- Copyright : (c) Junaid Rasheed, 2020-2023 +-- License : BSD-3 +-- Maintainer : jrasheed178@gmail.com +-- Stability : experimental +-- +-- Helper functions for performing the two-phase simplex method. module Linear.Simplex.Util where -import Prelude hiding (EQ); -import Linear.Simplex.Types -import Data.List +import Control.Lens +import Control.Monad.IO.Class (MonadIO (..)) +import Control.Monad.Logger (LogLevel (..), LogLine, MonadLogger, logDebug, logError, logInfo, logWarn) import Data.Bifunctor +import Data.Generics.Labels () +import Data.Generics.Product (field) +import Data.List +import qualified Data.Map as Map +import qualified Data.Map.Merge.Lazy as MapMerge +import Data.Maybe (fromMaybe) +import qualified Data.Text as T +import Data.Time (getCurrentTime) +import Data.Time.Format.ISO8601 (iso8601Show) +import Linear.Simplex.Types +import Prelude hiding (EQ) --- |Is the given 'ObjectiveFunction' to be 'Max'imized? +-- | Is the given 'ObjectiveFunction' to be 'Max'imized? isMax :: ObjectiveFunction -> Bool isMax (Max _) = True isMax (Min _) = False --- |Extract the objective ('VarConstMap') from an 'ObjectiveFunction' -getObjective :: ObjectiveFunction -> VarConstMap -getObjective (Max o) = o -getObjective (Min o) = o - --- |Simplifies a system of 'PolyConstraint's by first calling 'simplifyPolyConstraint', --- then reducing 'LEQ' and 'GEQ' with same LHS and RHS (and other similar situations) into 'EQ', --- and finally removing duplicate elements using 'nub'. +-- | Simplifies a system of 'PolyConstraint's by first calling 'simplifyPolyConstraint', +-- then reducing 'LEQ' and 'GEQ' with same LHS and RHS (and other similar situations) into 'EQ', +-- and finally removing duplicate elements using 'nub'. simplifySystem :: [PolyConstraint] -> [PolyConstraint] -simplifySystem = nub . reduceSystem . map simplifyPolyConstraint +simplifySystem = nub . reduceSystem where reduceSystem :: [PolyConstraint] -> [PolyConstraint] reduceSystem [] = [] -- Reduce LEQ with matching GEQ and EQ into EQ reduceSystem ((LEQ lhs rhs) : pcs) = - let - matchingConstraints = - filter - (\case - GEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - EQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - _ -> False - ) - pcs - in - if null matchingConstraints - then LEQ lhs rhs : reduceSystem pcs - else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) + let matchingConstraints = + filter + ( \case + GEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' + EQ lhs' rhs' -> lhs == lhs' && rhs == rhs' + _ -> False + ) + pcs + in if null matchingConstraints + then LEQ lhs rhs : reduceSystem pcs + else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) -- Reduce GEQ with matching LEQ and EQ into EQ reduceSystem ((GEQ lhs rhs) : pcs) = - let - matchingConstraints = - filter - (\case - LEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - EQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - _ -> False - ) - pcs - in - if null matchingConstraints - then GEQ lhs rhs : reduceSystem pcs - else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) + let matchingConstraints = + filter + ( \case + LEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' + EQ lhs' rhs' -> lhs == lhs' && rhs == rhs' + _ -> False + ) + pcs + in if null matchingConstraints + then GEQ lhs rhs : reduceSystem pcs + else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) -- Reduce EQ with matching LEQ and GEQ into EQ reduceSystem ((EQ lhs rhs) : pcs) = - let - matchingConstraints = - filter - (\case - LEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - GEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' - _ -> False - ) - pcs - in - if null matchingConstraints - then EQ lhs rhs : reduceSystem pcs - else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) + let matchingConstraints = + filter + ( \case + LEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' + GEQ lhs' rhs' -> lhs == lhs' && rhs == rhs' + _ -> False + ) + pcs + in if null matchingConstraints + then EQ lhs rhs : reduceSystem pcs + else EQ lhs rhs : reduceSystem (pcs \\ matchingConstraints) --- |Simplify an 'ObjectiveFunction' by first 'sort'ing and then calling 'foldSumVarConstMap' on the 'VarConstMap'. -simplifyObjectiveFunction :: ObjectiveFunction -> ObjectiveFunction -simplifyObjectiveFunction (Max varConstMap) = Max (foldSumVarConstMap (sort varConstMap)) -simplifyObjectiveFunction (Min varConstMap) = Min (foldSumVarConstMap (sort varConstMap)) +-- | Converts a 'Dict' to a 'Tableau' using 'dictEntryToTableauEntry'. +-- FIXME: maybe remove this line. The basic variables will have a coefficient of 1 in the 'Tableau'. +dictionaryFormToTableau :: Dict -> Tableau +dictionaryFormToTableau = + Map.mapWithKey + ( \basicVar (DictValue {..}) -> + TableauRow + { lhs = Map.insert basicVar 1 $ negate <$> varMapSum + , rhs = constant + } + ) --- |Simplify a 'PolyConstraint' by first 'sort'ing and then calling 'foldSumVarConstMap' on the 'VarConstMap'. -simplifyPolyConstraint :: PolyConstraint -> PolyConstraint -simplifyPolyConstraint (LEQ varConstMap rhs) = LEQ (foldSumVarConstMap (sort varConstMap)) rhs -simplifyPolyConstraint (GEQ varConstMap rhs) = GEQ (foldSumVarConstMap (sort varConstMap)) rhs -simplifyPolyConstraint (EQ varConstMap rhs) = EQ (foldSumVarConstMap (sort varConstMap)) rhs +-- | Converts a 'Tableau' to a 'Dict'. +-- We do this by isolating the basic variable on the LHS, ending up with all non basic variables and a 'SimplexNum' constant on the RHS. +tableauInDictionaryForm :: Tableau -> Dict +tableauInDictionaryForm = + Map.mapWithKey + ( \basicVar (TableauRow {..}) -> + let basicVarCoeff = fromMaybe 1 $ Map.lookup basicVar lhs + in DictValue + { varMapSum = + Map.map + (\c -> negate c / basicVarCoeff) + $ Map.delete basicVar lhs + , constant = rhs / basicVarCoeff + } + ) --- |Add a sorted list of 'VarConstMap's, folding where the variables are equal -foldSumVarConstMap :: [(Integer, Rational)] -> [(Integer, Rational)] -foldSumVarConstMap [] = [] -foldSumVarConstMap [(v, c)] = [(v, c)] -foldSumVarConstMap ((v1, c1) : (v2, c2) : vcm) = - if v1 == v2 - then - let newC = c1 + c2 - in - if newC == 0 - then foldSumVarConstMap vcm - else foldSumVarConstMap $ (v1, c1 + c2) : vcm - else (v1, c1) : foldSumVarConstMap ((v2, c2) : vcm) +-- | If this function is given 'Nothing', return 'Nothing'. +-- Otherwise, we 'lookup' the 'Integer' given in the first item of the pair in the map given in the second item of the pair. +-- This is typically used to extract the value of the 'ObjectiveFunction' after calling 'Linear.Simplex.Solver.TwoPhase.twoPhaseSimplex'. +extractObjectiveValue :: Maybe Result -> Maybe SimplexNum +extractObjectiveValue = fmap $ \result -> + case Map.lookup result.objectiveVar result.varValMap of + Nothing -> error "Objective not found in results when extracting objective value" + Just r -> r --- |Get a map of the value of every 'Integer' variable in a 'Tableau' -displayTableauResults :: Tableau -> [(Integer, Rational)] -displayTableauResults = map (\(basicVar, (_, rhs)) -> (basicVar, rhs)) +-- | Combines two 'VarLitMapSums together by summing values with matching keys +combineVarLitMapSums :: VarLitMapSum -> VarLitMapSum -> VarLitMapSum +combineVarLitMapSums = + MapMerge.merge + (MapMerge.mapMaybeMissing keepVal) + (MapMerge.mapMaybeMissing keepVal) + (MapMerge.zipWithMaybeMatched sumVals) + where + keepVal = const pure + sumVals k v1 v2 = Just $ v1 + v2 --- |Get a map of the value of every 'Integer' variable in a 'DictionaryForm' -displayDictionaryResults :: DictionaryForm -> [(Integer, Rational)] -displayDictionaryResults dict = displayTableauResults$ dictionaryFormToTableau dict +foldDictValue :: [DictValue] -> DictValue +foldDictValue [] = error "Empty list of DictValues given to foldDictValue" +foldDictValue [x] = x +foldDictValue (DictValue {varMapSum = vm1, constant = c1} : DictValue {varMapSum = vm2, constant = c2} : dvs) = + let combinedDictValue = + DictValue + { varMapSum = foldVarLitMap [vm1, vm2] + , constant = c1 + c2 + } + in foldDictValue $ combinedDictValue : dvs --- |Map the given 'Integer' variable to the given 'ObjectiveFunction', for entering into 'DictionaryForm'. -createObjectiveDict :: ObjectiveFunction -> Integer -> (Integer, VarConstMap) -createObjectiveDict (Max obj) objectiveVar = (objectiveVar, obj) -createObjectiveDict (Min obj) objectiveVar = (objectiveVar, map (second negate) obj) +foldVarLitMap :: [VarLitMap] -> VarLitMap +foldVarLitMap [] = error "Empty list of VarLitMaps given to foldVarLitMap" +foldVarLitMap [x] = x +foldVarLitMap (vm1 : vm2 : vms) = + let combinedVars = nub $ Map.keys vm1 <> Map.keys vm2 --- |Converts a 'Tableau' to 'DictionaryForm'. --- We do this by isolating the basic variable on the LHS, ending up with all non basic variables and a 'Rational' constant on the RHS. --- (-1) is used to represent the rational constant. -tableauInDictionaryForm :: Tableau -> DictionaryForm -tableauInDictionaryForm [] = [] -tableauInDictionaryForm ((basicVar, (vcm, r)) : rows) = - (basicVar, (-1, r / basicCoeff) : map (\(v, c) -> (v, negate c / basicCoeff)) nonBasicVars) : tableauInDictionaryForm rows - where - basicCoeff = if null basicVars then 1 else snd $ head basicVars - (basicVars, nonBasicVars) = partition (\(v, _) -> v == basicVar) vcm + combinedVarMap = + Map.fromList $ + map + ( \var -> + let mVm1VarVal = Map.lookup var vm1 + mVm2VarVal = Map.lookup var vm2 + in ( var + , case (mVm1VarVal, mVm2VarVal) of + (Just vm1VarVal, Just vm2VarVal) -> vm1VarVal + vm2VarVal + (Just vm1VarVal, Nothing) -> vm1VarVal + (Nothing, Just vm2VarVal) -> vm2VarVal + (Nothing, Nothing) -> error "Reached unreachable branch in foldDictValue" + ) + ) + combinedVars + in foldVarLitMap $ combinedVarMap : vms --- |Converts a 'DictionaryForm' to a 'Tableau'. --- This is done by moving all non-basic variables from the right to the left. --- The rational constant (represented by the 'Integer' variable -1) stays on the right. --- The basic variables will have a coefficient of 1 in the 'Tableau'. -dictionaryFormToTableau :: DictionaryForm -> Tableau -dictionaryFormToTableau [] = [] -dictionaryFormToTableau ((basicVar, row) : rows) = - (basicVar, ((basicVar, 1) : map (second negate) nonBasicVars, r)) : dictionaryFormToTableau rows - where - (rationalConstant, nonBasicVars) = partition (\(v,_) -> v == (-1)) row - r = if null rationalConstant then 0 else (snd . head) rationalConstant -- If there is no rational constant found in the right side, the rational constant is 0. +insertPivotObjectiveToDict :: PivotObjective -> Dict -> Dict +insertPivotObjectiveToDict objective = Map.insert objective.variable (DictValue {varMapSum = objective.function, constant = objective.constant}) --- |If this function is given 'Nothing', return 'Nothing'. --- Otherwise, we 'lookup' the 'Integer' given in the first item of the pair in the map given in the second item of the pair. --- This is typically used to extract the value of the 'ObjectiveFunction' after calling 'Linear.Simplex.Simplex.twoPhaseSimplex'. -extractObjectiveValue :: Maybe (Integer, [(Integer, Rational)]) -> Maybe Rational -extractObjectiveValue Nothing = Nothing -extractObjectiveValue (Just (objVar, results)) = - case lookup objVar results of - Nothing -> error "Objective not found in results when extracting objective value" - r -> r +showT :: (Show a) => a -> T.Text +showT = T.pack . show + +logMsg :: (MonadIO m, MonadLogger m) => LogLevel -> T.Text -> m () +logMsg lvl msg = do + currTime <- T.pack . iso8601Show <$> liftIO getCurrentTime + let msgToLog = currTime <> ": " <> msg + case lvl of + LevelDebug -> $logDebug msgToLog + LevelInfo -> $logInfo msgToLog + LevelWarn -> $logWarn msgToLog + LevelError -> $logError msgToLog + LevelOther otherLvl -> error "logMsg: LevelOther is not implemented" + +extractTableauValues :: Tableau -> Map.Map Var SimplexNum +extractTableauValues = Map.map (.rhs) + +extractDictValues :: Dict -> Map.Map Var SimplexNum +extractDictValues = Map.map (.constant) diff --git a/test/Spec.hs b/test/Spec.hs --- a/test/Spec.hs +++ b/test/Spec.hs @@ -1,28 +1,42 @@ module Main where -import Linear.Simplex.Simplex +import Control.Monad +import Control.Monad.IO.Class +import Control.Monad.Logger + import Linear.Simplex.Prettify +import Linear.Simplex.Solver.TwoPhase +import Linear.Simplex.Types import Linear.Simplex.Util + import TestFunctions main :: IO () -main = runTests testsList +main = runStdoutLoggingT $ filterLogger (\_logSource logLevel -> logLevel > LevelInfo) $ runTests testsList -runTests [] = putStrLn "All tests passed" +runTests :: (MonadLogger m, MonadFail m, MonadIO m) => [((ObjectiveFunction, [PolyConstraint]), Maybe Result)] -> m () +runTests [] = do + liftIO $ putStrLn "All tests passed" + pure () runTests (((testObjective, testConstraints), expectedResult) : tests) = - let testResult = twoPhaseSimplex testObjective testConstraints in - if testResult == expectedResult - then runTests tests - else do - putStrLn "The following test failed: \n" - putStrLn ("Objective Function (Non-prettified): " ++ show testObjective) - putStrLn ("Constraints (Non-prettified): " ++ show testConstraints) - putStrLn "====================================\n" - putStrLn ("Objective Function (Prettified): " ++ prettyShowObjectiveFunction testObjective) - putStrLn "Constraints (Prettified): " - putStrLn (concatMap ((\c -> "\t" ++ prettyShowPolyConstraint c ++ "\n")) testConstraints) - putStrLn "====================================\n" - putStrLn ("Expected Solution (Full): " ++ show expectedResult) - putStrLn ("Actual Solution (Full): " ++ show testResult) - putStrLn ("Expected Solution (Objective): " ++ show (extractObjectiveValue expectedResult)) - putStrLn ("Actual Solution (Objective): " ++ show (extractObjectiveValue testResult)) + do + testResult <- twoPhaseSimplex testObjective testConstraints + if testResult == expectedResult + then runTests tests + else do + let msg = + "\nThe following test failed: " + <> ("\nObjective Function (Non-prettified): " ++ show testObjective) + <> ("\nConstraints (Non-prettified): " ++ show testConstraints) + <> "\n====================================" + <> ("\nObjective Function (Prettified): " ++ prettyShowObjectiveFunction testObjective) + <> "\nConstraints (Prettified): " + <> "\n" + <> concatMap (\c -> "\t" ++ prettyShowPolyConstraint c ++ "\n") testConstraints + <> "\n====================================" + <> ("\nExpected Solution (Full): " ++ show expectedResult) + <> ("\nActual Solution (Full): " ++ show testResult) + <> ("\nExpected Solution (Objective): " ++ show (extractObjectiveValue expectedResult)) + <> ("\nActual Solution (Objective): " ++ show (extractObjectiveValue testResult)) + <> "\n" + fail msg diff --git a/test/TestFunctions.hs b/test/TestFunctions.hs --- a/test/TestFunctions.hs +++ b/test/TestFunctions.hs @@ -1,1078 +1,1048 @@ module TestFunctions where -import Prelude hiding (EQ) -import Linear.Simplex.Types -import Data.Ratio - -testsList :: [((ObjectiveFunction, [PolyConstraint]), Maybe (Integer, [(Integer, Rational)]))] -testsList = - [ - (test1, Just (7,[(7,29 % 1),(1,3 % 1),(2,4 % 1)])) - , (test2, Just (7,[(7,0 % 1)])) - , (test3, Nothing) - , (test4, Just (11,[(11,237 % 7),(1,24 % 7),(2,33 % 7)])) - , (test5, Just (9,[(9,3 % 5),(2,14 % 5),(3,17 % 5)])) - , (test6, Nothing) - , (test7, Just (8,[(8,1 % 1),(2,2 % 1),(1,3 % 1)])) - , (test8, Just (8,[(8,(-1) % 4),(2,9 % 2),(1,17 % 4)])) - , (test9, Just (7,[(7,5 % 1),(3,2 % 1),(4,1 % 1)])) - , (test10, Just (7,[(7,8 % 1),(1,2 % 1),(2,6 % 1)])) - , (test11, Just (8,[(8,20 % 1),(4,16 % 1),(3,6 % 1)])) - , (test12, Just (8,[(8,6 % 1),(4,2 % 1),(5,2 % 1)])) - , (test13, Just (6,[(6,150 % 1),(2,150 % 1)])) - , (test14, Just (6,[(6,40 % 3),(2,40 % 3)])) - , (test15, Nothing) - , (test16, Just (6,[(6,75 % 1),(1,75 % 2)])) - , (test17, Just (7,[(7,(-120) % 1),(1,20 % 1)])) - , (test18, Just (7,[(7,10 % 1),(3,5 % 1)])) - , (test19, Nothing) - , (test20, Nothing) - , (test21, Just (7,[(7,250 % 1),(2,50 % 1)])) - , (test22, Just (7,[(7,0 % 1)])) - , (test23, Nothing) - , (test24, Just (10,[(10,300 % 1),(3,150 % 1)])) - , (test25, Just (3,[(3,15 % 1),(1,15 % 1)])) - , (test26, Just (6,[(6,20 % 1),(1,10 % 1),(2,10 % 1)])) - , (test27, Just (3,[(3,0 % 1)])) - , (test28, Just (6,[(6,0 % 1),(2,10 % 1)])) - , (test29, Nothing) - , (test30, Nothing) - , (testPolyPaver1, Just (12,[(12,7 % 4),(2,5 % 2),(1,7 % 4),(3,0 % 1)])) - , (testPolyPaver2, Just (12,[(12,5 % 2),(2,5 % 3),(1,5 % 2),(3,0 % 1)])) - , (testPolyPaver3, Just (12,[(12,5 % 3),(2,5 % 3),(1,5 % 2),(3,0 % 1)])) - , (testPolyPaver4, Just (12,[(12,5 % 2),(2,5 % 2),(1,5 % 2),(3,0 % 1)])) - , (testPolyPaver5, Nothing) - , (testPolyPaver6, Nothing) - , (testPolyPaver7, Nothing) - , (testPolyPaver8, Nothing) - , (testPolyPaver9, Just (12,[(12,7 % 2),(2,5 % 9),(1,7 % 2),(3,0 % 1)])) - , (testPolyPaver10, Just (12,[(12,17 % 20),(2,7 % 2),(1,17 % 20),(3,0 % 1)])) - , (testPolyPaver11, Just (12,[(12,7 % 2),(2,7 % 2),(1,22 % 9)])) - , (testPolyPaver12, Just (12,[(12,5 % 9),(2,5 % 9),(1,7 % 2),(3,0 % 1)])) - , (testPolyPaverTwoFs1, Nothing) - , (testPolyPaverTwoFs2, Nothing) - , (testPolyPaverTwoFs3, Nothing) - , (testPolyPaverTwoFs4, Nothing) - , (testPolyPaverTwoFs5, Just (17,[(17,5 % 2),(2,45 % 22),(1,5 % 2),(4,0 % 1)])) - , (testPolyPaverTwoFs6, Just (17,[(17,45 % 22),(2,5 % 2),(1,45 % 22),(4,0 % 1)])) - , (testPolyPaverTwoFs7, Just (17,[(17,5 % 2),(2,5 % 2),(1,5 % 2),(4,0 % 1)])) - , (testPolyPaverTwoFs8, Just (17,[(17,45 % 22),(2,45 % 22),(1,5 % 2),(4,0 % 1)])) - , (testLeqGeqBugMin1, Just (5,[(5,3 % 1),(1,3 % 1),(2,3 % 1)])) - , (testLeqGeqBugMax1, Just (5,[(5,3 % 1),(1,3 % 1),(2,3 % 1)])) - , (testLeqGeqBugMin2, Just (5,[(5,3 % 1),(1,3 % 1),(2,3 % 1)])) - , (testLeqGeqBugMax2, Just (5,[(5,3 % 1),(1,3 % 1),(2,3 % 1)])) - , (testQuickCheck1, Just (10,[(10,(-370) % 1),(2,26 % 1),(1,5 % 3)])) - , (testQuickCheck2, Just (8,[(8,(-2) % 9),(1,14 % 9),(2,8 % 9)])) - , (testQuickCheck3, Just (7,[(7,(-8) % 1),(2,2 % 1)])) - ] - -testLeqGeqBugMin1 = - ( - Min [(1, 1)], - [ - GEQ [(1,1 % 1)] (3 % 1), - LEQ [(1,1 % 1)] (3 % 1), - GEQ [(2,1 % 1)] (3 % 1), - LEQ [(2,1 % 1)] (3 % 1) - ] - ) - -testLeqGeqBugMax1 = - ( - Min [(1, 1)], - [ - GEQ [(1,1 % 1)] (3 % 1), - LEQ [(1,1 % 1)] (3 % 1), - GEQ [(2,1 % 1)] (3 % 1), - LEQ [(2,1 % 1)] (3 % 1) - ] - ) - -testLeqGeqBugMin2 = - ( - Min [(1, 1)], - [ - GEQ [(1,1 % 1)] (3 % 1), - LEQ [(1,1 % 1)] (3 % 1), - GEQ [(2,1 % 1)] (3 % 1), - LEQ [(2,1 % 1)] (3 % 1) - ] - ) - -testLeqGeqBugMax2 = - ( - Min [(1, 1)], - [ - GEQ [(1,1 % 1)] (3 % 1), - LEQ [(1,1 % 1)] (3 % 1), - GEQ [(2,1 % 1)] (3 % 1), - LEQ [(2,1 % 1)] (3 % 1) - ] - ) - --- From page 50 of 'Linear and Integer Programming Made Easy' --- Solution: obj = 29, 1 = 3, 2 = 4, -test1 :: (ObjectiveFunction, [PolyConstraint]) -test1 = - ( - Max [(1, 3), (2, 5)], - [ - LEQ [(1, 3), (2, 1)] 15, - LEQ [(1, 1), (2, 1)] 7, - LEQ [(2, 1)] 4, - LEQ [(1, -1), (2, 2)] 6 - ] - ) - -test2 :: (ObjectiveFunction, [PolyConstraint]) -test2 = - ( - Min [(1, 3), (2, 5)], - [ - LEQ [(1, 3), (2, 1)] 15, - LEQ [(1, 1), (2, 1)] 7, - LEQ [(2, 1)] 4, - LEQ [(1, -1), (2, 2)] 6 - ] - ) - -test3 :: (ObjectiveFunction, [PolyConstraint]) -test3 = - ( - Max [(1, 3), (2, 5)], - [ - GEQ [(1, 3), (2, 1)] 15, - GEQ [(1, 1), (2, 1)] 7, - GEQ [(2, 1)] 4, - GEQ [(1, -1), (2, 2)] 6 - ] - ) - -test4 :: (ObjectiveFunction, [PolyConstraint]) -test4 = - ( - Min [(1, 3), (2, 5)], - [ - GEQ [(1, 3), (2, 1)] 15, - GEQ [(1, 1), (2, 1)] 7, - GEQ [(2, 1)] 4, - GEQ [(1, -1), (2, 2)] 6 - ] - ) - --- From https://www.eng.uwaterloo.ca/~syde05/phase1.pdf --- Solution: obj = 3/5, 2 = 14/5, 3 = 17/5 --- requires two phases -test5 :: (ObjectiveFunction, [PolyConstraint]) -test5 = - ( - Max [(1, 1), (2, -1), (3, 1)], - [ - LEQ [(1, 2), (2, -1), (3, 2)] 4, - LEQ [(1, 2), (2, -3), (3, 1)] (-5), - LEQ [(1, -1), (2, 1), (3, -2)] (-1) - ] - ) - -test6 :: (ObjectiveFunction, [PolyConstraint]) -test6 = - ( - Min [(1, 1), (2, -1), (3, 1)], - [ - LEQ [(1, 2), (2, -1), (3, 2)] 4, - LEQ [(1, 2), (2, -3), (3, 1)] (-5), - LEQ [(1, -1), (2, 1), (3, -2)] (-1) - ] - ) -test7 :: (ObjectiveFunction, [PolyConstraint]) -test7 = - ( - Max [(1, 1), (2, -1), (3, 1)], - [ - GEQ [(1, 2), (2, -1), (3, 2)] 4, - GEQ [(1, 2), (2, -3), (3, 1)] (-5), - GEQ [(1, -1), (2, 1), (3, -2)] (-1) - ] - ) -test8 :: (ObjectiveFunction, [PolyConstraint]) -test8 = - ( - Min [(1, 1), (2, -1), (3, 1)], - [ - GEQ [(1, 2), (2, -1), (3, 2)] 4, - GEQ [(1, 2), (2, -3), (3, 1)] (-5), - GEQ [(1, -1), (2, 1), (3, -2)] (-1) - ] - ) - --- From page 49 of 'Linear and Integer Programming Made Easy' --- Solution: obj = -5, 3 = 2, 4 = 1, objVar was negated so actual val is 5 wa --- requires two phases -test9 :: (ObjectiveFunction, [PolyConstraint]) -test9 = - ( - Min [(1, 1), (2, 1), (3, 2), (4, 1)], - [ - EQ [(1, 1), (3, 2), (4, -2)] 2, - EQ [(2, 1), (3, 1), (4, 4)] 6 - ] - ) - -test10 :: (ObjectiveFunction, [PolyConstraint]) -test10 = - ( - Max [(1, 1), (2, 1), (3, 2), (4, 1)], - [ - EQ [(1, 1), (3, 2), (4, -2)] 2, - EQ [(2, 1), (3, 1), (4, 4)] 6 - ] - ) - --- Adapted from page 52 of 'Linear and Integer Programming Made Easy' --- Removed variables which do not appear in the system (these should be artificial variables) --- Solution: obj = 20, 3 = 6, 4 = 16 wq -test11 :: (ObjectiveFunction, [PolyConstraint]) -test11 = - ( - Max [(3, -2), (4, 2), (5, 1)], - [ - EQ [(3, -2), (4, 1), (5, 1)] 4, - EQ [(3, 3), (4, -1), (5, 2)] 2 - ] - ) - -test12 :: (ObjectiveFunction, [PolyConstraint]) -test12 = - ( - Min [(3, -2), (4, 2), (5, 1)], - [ - EQ [(3, -2), (4, 1), (5, 1)] 4, - EQ [(3, 3), (4, -1), (5, 2)] 2 - ] - ) - --- From page 59 of 'Linear and Integer Programming Made Easy' --- Solution: obj = 150, 1 = 0, 2 = 150 --- requires two phases -test13 :: (ObjectiveFunction, [PolyConstraint]) -test13 = - ( - Max [(1, 2), (2, 1)], - [ - LEQ [(1, 4), (2, 1)] 150, - LEQ [(1, 2), (2, -3)] (-40) - ] - ) - -test14 :: (ObjectiveFunction, [PolyConstraint]) -test14 = - ( - Min [(1, 2), (2, 1)], - [ - LEQ [(1, 4), (2, 1)] 150, - LEQ [(1, 2), (2, -3)] (-40) - ] - ) - -test15 :: (ObjectiveFunction, [PolyConstraint]) -test15 = - ( - Max [(1, 2), (2, 1)], - [ - GEQ [(1, 4), (2, 1)] 150, - GEQ [(1, 2), (2, -3)] (-40) - ] - ) - -test16 :: (ObjectiveFunction, [PolyConstraint]) -test16 = - ( - Min [(1, 2), (2, 1)], - [ - GEQ [(1, 4), (2, 1)] 150, - GEQ [(1, 2), (2, -3)] (-40) - ] - ) - --- From page 59 of 'Linear and Integer Programming Made Easy' --- Solution: obj = 120, 1 = 20, 2 = 0, 3 = 0, objVar was negated so actual val is -120 -test17 :: (ObjectiveFunction, [PolyConstraint]) -test17 = - ( - Min [(1, -6), (2, -4), (3, 2)], - [ - LEQ [(1, 1), (2, 1), (3, 4)] 20, - LEQ [(2, -5), (3, 5)] 100, - LEQ [(1, 1), (3, 1), (1, 1)] 400 - ] - ) - -test18 :: (ObjectiveFunction, [PolyConstraint]) -test18 = - ( - Max [(1, -6), (2, -4), (3, 2)], - [ - LEQ [(1, 1), (2, 1), (3, 4)] 20, - LEQ [(2, -5), (3, 5)] 100, - LEQ [(1, 1), (3, 1), (1, 1)] 400 - ] - ) - -test19 :: (ObjectiveFunction, [PolyConstraint]) -test19 = - ( - Min [(1, -6), (2, -4), (3, 2)], - [ - GEQ [(1, 1), (2, 1), (3, 4)] 20, - GEQ [(2, -5), (3, 5)] 100, - GEQ [(1, 1), (3, 1), (1, 1)] 400 - ] - ) - -test20 :: (ObjectiveFunction, [PolyConstraint]) -test20 = - ( - Max [(1, -6), (2, -4), (3, 2)], - [ - GEQ [(1, 1), (2, 1), (3, 4)] 20, - GEQ [(2, -5), (3, 5)] 100, - GEQ [(1, 1), (3, 1), (1, 1)] 400 - ] - ) - --- From page 59 of 'Linear and Integer Programming Made Easy' --- Solution: obj = 250, 1 = 0, 2 = 50, 3 = 0 -test21 :: (ObjectiveFunction, [PolyConstraint]) -test21 = - ( - Max [(1, 3), (2, 5), (3, 2)], - [ - LEQ [(1, 5), (2, 1), (3, 4)] 50, - LEQ [(1, 1), (2, -1), (3, 1)] 150, - LEQ [(1, 2), (2, 1), (3, 2)] 100 - ] - ) - -test22 :: (ObjectiveFunction, [PolyConstraint]) -test22 = - ( - Min [(1, 3), (2, 5), (3, 2)], - [ - LEQ [(1, 5), (2, 1), (3, 4)] 50, - LEQ [(1, 1), (2, -1), (3, 1)] 150, - LEQ [(1, 2), (2, 1), (3, 2)] 100 - ] - ) - -test23 :: (ObjectiveFunction, [PolyConstraint]) -test23 = - ( - Max [(1, 3), (2, 5), (3, 2)], - [ - GEQ [(1, 5), (2, 1), (3, 4)] 50, - GEQ [(1, 1), (2, -1), (3, 1)] 150, - GEQ [(1, 2), (2, 1), (3, 2)] 100 - ] - ) - -test24 :: (ObjectiveFunction, [PolyConstraint]) -test24 = - ( - Min [(1, 3), (2, 5), (3, 2)], - [ - GEQ [(1, 5), (2, 1), (3, 4)] 50, - GEQ [(1, 1), (2, -1), (3, 1)] 150, - GEQ [(1, 2), (2, 1), (3, 2)] 100 - ] - ) - -test25 :: (ObjectiveFunction, [PolyConstraint]) -test25 = - ( - Max [(1, 1)], - [ - LEQ [(1, 1)] 15 - ] - ) - -test26 :: (ObjectiveFunction, [PolyConstraint]) -test26 = - ( - Max [(1, 2)], - [ - LEQ [(1, 2)] 20, - GEQ [(2, 1)] 10 - ] - ) - -test27 :: (ObjectiveFunction, [PolyConstraint]) -test27 = - ( - Min [(1, 1)], - [ - LEQ [(1, 1)] 15 - ] - ) - -test28 :: (ObjectiveFunction, [PolyConstraint]) -test28 = - ( - Min [(1, 2)], - [ - LEQ [(1, 2)] 20, - GEQ [(2, 1)] 10 - ] - ) - -test29 :: (ObjectiveFunction, [PolyConstraint]) -test29 = - ( - Max [(1, 1)], - [ - LEQ [(1, 1)] 15, - GEQ [(1, 1)] 15.01 - ] - ) - -test30 :: (ObjectiveFunction, [PolyConstraint]) -test30 = - ( - Max [(1, 1)], - [ - LEQ [(1, 1)] 15, - GEQ [(1, 1)] 15.01, - GEQ [(2, 1)] 10 - ] - ) - --- Tests for systems similar to those from PolyPaver2 -testPolyPaver1 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver1 = - ( - Min [(1 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver2 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver2 = - ( - Max [(1 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver3 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver3 = - ( - Min [(2 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver4 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver4 = - ( - Max [(2 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver5 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver5 = - ( - Max [(1 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 1.5 - x2l = 0.0 - x2r = 1.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver6 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver6 = - ( - Min [(1 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 1.5 - x2l = 0.0 - x2r = 1.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver7 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver7 = - ( - Max [(2 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 1.5 - x2l = 0.0 - x2r = 1.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver8 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver8 = - ( - Min [(2 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 1.5 - x2l = 0.0 - x2r = 1.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver9 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver9 = - ( - Max [(1 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 3.5 - x2l = 0.0 - x2r = 3.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver10 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver10 = - ( - Min [(1 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 3.5 - x2l = 0.0 - x2r = 3.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver11 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver11 = - ( - Max [(2 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 3.5 - x2l = 0.0 - x2r = 3.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaver12 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaver12 = - ( - Min [(2 , 1)], - [ - LEQ [(1, dx1l), (2, dx2l), (3, (-1))] ((-yl) + (dx1l * x1l) + (dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, dx1r), (2, dx2r), (3, (-1))] ((-yr) + (dx1r * x1l) + (dx2r * x2l)), -- -5 - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 3.5 - x2l = 0.0 - x2r = 3.5 - dx1l = -1 - dx1r = -0.9 - dx2l = -0.9 - dx2r = -0.8 - yl = 4 - yr = 5 - -testPolyPaverTwoFs1 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs1 = - ( - Max [(1 , 1)], - [ - LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)), - LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)), - GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)), - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0, - LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -1 - f2dx1r = -0.9 - f2dx2l = -0.9 - f2dx2r = -0.8 - f2yl = 1 - f2yr = 2 - -testPolyPaverTwoFs2 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs2 = - ( - Min [(1 , 1)], - [ - LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)), - LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)), - GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)), - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0, - LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -1 - f2dx1r = -0.9 - f2dx2l = -0.9 - f2dx2r = -0.8 - f2yl = 1 - f2yr = 2 - -testPolyPaverTwoFs3 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs3 = - ( - Max [(2 , 1)], - [ - LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)), - LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)), - GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)), - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0, - LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -1 - f2dx1r = -0.9 - f2dx2l = -0.9 - f2dx2r = -0.8 - f2yl = 1 - f2yr = 2 - -testPolyPaverTwoFs4 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs4 = - ( - Min [(2 , 1)], - [ - LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)), - LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)), - GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)), - GEQ [(1, 1)] x1l, - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0, - LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -1 - f2dx1r = -0.9 - f2dx2l = -0.9 - f2dx2r = -0.8 - f2yl = 1 - f2yr = 2 - -testPolyPaverTwoFs5 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs5 = - ( - Max [(1 , 1)], - [ - LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)), - LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)), - GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)), - GEQ [(1, 1)] x1l, -- don't need variable >= 0, already assumed - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0, - LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -0.66 - f2dx1r = -0.66 - f2dx2l = -0.66 - f2dx2r = -0.66 - f2yl = 3 - f2yr = 4 - -testPolyPaverTwoFs6 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs6 = - ( - Min [(1 , 1)], - [ - LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)), - LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)), - GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)), - GEQ [(1, 1)] x1l, -- don't need variable >= 0, already assumed - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0, - LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -0.66 - f2dx1r = -0.66 - f2dx2l = -0.66 - f2dx2r = -0.66 - f2yl = 3 - f2yr = 4 - -testPolyPaverTwoFs7 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs7 = - ( - Max [(2 , 1)], - [ - LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)), - LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)), - GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)), - GEQ [(1, 1)] x1l, -- don't need variable >= 0, already assumed - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0, - LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -0.66 - f2dx1r = -0.66 - f2dx2l = -0.66 - f2dx2r = -0.66 - f2yl = 3 - f2yr = 4 - -testPolyPaverTwoFs8 :: (ObjectiveFunction, [PolyConstraint]) -testPolyPaverTwoFs8 = - ( - Min [(2 , 1)], - [ - LEQ [(1, f1dx1l), (2, f1dx2l), (3, (-1))] ((-f1yl) + (f1dx1l * x1l) + (f1dx2l * x2l)), -- -4, This will need an artificial variable - GEQ [(1, f1dx1r), (2, f1dx2r), (3, (-1))] ((-f1yr) + (f1dx1r * x1l) + (f1dx2r * x2l)), - LEQ [(1, f2dx1l), (2, f2dx2l), (4, (-1))] ((-f2yl) + (f2dx1l * x1l) + (f2dx2l * x2l)), - GEQ [(1, f2dx1r), (2, f2dx2r), (4, (-1))] ((-f2yr) + (f2dx1r * x1l) + (f2dx2r * x2l)), - GEQ [(1, 1)] x1l, -- don't need variable >= 0, already assumed - LEQ [(1, 1)] x1r, - GEQ [(2, 1)] x2l, - LEQ [(2, 1)] x2r, - LEQ [(3, 1)] 0, - LEQ [(4, 1)] 0 - ] - ) - where - x1l = 0.0 - x1r = 2.5 - x2l = 0.0 - x2r = 2.5 - f1dx1l = -1 - f1dx1r = -0.9 - f1dx2l = -0.9 - f1dx2r = -0.8 - f1yl = 4 - f1yr = 5 - f2dx1l = -0.66 - f2dx1r = -0.66 - f2dx2l = -0.66 - f2dx2r = -0.66 - f2yl = 3 - f2yr = 4 - --- Test cases produced by old simplex-haskell/SoPlex QuickCheck prop - --- SoPlex gives -400 for the following system but -370 is the optimized solution --- simplex-haskell gives -370 --- SoPlex gives -370 if we simplify the system before sending it to SoPlex -testQuickCheck1 = - ( - Max [(1, -6), (1, -8), (1, 9), (1, 10), (1, 8), (2, -15), (1, 13), (1, -14), (2, 0)], - [ - EQ [(1, 5), (1, 6), (2, -2), (1, 7), (1, 6), (2, 0)] (-12), - GEQ [(1, 11), (1, 0), (1, -5), (1, -12), (1, -14), (2, 11)] (-7), - GEQ [(1, -12), (1, -7), (1, -2), (2, -9), (1, 3), (1, 5), (1, -15), (2, 14)] (-8), GEQ [(1, 13), (1, 1), (1, -11), (2, 0)] 5, - LEQ [(1, -10), (1, -14), (1, 4), (1, -2), (1, -10), (1, -5), (1, -11)] (-1) - ] - ) - --- If we do not call simplifyPolyConstraints before we start the simplex algorithm, the following return a wrong solution --- Correct solution is -2/9 -testQuickCheck2 = - ( - Max [(1, -3), (2, 5)], - [ - LEQ [(2, -1), (1, -6), (2, 7)] 4, - LEQ [(1, 1), (2, -4), (3, 3)] (-2), - LEQ [(2, 6), (1, -4), (2, 1)] 0] - ) - --- This test will fail if the objective function is not simplified -testQuickCheck3 = - ( - Min [(2, 0), (2, -4)], - [ - GEQ [(1, 5), (2, 4)] (-4), - LEQ [(1, -1), (2, -1)] 2, - LEQ [(2, 1)] 2, - GEQ [(1, -5), (2, -1), (2, 1)] (-5) +import qualified Data.Map as M +import Data.Ratio +import Linear.Simplex.Types +import Prelude hiding (EQ) + +testsList :: [((ObjectiveFunction, [PolyConstraint]), Maybe Result)] +testsList = + [ (test1, Just (Result 7 (M.fromList [(7, 29), (1, 3), (2, 4)]))) + , (test2, Just (Result 7 (M.fromList [(7, 0)]))) + , (test3, Nothing) + , (test4, Just (Result 11 (M.fromList [(11, 237 % 7), (1, 24 % 7), (2, 33 % 7)]))) + , (test5, Just (Result 9 (M.fromList [(9, 3 % 5), (2, 14 % 5), (3, 17 % 5)]))) + , (test6, Nothing) + , (test7, Just (Result 8 (M.fromList [(8, 1), (2, 2), (1, 3)]))) + , (test8, Just (Result 8 (M.fromList [(8, (-1) % 4), (2, 9 % 2), (1, 17 % 4)]))) + , (test9, Just (Result 7 (M.fromList [(7, 5), (3, 2), (4, 1)]))) + , (test10, Just (Result 7 (M.fromList [(7, 8), (1, 2), (2, 6)]))) + , (test11, Just (Result 8 (M.fromList [(8, 20), (4, 16), (3, 6)]))) + , (test12, Just (Result 8 (M.fromList [(8, 6), (4, 2), (5, 2)]))) + , (test13, Just (Result 6 (M.fromList [(6, 150), (2, 150)]))) + , (test14, Just (Result 6 (M.fromList [(6, 40 % 3), (2, 40 % 3)]))) + , (test15, Nothing) + , (test16, Just (Result 6 (M.fromList [(6, 75), (1, 75 % 2)]))) + , (test17, Just (Result 7 (M.fromList [(7, (-120)), (1, 20)]))) + , (test18, Just (Result 7 (M.fromList [(7, 10), (3, 5)]))) + , (test19, Nothing) + , (test20, Nothing) + , (test21, Just (Result 7 (M.fromList [(7, 250), (2, 50)]))) + , (test22, Just (Result 7 (M.fromList [(7, 0)]))) + , (test23, Nothing) + , (test24, Just (Result 10 (M.fromList [(10, 300), (3, 150)]))) + , (test25, Just (Result 3 (M.fromList [(3, 15), (1, 15)]))) + , (test26, Just (Result 6 (M.fromList [(6, 20), (1, 10), (2, 10)]))) + , (test27, Just (Result 3 (M.fromList [(3, 0)]))) + , (test28, Just (Result 6 (M.fromList [(6, 0), (2, 10)]))) + , (test29, Nothing) + , (test30, Nothing) + , (test31, Just (Result 5 (M.fromList [(2, 1 % 1), (5, 0 % 1)]))) + , (test32, Nothing) + , (testPolyPaver1, Just (Result 12 (M.fromList [(12, 7 % 4), (2, 5 % 2), (1, 7 % 4), (3, 0)]))) + , (testPolyPaver2, Just (Result 12 (M.fromList [(12, 5 % 2), (2, 5 % 3), (1, 5 % 2), (3, 0)]))) + , (testPolyPaver3, Just (Result 12 (M.fromList [(12, 5 % 3), (2, 5 % 3), (1, 5 % 2), (3, 0)]))) + , (testPolyPaver4, Just (Result 12 (M.fromList [(12, 5 % 2), (2, 5 % 2), (1, 5 % 2), (3, 0)]))) + , (testPolyPaver5, Nothing) + , (testPolyPaver6, Nothing) + , (testPolyPaver7, Nothing) + , (testPolyPaver8, Nothing) + , (testPolyPaver9, Just (Result 12 (M.fromList [(12, 7 % 2), (2, 5 % 9), (1, 7 % 2), (3, 0)]))) + , (testPolyPaver10, Just (Result 12 (M.fromList [(12, 17 % 20), (2, 7 % 2), (1, 17 % 20), (3, 0)]))) + , (testPolyPaver11, Just (Result 12 (M.fromList [(12, 7 % 2), (2, 7 % 2), (1, 22 % 9)]))) + , (testPolyPaver12, Just (Result 12 (M.fromList [(12, 5 % 9), (2, 5 % 9), (1, 7 % 2), (3, 0)]))) + , (testPolyPaverTwoFs1, Nothing) + , (testPolyPaverTwoFs2, Nothing) + , (testPolyPaverTwoFs3, Nothing) + , (testPolyPaverTwoFs4, Nothing) + , (testPolyPaverTwoFs5, Just (Result 17 (M.fromList [(17, 5 % 2), (2, 45 % 22), (1, 5 % 2), (4, 0)]))) + , (testPolyPaverTwoFs6, Just (Result 17 (M.fromList [(17, 45 % 22), (2, 5 % 2), (1, 45 % 22), (4, 0)]))) + , (testPolyPaverTwoFs7, Just (Result 17 (M.fromList [(17, 5 % 2), (2, 5 % 2), (1, 5 % 2), (4, 0)]))) + , (testPolyPaverTwoFs8, Just (Result 17 (M.fromList [(17, 45 % 22), (2, 45 % 22), (1, 5 % 2), (4, 0)]))) + , (testLeqGeqBugMin1, Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)]))) + , (testLeqGeqBugMax1, Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)]))) + , (testLeqGeqBugMin2, Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)]))) + , (testLeqGeqBugMax2, Just (Result 5 (M.fromList [(5, 3), (1, 3), (2, 3)]))) + , (testQuickCheck1, Just (Result 10 (M.fromList [(10, (-370)), (2, 26), (1, 5 % 3)]))) + , (testQuickCheck2, Just (Result 8 (M.fromList [(8, (-2) % 9), (1, 14 % 9), (2, 8 % 9)]))) + , (testQuickCheck3, Just (Result 7 (M.fromList [(7, (-8)), (2, 2)]))) + ] + +testLeqGeqBugMin1 :: (ObjectiveFunction, [PolyConstraint]) +testLeqGeqBugMin1 = + ( Min (M.fromList [(1, 1)]) + , + [ GEQ (M.fromList [(1, 1)]) 3 + , LEQ (M.fromList [(1, 1)]) 3 + , GEQ (M.fromList [(2, 1)]) 3 + , LEQ (M.fromList [(2, 1)]) 3 + ] + ) + +testLeqGeqBugMax1 :: (ObjectiveFunction, [PolyConstraint]) +testLeqGeqBugMax1 = + ( Min (M.fromList [(1, 1)]) + , + [ GEQ (M.fromList [(1, 1)]) 3 + , LEQ (M.fromList [(1, 1)]) 3 + , GEQ (M.fromList [(2, 1)]) 3 + , LEQ (M.fromList [(2, 1)]) 3 + ] + ) + +testLeqGeqBugMin2 :: (ObjectiveFunction, [PolyConstraint]) +testLeqGeqBugMin2 = + ( Min (M.fromList [(1, 1)]) + , + [ GEQ (M.fromList [(1, 1)]) 3 + , LEQ (M.fromList [(1, 1)]) 3 + , GEQ (M.fromList [(2, 1)]) 3 + , LEQ (M.fromList [(2, 1)]) 3 + ] + ) + +testLeqGeqBugMax2 :: (ObjectiveFunction, [PolyConstraint]) +testLeqGeqBugMax2 = + ( Min (M.fromList [(1, 1)]) + , + [ GEQ (M.fromList [(1, 1)]) 3 + , LEQ (M.fromList [(1, 1)]) 3 + , GEQ (M.fromList [(2, 1)]) 3 + , LEQ (M.fromList [(2, 1)]) 3 + ] + ) + +-- From page 50 of 'Linear and Integer Programming Made Easy' +-- Solution: obj = 29, 1 = 3, 2 = 4, +test1 :: (ObjectiveFunction, [PolyConstraint]) +test1 = + ( Max (M.fromList [(1, 3), (2, 5)]) + , + [ LEQ (M.fromList [(1, 3), (2, 1)]) 15 + , LEQ (M.fromList [(1, 1), (2, 1)]) 7 + , LEQ (M.fromList [(2, 1)]) 4 + , LEQ (M.fromList [(1, -1), (2, 2)]) 6 + ] + ) + +test2 :: (ObjectiveFunction, [PolyConstraint]) +test2 = + ( Min (M.fromList [(1, 3), (2, 5)]) + , + [ LEQ (M.fromList [(1, 3), (2, 1)]) 15 + , LEQ (M.fromList [(1, 1), (2, 1)]) 7 + , LEQ (M.fromList [(2, 1)]) 4 + , LEQ (M.fromList [(1, -1), (2, 2)]) 6 + ] + ) + +test3 :: (ObjectiveFunction, [PolyConstraint]) +test3 = + ( Max (M.fromList [(1, 3), (2, 5)]) + , + [ GEQ (M.fromList [(1, 3), (2, 1)]) 15 + , GEQ (M.fromList [(1, 1), (2, 1)]) 7 + , GEQ (M.fromList [(2, 1)]) 4 + , GEQ (M.fromList [(1, -1), (2, 2)]) 6 + ] + ) + +test4 :: (ObjectiveFunction, [PolyConstraint]) +test4 = + ( Min (M.fromList [(1, 3), (2, 5)]) + , + [ GEQ (M.fromList [(1, 3), (2, 1)]) 15 + , GEQ (M.fromList [(1, 1), (2, 1)]) 7 + , GEQ (M.fromList [(2, 1)]) 4 + , GEQ (M.fromList [(1, -1), (2, 2)]) 6 + ] + ) + +-- From https://www.eng.uwaterloo.ca/~syde05/phase1.pdf +-- Solution: obj = 3/5, 2 = 14/5, 3 = 17/5 +-- requires two phases +test5 :: (ObjectiveFunction, [PolyConstraint]) +test5 = + ( Max (M.fromList [(1, 1), (2, -1), (3, 1)]) + , + [ LEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 + , LEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) + , LEQ (M.fromList [(1, -1), (2, 1), (3, -2)]) (-1) + ] + ) + +test6 :: (ObjectiveFunction, [PolyConstraint]) +test6 = + ( Min (M.fromList [(1, 1), (2, -1), (3, 1)]) + , + [ LEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 + , LEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) + , LEQ (M.fromList [(1, -1), (2, 1), (3, -2)]) (-1) + ] + ) + +test7 :: (ObjectiveFunction, [PolyConstraint]) +test7 = + ( Max (M.fromList [(1, 1), (2, -1), (3, 1)]) + , + [ GEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 + , GEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) + , GEQ (M.fromList [(1, -1), (2, 1), (3, -2)]) (-1) + ] + ) + +test8 :: (ObjectiveFunction, [PolyConstraint]) +test8 = + ( Min (M.fromList [(1, 1), (2, -1), (3, 1)]) + , + [ GEQ (M.fromList [(1, 2), (2, -1), (3, 2)]) 4 + , GEQ (M.fromList [(1, 2), (2, -3), (3, 1)]) (-5) + , GEQ (M.fromList [(1, -1), (2, 1), (3, -2)]) (-1) + ] + ) + +-- From page 49 of 'Linear and Integer Programming Made Easy' +-- Solution: obj = -5, 3 = 2, 4 = 1, objVar was negated so actual val is 5 wa +-- requires two phases +test9 :: (ObjectiveFunction, [PolyConstraint]) +test9 = + ( Min (M.fromList [(1, 1), (2, 1), (3, 2), (4, 1)]) + , + [ EQ (M.fromList [(1, 1), (3, 2), (4, -2)]) 2 + , EQ (M.fromList [(2, 1), (3, 1), (4, 4)]) 6 + ] + ) + +test10 :: (ObjectiveFunction, [PolyConstraint]) +test10 = + ( Max (M.fromList [(1, 1), (2, 1), (3, 2), (4, 1)]) + , + [ EQ (M.fromList [(1, 1), (3, 2), (4, -2)]) 2 + , EQ (M.fromList [(2, 1), (3, 1), (4, 4)]) 6 + ] + ) + +-- Adapted from page 52 of 'Linear and Integer Programming Made Easy' +-- Removed variables which do not appear in the system (these should be artificial variables) +-- Solution: obj = 20, 3 = 6, 4 = 16 wq +test11 :: (ObjectiveFunction, [PolyConstraint]) +test11 = + ( Max (M.fromList [(3, -2), (4, 2), (5, 1)]) + , + [ EQ (M.fromList [(3, -2), (4, 1), (5, 1)]) 4 + , EQ (M.fromList [(3, 3), (4, -1), (5, 2)]) 2 + ] + ) + +test12 :: (ObjectiveFunction, [PolyConstraint]) +test12 = + ( Min (M.fromList [(3, -2), (4, 2), (5, 1)]) + , + [ EQ (M.fromList [(3, -2), (4, 1), (5, 1)]) 4 + , EQ (M.fromList [(3, 3), (4, -1), (5, 2)]) 2 + ] + ) + +-- From page 59 of 'Linear and Integer Programming Made Easy' +-- Solution: obj = 150, 1 = 0, 2 = 150 +-- requires two phases +test13 :: (ObjectiveFunction, [PolyConstraint]) +test13 = + ( Max (M.fromList [(1, 2), (2, 1)]) + , + [ LEQ (M.fromList [(1, 4), (2, 1)]) 150 + , LEQ (M.fromList [(1, 2), (2, -3)]) (-40) + ] + ) + +test14 :: (ObjectiveFunction, [PolyConstraint]) +test14 = + ( Min (M.fromList [(1, 2), (2, 1)]) + , + [ LEQ (M.fromList [(1, 4), (2, 1)]) 150 + , LEQ (M.fromList [(1, 2), (2, -3)]) (-40) + ] + ) + +test15 :: (ObjectiveFunction, [PolyConstraint]) +test15 = + ( Max (M.fromList [(1, 2), (2, 1)]) + , + [ GEQ (M.fromList [(1, 4), (2, 1)]) 150 + , GEQ (M.fromList [(1, 2), (2, -3)]) (-40) + ] + ) + +test16 :: (ObjectiveFunction, [PolyConstraint]) +test16 = + ( Min (M.fromList [(1, 2), (2, 1)]) + , + [ GEQ (M.fromList [(1, 4), (2, 1)]) 150 + , GEQ (M.fromList [(1, 2), (2, -3)]) (-40) + ] + ) + +-- From page 59 of 'Linear and Integer Programming Made Easy' +-- Solution: obj = 120, 1 = 20, 2 = 0, 3 = 0, objVar was negated so actual val is -120 +test17 :: (ObjectiveFunction, [PolyConstraint]) +test17 = + ( Min (M.fromList [(1, -6), (2, -4), (3, 2)]) + , + [ LEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 + , LEQ (M.fromList [(2, -5), (3, 5)]) 100 + , LEQ (M.fromList [(1, 1), (3, 1), (1, 1)]) 400 + ] + ) + +test18 :: (ObjectiveFunction, [PolyConstraint]) +test18 = + ( Max (M.fromList [(1, -6), (2, -4), (3, 2)]) + , + [ LEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 + , LEQ (M.fromList [(2, -5), (3, 5)]) 100 + , LEQ (M.fromList [(1, 1), (3, 1), (1, 1)]) 400 + ] + ) + +test19 :: (ObjectiveFunction, [PolyConstraint]) +test19 = + ( Min (M.fromList [(1, -6), (2, -4), (3, 2)]) + , + [ GEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 + , GEQ (M.fromList [(2, -5), (3, 5)]) 100 + , GEQ (M.fromList [(1, 1), (3, 1), (1, 1)]) 400 + ] + ) + +test20 :: (ObjectiveFunction, [PolyConstraint]) +test20 = + ( Max (M.fromList [(1, -6), (2, -4), (3, 2)]) + , + [ GEQ (M.fromList [(1, 1), (2, 1), (3, 4)]) 20 + , GEQ (M.fromList [(2, -5), (3, 5)]) 100 + , GEQ (M.fromList [(1, 1), (3, 1), (1, 1)]) 400 + ] + ) + +-- From page 59 of 'Linear and Integer Programming Made Easy' +-- Solution: obj = 250, 1 = 0, 2 = 50, 3 = 0 +test21 :: (ObjectiveFunction, [PolyConstraint]) +test21 = + ( Max (M.fromList [(1, 3), (2, 5), (3, 2)]) + , + [ LEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 + , LEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 + , LEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 + ] + ) + +test22 :: (ObjectiveFunction, [PolyConstraint]) +test22 = + ( Min (M.fromList [(1, 3), (2, 5), (3, 2)]) + , + [ LEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 + , LEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 + , LEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 + ] + ) + +test23 :: (ObjectiveFunction, [PolyConstraint]) +test23 = + ( Max (M.fromList [(1, 3), (2, 5), (3, 2)]) + , + [ GEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 + , GEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 + , GEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 + ] + ) + +test24 :: (ObjectiveFunction, [PolyConstraint]) +test24 = + ( Min (M.fromList [(1, 3), (2, 5), (3, 2)]) + , + [ GEQ (M.fromList [(1, 5), (2, 1), (3, 4)]) 50 + , GEQ (M.fromList [(1, 1), (2, -1), (3, 1)]) 150 + , GEQ (M.fromList [(1, 2), (2, 1), (3, 2)]) 100 + ] + ) + +test25 :: (ObjectiveFunction, [PolyConstraint]) +test25 = + ( Max (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, 1)]) 15 + ] + ) + +test26 :: (ObjectiveFunction, [PolyConstraint]) +test26 = + ( Max (M.fromList [(1, 2)]) + , + [ LEQ (M.fromList [(1, 2)]) 20 + , GEQ (M.fromList [(2, 1)]) 10 + ] + ) + +test27 :: (ObjectiveFunction, [PolyConstraint]) +test27 = + ( Min (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, 1)]) 15 + ] + ) + +test28 :: (ObjectiveFunction, [PolyConstraint]) +test28 = + ( Min (M.fromList [(1, 2)]) + , + [ LEQ (M.fromList [(1, 2)]) 20 + , GEQ (M.fromList [(2, 1)]) 10 + ] + ) + +test29 :: (ObjectiveFunction, [PolyConstraint]) +test29 = + ( Max (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, 1)]) 15 + , GEQ (M.fromList [(1, 1)]) 15.01 + ] + ) + +test30 :: (ObjectiveFunction, [PolyConstraint]) +test30 = + ( Max (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, 1)]) 15 + , GEQ (M.fromList [(1, 1)]) 15.01 + , GEQ (M.fromList [(2, 1)]) 10 + ] + ) + +test31 :: (ObjectiveFunction, [PolyConstraint]) +test31 = + ( Min (M.fromList [(1, 1)]) + , + [ GEQ (M.fromList [(1, 1), (2, 1)]) 1 + , GEQ (M.fromList [(1, 1), (2, 1)]) 1 + ] + ) + +test32 :: (ObjectiveFunction, [PolyConstraint]) +test32 = + ( Min (M.fromList [(1, 1)]) + , + [ GEQ (M.fromList [(1, 1), (2, 1)]) 2 + , LEQ (M.fromList [(1, 1), (2, 1)]) 1 + ] + ) + +-- Tests for systems similar to those from PolyPaver2 +testPolyPaver1 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver1 = + ( Min (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaver2 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver2 = + ( Max (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaver3 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver3 = + ( Min (M.fromList [(2, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaver4 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver4 = + ( Max (M.fromList [(2, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaver5 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver5 = + ( Max (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 1.5 + x2l = 0.0 + x2r = 1.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaver6 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver6 = + ( Min (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 1.5 + x2l = 0.0 + x2r = 1.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaver7 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver7 = + ( Max (M.fromList [(2, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 1.5 + x2l = 0.0 + x2r = 1.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaver8 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver8 = + ( Min (M.fromList [(2, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 1.5 + x2l = 0.0 + x2r = 1.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaver9 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver9 = + ( Max (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 3.5 + x2l = 0.0 + x2r = 3.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaver10 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver10 = + ( Min (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 3.5 + x2l = 0.0 + x2r = 3.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaver11 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver11 = + ( Max (M.fromList [(2, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 3.5 + x2l = 0.0 + x2r = 3.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaver12 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaver12 = + ( Min (M.fromList [(2, 1)]) + , + [ LEQ (M.fromList [(1, dx1l), (2, dx2l), (3, (-1))]) (-yl + dx1l * x1l + dx2l * x2l) + , GEQ (M.fromList [(1, dx1r), (2, dx2r), (3, (-1))]) (-yr + dx1r * x1l + dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 3.5 + x2l = 0.0 + x2r = 3.5 + dx1l = -1 + dx1r = -0.9 + dx2l = -0.9 + dx2r = -0.8 + yl = 4 + yr = 5 + +testPolyPaverTwoFs1 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaverTwoFs1 = + ( Max (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) + , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) + , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) + , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + , LEQ (M.fromList [(4, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + f1dx1l = -1 + f1dx1r = -0.9 + f1dx2l = -0.9 + f1dx2r = -0.8 + f1yl = 4 + f1yr = 5 + f2dx1l = -1 + f2dx1r = -0.9 + f2dx2l = -0.9 + f2dx2r = -0.8 + f2yl = 1 + f2yr = 2 + +testPolyPaverTwoFs2 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaverTwoFs2 = + ( Min (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) + , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) + , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) + , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + , LEQ (M.fromList [(4, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + f1dx1l = -1 + f1dx1r = -0.9 + f1dx2l = -0.9 + f1dx2r = -0.8 + f1yl = 4 + f1yr = 5 + f2dx1l = -1 + f2dx1r = -0.9 + f2dx2l = -0.9 + f2dx2r = -0.8 + f2yl = 1 + f2yr = 2 + +testPolyPaverTwoFs3 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaverTwoFs3 = + ( Max (M.fromList [(2, 1)]) + , + [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) + , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) + , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) + , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + , LEQ (M.fromList [(4, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + f1dx1l = -1 + f1dx1r = -0.9 + f1dx2l = -0.9 + f1dx2r = -0.8 + f1yl = 4 + f1yr = 5 + f2dx1l = -1 + f2dx1r = -0.9 + f2dx2l = -0.9 + f2dx2r = -0.8 + f2yl = 1 + f2yr = 2 + +testPolyPaverTwoFs4 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaverTwoFs4 = + ( Min (M.fromList [(2, 1)]) + , + [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) + , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) + , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) + , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + , LEQ (M.fromList [(4, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + f1dx1l = -1 + f1dx1r = -0.9 + f1dx2l = -0.9 + f1dx2r = -0.8 + f1yl = 4 + f1yr = 5 + f2dx1l = -1 + f2dx1r = -0.9 + f2dx2l = -0.9 + f2dx2r = -0.8 + f2yl = 1 + f2yr = 2 + +testPolyPaverTwoFs5 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaverTwoFs5 = + ( Max (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) + , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) + , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) + , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + , LEQ (M.fromList [(4, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + f1dx1l = -1 + f1dx1r = -0.9 + f1dx2l = -0.9 + f1dx2r = -0.8 + f1yl = 4 + f1yr = 5 + f2dx1l = -0.66 + f2dx1r = -0.66 + f2dx2l = -0.66 + f2dx2r = -0.66 + f2yl = 3 + f2yr = 4 + +testPolyPaverTwoFs6 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaverTwoFs6 = + ( Min (M.fromList [(1, 1)]) + , + [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) + , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) + , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) + , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + , LEQ (M.fromList [(4, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + f1dx1l = -1 + f1dx1r = -0.9 + f1dx2l = -0.9 + f1dx2r = -0.8 + f1yl = 4 + f1yr = 5 + f2dx1l = -0.66 + f2dx1r = -0.66 + f2dx2l = -0.66 + f2dx2r = -0.66 + f2yl = 3 + f2yr = 4 + +testPolyPaverTwoFs7 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaverTwoFs7 = + ( Max (M.fromList [(2, 1)]) + , + [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) + , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) + , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) + , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + , LEQ (M.fromList [(4, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + f1dx1l = -1 + f1dx1r = -0.9 + f1dx2l = -0.9 + f1dx2r = -0.8 + f1yl = 4 + f1yr = 5 + f2dx1l = -0.66 + f2dx1r = -0.66 + f2dx2l = -0.66 + f2dx2r = -0.66 + f2yl = 3 + f2yr = 4 + +testPolyPaverTwoFs8 :: (ObjectiveFunction, [PolyConstraint]) +testPolyPaverTwoFs8 = + ( Min (M.fromList [(2, 1)]) + , + [ LEQ (M.fromList [(1, f1dx1l), (2, f1dx2l), (3, (-1))]) (-f1yl + f1dx1l * x1l + f1dx2l * x2l) + , GEQ (M.fromList [(1, f1dx1r), (2, f1dx2r), (3, (-1))]) (-f1yr + f1dx1r * x1l + f1dx2r * x2l) + , LEQ (M.fromList [(1, f2dx1l), (2, f2dx2l), (4, (-1))]) (-f2yl + f2dx1l * x1l + f2dx2l * x2l) + , GEQ (M.fromList [(1, f2dx1r), (2, f2dx2r), (4, (-1))]) (-f2yr + f2dx1r * x1l + f2dx2r * x2l) + , GEQ (M.fromList [(1, 1)]) x1l + , LEQ (M.fromList [(1, 1)]) x1r + , GEQ (M.fromList [(2, 1)]) x2l + , LEQ (M.fromList [(2, 1)]) x2r + , LEQ (M.fromList [(3, 1)]) 0 + , LEQ (M.fromList [(4, 1)]) 0 + ] + ) + where + x1l = 0.0 + x1r = 2.5 + x2l = 0.0 + x2r = 2.5 + f1dx1l = -1 + f1dx1r = -0.9 + f1dx2l = -0.9 + f1dx2r = -0.8 + f1yl = 4 + f1yr = 5 + f2dx1l = -0.66 + f2dx1r = -0.66 + f2dx2l = -0.66 + f2dx2r = -0.66 + f2yl = 3 + f2yr = 4 + +-- Test cases produced by old simplex-haskell/SoPlex QuickCheck prop + +testQuickCheck1 :: (ObjectiveFunction, [PolyConstraint]) +testQuickCheck1 = + ( Max (M.fromList [(1, 12), (2, -15)]) + , + [ EQ (M.fromList [(1, 24), (2, -2)]) (-12) + , GEQ (M.fromList [(1, -20), (2, 11)]) (-7) + , GEQ (M.fromList [(1, -28), (2, 5)]) (-8) + , GEQ (M.fromList [(1, 3), (2, 0)]) 5 + , LEQ (M.fromList [(1, -48)]) (-1) + ] + ) + +-- Correct solution is -2/9 +testQuickCheck2 :: (ObjectiveFunction, [PolyConstraint]) +testQuickCheck2 = + ( Max (M.fromList [(1, -3), (2, 5)]) + , + [ LEQ (M.fromList [(1, -6), (2, 6)]) 4 + , LEQ (M.fromList [(1, 1), (2, -4), (3, 3)]) (-2) + , LEQ (M.fromList [(2, 7), (1, -4)]) 0 + ] + ) + +-- This test will fail if the objective function is not simplified +testQuickCheck3 :: (ObjectiveFunction, [PolyConstraint]) +testQuickCheck3 = + ( Min (M.fromList [(2, 0), (2, -4)]) + , + [ GEQ (M.fromList [(1, 5), (2, 4)]) (-4) + , LEQ (M.fromList [(1, -1), (2, -1)]) 2 + , LEQ (M.fromList [(2, 1)]) 2 + , GEQ (M.fromList [(1, -5), (2, -1), (2, 1)]) (-5) ] )