diff --git a/ChangeLog.md b/ChangeLog.md index d47e0a4..dd86ccb 100644 --- a/ChangeLog.md +++ b/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` diff --git a/src/Linear/Simplex/Solver/TwoPhase.hs b/src/Linear/Simplex/Solver/TwoPhase.hs index ca148f0..40408dd 100644 --- a/src/Linear/Simplex/Solver/TwoPhase.hs +++ b/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" diff --git a/test/Linear/Simplex/Solver/TwoPhaseSpec.hs b/test/Linear/Simplex/Solver/TwoPhaseSpec.hs index 9ca2246..f431fec 100644 --- a/test/Linear/Simplex/Solver/TwoPhaseSpec.hs +++ b/test/Linear/Simplex/Solver/TwoPhaseSpec.hs @@ -5,7 +5,7 @@ module Linear.Simplex.Solver.TwoPhaseSpec where 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 @@ spec = do (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) @@ -3084,6 +3115,30 @@ spec = do 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