fix(simplex): handle phase-one variables and constant constraints

ChangeLog.md

@@ -2,6 +2,7 @@
 
 ## Unreleased changes
 
+- Fix phase-one variable allocation to reserve objective-only variables, and evaluate constant constraints without crashing. ([#17](https://github.com/rasheedja/simplex-method/pull/17))
 - `twoPhaseSimplex` now takes a `VarDomainMap` as its first argument
   - Specify each variable's domain using smart constructors: `nonNegative`, `unbounded`, `lowerBoundOnly`, `upperBoundOnly`, or `boundedRange`
   - Variables not in the `VarDomainMap` are assumed to be `unbounded`

src/Linear/Simplex/Solver/TwoPhase.hs

@@ -85,7 +85,28 @@ import Linear.Simplex.Util
 --  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
+findFeasibleSolution = findFeasibleSolutionWithReservedVars Set.empty
+
+-- Reserve objective variables before allocating slack, artificial, and objective IDs.
+findFeasibleSolutionWithReservedVars ::
+  (MonadIO m, MonadLogger m) => Set Var -> [PolyConstraint] -> m (Maybe FeasibleSystem)
+findFeasibleSolutionWithReservedVars reservedVars constraints
+  | any isContradictoryConstant constraints = pure Nothing
+  | otherwise = findFeasibleVariableSystem maxVar variableConstraints
+  where
+    variableConstraints = filter (any (/= 0) . (.lhs)) constraints
+    maxVar = fromMaybe 0 $ Set.lookupMax (reservedVars <> collectAllVars [] variableConstraints)
+
+    isContradictoryConstant constraint =
+      all (== 0) constraint.lhs && case constraint of
+        LEQ _ rhs -> 0 > rhs
+        GEQ _ rhs -> 0 < rhs
+        EQ _ rhs -> 0 /= rhs
+
+-- Constant rows have already been evaluated, and maxVar covers every reserved ID.
+findFeasibleVariableSystem ::
+  (MonadIO m, MonadLogger m) => Var -> [PolyConstraint] -> m (Maybe FeasibleSystem)
+findFeasibleVariableSystem maxVar unsimplifiedSystem = do
   logMsg LevelInfo $ "findFeasibleSolution: Looking for solution for " <> showT unsimplifiedSystem
   if null artificialVars -- No artificial vars, we have a feasible system
     then do
@@ -175,19 +196,6 @@ findFeasibleSolution unsimplifiedSystem = do
   where
     system = simplifySystem unsimplifiedSystem
 
-    maxVar =
-      if null system
-        then 0
-        else
-          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
@@ -455,7 +463,8 @@ twoPhaseSimplex domainMap objFunctions constraints = do
       <> showT transformedObjs
       <> "; Transformed constraints: "
       <> showT transformedConstraints
-  mFeasibleSystem <- findFeasibleSolution transformedConstraints
+  mFeasibleSystem <-
+    findFeasibleSolutionWithReservedVars (collectAllVars transformedObjs transformedConstraints) transformedConstraints
   case mFeasibleSystem of
     Nothing -> do
       logMsg LevelInfo "twoPhaseSimplex: No feasible solution found in phase 1"

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 ((%))
@@ -2221,6 +2221,37 @@
             (x1 + x2) `shouldSatisfy` (>= 0)
           _ -> expectationFailure "Unexpected result format"
 
+  describe "constant constraints" $ do
+    let tautologies = [LEQ M.empty 1, LEQ M.empty 0, GEQ M.empty (-1), GEQ M.empty 0, EQ M.empty 0]
+        contradictions = [LEQ M.empty (-1), GEQ M.empty 1, EQ M.empty 1, EQ M.empty (-1)]
+    mapM_
+      ( \constraint ->
+          it ("accepts constant tautology " ++ show constraint) $ do
+            result <- runNoLoggingT $ twoPhaseSimplex (VarDomainMap M.empty) [] [constraint]
+            result.feasibleSystem `shouldSatisfy` isJust
+      )
+      tautologies
+    mapM_
+      ( \constraint ->
+          it ("rejects constant contradiction " ++ show constraint) $ do
+            result <- runNoLoggingT $ twoPhaseSimplex (VarDomainMap M.empty) [] [constraint]
+            result `shouldBe` SimplexResult Nothing []
+      )
+      contradictions
+
+    it "optimizes variable constraints alongside constant tautologies" $ do
+      let obj = Max (M.singleton 1 1)
+          domains = VarDomainMap $ M.singleton 1 (boundedRange 0 3)
+      result <- runNoLoggingT $ twoPhaseSimplex domains [obj] (LEQ (M.singleton 1 1) 2 : tautologies)
+      case result.objectiveResults of
+        [ObjectiveResult _ (Optimal values)] -> computeObjective obj values `shouldBe` 2
+        _ -> expectationFailure $ "Unexpected result: " ++ show result
+
+    it "recognizes explicitly zero coefficients as constant constraints" $ do
+      let domains = VarDomainMap $ M.singleton 1 nonNegative
+      result <- runNoLoggingT $ twoPhaseSimplex domains [] [LEQ (M.singleton 1 0) (-1)]
+      result `shouldBe` SimplexResult Nothing []
+
   describe "twoPhaseSimplex edge cases and infeasibility" $ do
     it "Infeasible: negative lower bound conflicts with GEQ constraint" $ do
       -- x₁ ≥ -5 (domain), but x₁ ≥ 10 and x₁ ≤ 5 (constraints conflict)
@@ -2615,7 +2646,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)
@@ -2627,7 +2658,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"
 
@@ -3084,6 +3115,30 @@
     it "shows Min objective" $ do
       prettyShowObjectiveFunction (Min (M.fromList [(1, 5)])) `shouldBe` "min: 5 * 1 + "
 
+  describe "objective-only variable allocation" $ do
+    it "keeps objective-only variables distinct from slack variables" $ do
+      let obj = Max (M.fromList [(1, 1), (2, 1)])
+          domains = VarDomainMap $ M.fromList [(1, nonNegative), (2, nonNegative)]
+      result <- runNoLoggingT $ twoPhaseSimplex domains [obj] [LEQ (M.singleton 1 1) 1]
+      result.feasibleSystem `shouldSatisfy` isJust
+      result.objectiveResults `shouldBe` [ObjectiveResult obj Unbounded]
+
+    it "reserves variables from later objectives when sharing phase one" $ do
+      let boundedObj = Max (M.singleton 1 1)
+          unboundedObj = Max (M.fromList [(1, 1), (2, 1)])
+          domains = VarDomainMap $ M.fromList [(1, nonNegative), (2, nonNegative)]
+      result <- runNoLoggingT $ twoPhaseSimplex domains [boundedObj, unboundedObj] [LEQ (M.singleton 1 1) 1]
+      case result.objectiveResults of
+        [ObjectiveResult _ (Optimal values), ObjectiveResult _ Unbounded] ->
+          computeObjective boundedObj values `shouldBe` 1
+        _ -> expectationFailure $ "Unexpected results: " ++ show result
+
+    it "reserves transformed objective-only variables" $ do
+      let obj = Max (M.fromList [(1, 1), (2, 1)])
+          domains = VarDomainMap $ M.fromList [(1, nonNegative), (2, unbounded)]
+      result <- runNoLoggingT $ twoPhaseSimplex domains [obj] [LEQ (M.singleton 1 1) 1]
+      result.objectiveResults `shouldBe` [ObjectiveResult obj Unbounded]
+
   describe "twoPhaseSimplex with multiple objectives" $ do
     it "optimizes two objectives over the same feasible region" $ do
       -- Feasible region: x₁ + x₂ ≤ 10, x₁ ≤ 6, x₂ ≤ 8, x₁,x₂ ≥ 0