fix(simplex): support zero objectives after phase one

ChangeLog.md

@@ -2,6 +2,7 @@
 
 ## Unreleased changes
 
+- Support zero objectives after the artificial-variable phase and replace custom coefficient-map merging with standard Data.Map operations. ([#18](https://github.com/rasheedja/simplex-method/pull/18))
 - Correct fractional coefficients in pretty-printed expressions and remove trailing plus separators. ([#19](https://github.com/rasheedja/simplex-method/pull/19))
 - Share shift and split coefficient-map transformations between objectives and constraints, and simplify access to their fields. ([#23](https://github.com/rasheedja/simplex-method/pull/23))
   - Avoid ambiguous constraint record updates using the existing field lenses, group the substitution helpers together, and test them directly.

src/Linear/Simplex/Util.hs

@@ -14,7 +14,6 @@ import Control.Monad.Logger (LogLevel (..), MonadLogger, logDebug, logError, log
 import Data.Generics.Labels ()
 import Data.List (nub, (\\))
 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)
@@ -117,37 +116,11 @@ tableauInDictionaryForm =
 
 -- | 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
+combineVarLitMapSums = Map.unionWith (+)
 
+-- | Sum coefficient maps, treating an empty list as the zero expression.
 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
-
-      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 foldVarLitMap"
-                    )
-            )
-            combinedVars
-  in  foldVarLitMap $ combinedVarMap : vms
+foldVarLitMap = Map.unionsWith (+)
 
 insertPivotObjectiveToDict :: PivotObjective -> Dict -> Dict
 insertPivotObjectiveToDict objective = Map.insert objective.variable (DictValue {varMapSum = objective.function, constant = objective.constant})

test/Linear/Simplex/Solver/TwoPhaseSpec.hs

@@ -5,7 +5,7 @@
 
 import Prelude hiding (EQ)
 
-import Control.Monad.Logger (LogLevel (LevelInfo), filterLogger, runStdoutLoggingT)
+import Control.Monad.Logger (LogLevel (LevelInfo), filterLogger, runNoLoggingT, runStdoutLoggingT)
 import qualified Data.Map as M
 import Data.Maybe (isJust)
 import Data.Ratio ((%))
@@ -1238,6 +1238,21 @@
           SimplexResult Nothing _ -> expectationFailure "Expected optimal but got infeasible"
           _ -> expectationFailure "Unexpected result"
 
+  describe "zero objectives after phase one" $ do
+    mapM_
+      ( \obj ->
+          it ("optimizes " ++ show obj ++ " after introducing artificial variables") $ do
+            let domains = VarDomainMap $ M.singleton 1 nonNegative
+            result <- runNoLoggingT $ twoPhaseSimplex domains [obj] [GEQ (M.singleton 1 1) 1]
+            result.feasibleSystem `shouldSatisfy` isJust
+            case result.objectiveResults of
+              [ObjectiveResult _ (Optimal values)] -> do
+                computeObjective obj values `shouldBe` 0
+                M.findWithDefault 0 1 values `shouldSatisfy` (>= 1)
+              _ -> expectationFailure $ "Unexpected result: " ++ show result
+      )
+      [Max M.empty, Min M.empty]
+
   describe "twoPhaseSimplex with empty constraint system" $ do
     describe "Single variable with boundedRange" $ do
       it "Max x₁ with 0 ≤ x₁ ≤ 10: optimal at x₁=10" $ do
@@ -2677,7 +2692,7 @@
             domainMap = VarDomainMap $ M.fromList [(1, lowerBoundOnly (-5))]
         let ([newObj], newConstraints, transforms) = preprocess [obj] domainMap constraints
         length transforms `shouldBe` 1
         case head transforms of
           Shift {..} -> do
             originalVar `shouldBe` 1
             shiftBy `shouldBe` (-5)
@@ -2689,7 +2704,7 @@
             domainMap = VarDomainMap $ M.fromList [(1, unbounded)]
         let (_, _, transforms) = preprocess [obj] domainMap constraints
         length transforms `shouldBe` 1
         case head transforms of
           Split {..} -> originalVar `shouldBe` 1
           _ -> expectationFailure "Expected Split transform"
 

test/Linear/Simplex/UtilSpec.hs

@@ -4,9 +4,8 @@ module Linear.Simplex.UtilSpec where
 
 import Prelude hiding (EQ)
 
-import Control.Exception (evaluate)
 import qualified Data.Map as M
-import Test.Hspec (Spec, anyErrorCall, describe, expectationFailure, it, shouldBe, shouldThrow)
+import Test.Hspec (Spec, describe, expectationFailure, it, shouldBe)
 import Test.QuickCheck (Positive (..), property)
 
 import Linear.Simplex.Types
@@ -243,8 +242,8 @@ spec = do
             m3 = M.fromList [(2, 7), (3, 5)]
         foldVarLitMap [m1, m2, m3] `shouldBe` M.fromList [(1, 3), (2, 10), (3, 5)]
 
-      it "throws error on empty list" $ do
-        evaluate (foldVarLitMap []) `shouldThrow` anyErrorCall
+      it "returns the zero expression for an empty list" $ do
+        foldVarLitMap [] `shouldBe` M.empty
 
     describe "Properties" $ do
       it "folding a singleton list is identity" $