@@ -34,6 +34,7 @@ module Linear.Simplex.Solver.TwoPhase
3434
3535import Prelude hiding (EQ )
3636
37+ import qualified Control.Applicative as LPPaver
3738import Control.Lens
3839import Control.Monad (unless )
3940import Control.Monad.IO.Class (MonadIO )
@@ -43,13 +44,12 @@ import Data.List
4344import qualified Data.Map as M
4445import Data.Maybe (fromJust , fromMaybe , mapMaybe )
4546import Data.Ratio (denominator , numerator , (%) )
46- import qualified Data.Text as Text
4747import Data.Set (Set )
4848import qualified Data.Set as Set
49+ import qualified Data.Text as Text
4950import GHC.Real (Ratio )
5051import Linear.Simplex.Types
5152import Linear.Simplex.Util
52- import qualified Control.Applicative as LPPaver
5353
5454-- | Find a feasible solution for the given system of 'PolyConstraint's by performing the first phase of the two-phase simplex method
5555-- All variables in the 'PolyConstraint' must be positive.
@@ -226,7 +226,6 @@ findFeasibleSolution unsimplifiedSystem = do
226226 , newArtificialVar : artificialVarsWithNewMaxVar -- Slack var is negative, r is positive (when original constraint was GEQ)
227227 )
228228 else -- r < 0
229-
230229 if basicVarCoeff <= 0 -- Should only be -1 in the standard call path
231230 then (M. insert basicVar (TableauRow {lhs = v, rhs = r}) newSystemWithoutNewMaxVar, artificialVarsWithoutNewMaxVar)
232231 else
@@ -279,21 +278,22 @@ optimizeFeasibleSystem :: (MonadIO m, MonadLogger m) => ObjectiveFunction -> Fea
279278optimizeFeasibleSystem objFunction fsys@ (FeasibleSystem {dict = phase1Dict, .. }) = do
280279 logMsg LevelInfo $
281280 " optimizeFeasibleSystem: Optimizing feasible system " <> showT fsys <> " with objective " <> showT objFunction
282- mResult <- if null artificialVars
283- then do
284- logMsg LevelInfo $
285- " optimizeFeasibleSystem: No artificial vars, system is feasible. Pivoting system (in dict form) "
286- <> showT phase1Dict
287- <> " with objective "
288- <> showT normalObjective
289- simplexPivot normalObjective phase1Dict
290- else do
291- logMsg LevelInfo $
292- " optimizeFeasibleSystem: Artificial vars present. Pivoting system (in dict form) "
293- <> showT phase1Dict
294- <> " with objective "
295- <> showT adjustedObjective
296- simplexPivot adjustedObjective phase1Dict
281+ mResult <-
282+ if null artificialVars
283+ then do
284+ logMsg LevelInfo $
285+ " optimizeFeasibleSystem: No artificial vars, system is feasible. Pivoting system (in dict form) "
286+ <> showT phase1Dict
287+ <> " with objective "
288+ <> showT normalObjective
289+ simplexPivot normalObjective phase1Dict
290+ else do
291+ logMsg LevelInfo $
292+ " optimizeFeasibleSystem: Artificial vars present. Pivoting system (in dict form) "
293+ <> showT phase1Dict
294+ <> " with objective "
295+ <> showT adjustedObjective
296+ simplexPivot adjustedObjective phase1Dict
297297 case mResult of
298298 Nothing -> do
299299 logMsg LevelInfo " optimizeFeasibleSystem: Objective is unbounded (ratio test failed)"
@@ -408,17 +408,21 @@ optimizeFeasibleSystem objFunction fsys@(FeasibleSystem {dict = phase1Dict, ..})
408408-- The returned SimplexResult contains:
409409-- * The feasible system (Nothing if infeasible)
410410-- * Results for each objective function (empty if infeasible)
411- twoPhaseSimplex :: (MonadIO m , MonadLogger m ) => VarDomainMap -> [ObjectiveFunction ] -> [PolyConstraint ] -> m SimplexResult
411+ twoPhaseSimplex ::
412+ (MonadIO m , MonadLogger m ) => VarDomainMap -> [ObjectiveFunction ] -> [PolyConstraint ] -> m SimplexResult
412413twoPhaseSimplex domainMap objFunctions constraints = do
413414 logMsg LevelInfo $
414415 " twoPhaseSimplex: Solving system with domain map " <> showT domainMap
415416 -- Collect original variables before any transformations
416417 let originalVars = collectAllVars objFunctions constraints
417418 let (transformedObjs, transformedConstraints, transforms) = preprocess objFunctions domainMap constraints
418419 logMsg LevelInfo $
419- " twoPhaseSimplex: Applied transforms " <> showT transforms
420- <> " ; Transformed objectives: " <> showT transformedObjs
421- <> " ; Transformed constraints: " <> showT transformedConstraints
420+ " twoPhaseSimplex: Applied transforms "
421+ <> showT transforms
422+ <> " ; Transformed objectives: "
423+ <> showT transformedObjs
424+ <> " ; Transformed constraints: "
425+ <> showT transformedConstraints
422426 mFeasibleSystem <- findFeasibleSolution transformedConstraints
423427 case mFeasibleSystem of
424428 Nothing -> do
@@ -435,13 +439,16 @@ twoPhaseSimplex domainMap objFunctions constraints = do
435439-- | Optimize all objective functions over the given feasible system.
436440-- Returns a list of ObjectiveResult, one for each objective function.
437441-- The originalVars set is used to filter the result to only include original decision variables.
438- optimizeAllObjectives :: (MonadIO m , MonadLogger m )
439- => Set. Set Var -- ^ Original decision variables
440- -> [VarTransform ]
441- -> FeasibleSystem
442- -> [(ObjectiveFunction , ObjectiveFunction )] -- ^ (original, transformed) objective pairs
443- -> m [ObjectiveResult ]
444- optimizeAllObjectives originalVars transforms feasibleSystem objPairs =
442+ optimizeAllObjectives ::
443+ (MonadIO m , MonadLogger m ) =>
444+ -- | Original decision variables
445+ Set. Set Var ->
446+ [VarTransform ] ->
447+ FeasibleSystem ->
448+ -- | (original, transformed) objective pairs
449+ [(ObjectiveFunction , ObjectiveFunction )] ->
450+ m [ObjectiveResult ]
451+ optimizeAllObjectives originalVars transforms feasibleSystem objPairs =
445452 mapM optimizeOne objPairs
446453 where
447454 optimizeOne (origObj, transformedObj) = do
@@ -460,22 +467,23 @@ postprocess originalVars transforms (Optimal varVals) =
460467 unappliedVarVals = unapplyTransformsToVarMap transforms varVals
461468 -- Filter to only include original decision variables
462469 filteredVarVals = M. filterWithKey (\ k _ -> Set. member k originalVars) unappliedVarVals
463- in Optimal filteredVarVals
470+ in Optimal filteredVarVals
464471
465472-- | Compute the value of an objective function given variable values.
466473computeObjective :: ObjectiveFunction -> M. Map Var SimplexNum -> SimplexNum
467474computeObjective objFunction varVals =
468475 let coeffs = case objFunction of
469476 Max m -> m
470477 Min m -> m
471- in sum $ map (\ (var, coeff) -> coeff * M. findWithDefault 0 var varVals) (M. toList coeffs)
478+ in sum $ map (\ (var, coeff) -> coeff * M. findWithDefault 0 var varVals) (M. toList coeffs)
472479
473480-- | Preprocess the system by applying variable transformations based on domain information.
474481-- Returns the transformed objectives, constraints, and the list of transforms applied.
475- preprocess :: [ObjectiveFunction ]
476- -> VarDomainMap
477- -> [PolyConstraint ]
478- -> ([ObjectiveFunction ], [PolyConstraint ], [VarTransform ])
482+ preprocess ::
483+ [ObjectiveFunction ] ->
484+ VarDomainMap ->
485+ [PolyConstraint ] ->
486+ ([ObjectiveFunction ], [PolyConstraint ], [VarTransform ])
479487preprocess objFunctions (VarDomainMap domainMap) constraints =
480488 let -- Collect all variables in the system (from all objectives and constraints)
481489 allVars = collectAllVars objFunctions constraints
@@ -488,8 +496,8 @@ preprocess objFunctions (VarDomainMap domainMap) constraints =
488496 transformedObjs = map (\ obj -> fst $ applyTransforms transforms obj constraints) objFunctions
489497 (_, transformedConstraints) = case objFunctions of
490498 [] -> (Max M. empty, applyTransformsToConstraints transforms constraints)
491- (obj: _) -> applyTransforms transforms obj constraints
492- in (transformedObjs, transformedConstraints, transforms)
499+ (obj : _) -> applyTransforms transforms obj constraints
500+ in (transformedObjs, transformedConstraints, transforms)
493501
494502-- | Apply transforms to constraints only (used when there are no objectives)
495503applyTransformsToConstraints :: [VarTransform ] -> [PolyConstraint ] -> [PolyConstraint ]
@@ -501,12 +509,12 @@ collectAllVars :: [ObjectiveFunction] -> [PolyConstraint] -> Set Var
501509collectAllVars objFunctions constraints =
502510 let objVars = Set. unions $ map getObjVars objFunctions
503511 constraintVars = Set. unions $ map getConstraintVars constraints
504- in Set. union objVars constraintVars
512+ in Set. union objVars constraintVars
505513 where
506514 getObjVars :: ObjectiveFunction -> Set Var
507515 getObjVars (Max m) = M. keysSet m
508516 getObjVars (Min m) = M. keysSet m
509-
517+
510518 getConstraintVars :: PolyConstraint -> Set Var
511519 getConstraintVars (LEQ m _) = M. keysSet m
512520 getConstraintVars (GEQ m _) = M. keysSet m
@@ -519,24 +527,24 @@ generateTransform :: M.Map Var VarDomain -> Var -> ([VarTransform], Var) -> ([Va
519527generateTransform domainMap var (transforms, nextFreshVar) =
520528 let domain = M. findWithDefault unbounded var domainMap
521529 (newTransforms, varOffset) = getTransform nextFreshVar var domain
522- in (newTransforms ++ transforms, nextFreshVar + varOffset)
530+ in (newTransforms ++ transforms, nextFreshVar + varOffset)
523531
524532-- | Determine what transforms are needed for a variable given its domain.
525533-- Returns a list of transforms and the number of fresh variables consumed.
526534getTransform :: Var -> Var -> VarDomain -> ([VarTransform ], Var )
527535getTransform nextFreshVar var (Bounded mLower mUpper) =
528536 let -- Handle lower bound
529537 (lowerTransforms, varOffset) = case mLower of
530- Nothing -> ([] , 0 ) -- No lower bound: will need Split
538+ Nothing -> ([] , 0 ) -- No lower bound: will need Split
531539 Just l
532- | l == 0 -> ([] , 0 ) -- NonNegative: no transform needed
533- | l > 0 -> ([AddLowerBound var l], 0 ) -- Positive lower bound: add constraint
534- | otherwise -> ([Shift var nextFreshVar l], 1 ) -- Negative lower bound: shift
540+ | l == 0 -> ([] , 0 ) -- NonNegative: no transform needed
541+ | l > 0 -> ([AddLowerBound var l], 0 ) -- Positive lower bound: add constraint
542+ | otherwise -> ([Shift var nextFreshVar l], 1 ) -- Negative lower bound: shift
535543
536544 -- Handle upper bound (if present)
537545 upperTransforms = case mUpper of
538546 Nothing -> []
539- Just u -> [AddUpperBound var u]
547+ Just u -> [AddUpperBound var u]
540548
541549 -- If no lower bound (Nothing), we need Split transformation
542550 -- Split replaces the variable, so upper bound would apply to the original var
@@ -549,8 +557,7 @@ getTransform nextFreshVar var (Bounded mLower mUpper) =
549557 (Split var nextFreshVar (nextFreshVar + 1 ) : upperTransforms, 2 )
550558 Just _ ->
551559 (lowerTransforms ++ upperTransforms, varOffset)
552-
553- in (finalTransforms, finalOffset)
560+ in (finalTransforms, finalOffset)
554561
555562-- | Apply all transforms to the objective function and constraints.
556563applyTransforms :: [VarTransform ] -> ObjectiveFunction -> [PolyConstraint ] -> (ObjectiveFunction , [PolyConstraint ])
@@ -564,21 +571,18 @@ applyTransform transform (objFunction, constraints) =
564571 -- AddLowerBound: Add a GEQ constraint for the variable
565572 AddLowerBound v bound ->
566573 (objFunction, GEQ (M. singleton v 1 ) bound : constraints)
567-
568574 -- AddUpperBound: Add a LEQ constraint for the variable
569575 AddUpperBound v bound ->
570576 (objFunction, LEQ (M. singleton v 1 ) bound : constraints)
571-
572577 -- Shift: originalVar = shiftedVar + shiftBy (where shiftBy < 0)
573578 -- Substitute: wherever we see originalVar, replace with shiftedVar
574579 -- and adjust the RHS by -coeff * shiftBy
575580 Shift origVar shiftedVar shiftBy ->
576581 ( applyShiftToObjective origVar shiftedVar shiftBy objFunction
577582 , map (applyShiftToConstraint origVar shiftedVar shiftBy) constraints
578583 )
579-
580584 -- Split: originalVar = posVar - negVar
581- -- Substitute: wherever we see originalVar with coeff c,
585+ -- Substitute: wherever we see originalVar with coeff c,
582586 -- replace with posVar with coeff c and negVar with coeff -c
583587 Split origVar posVar negVar ->
584588 ( applySplitToObjective origVar posVar negVar objFunction
@@ -612,13 +616,13 @@ applyShiftToConstraint origVar shiftedVar shiftBy constraint =
612616 case constraint of
613617 LEQ m rhs ->
614618 let (newMap, rhsAdjust) = substituteVarInMap origVar shiftedVar shiftBy m
615- in LEQ newMap (rhs - rhsAdjust)
619+ in LEQ newMap (rhs - rhsAdjust)
616620 GEQ m rhs ->
617621 let (newMap, rhsAdjust) = substituteVarInMap origVar shiftedVar shiftBy m
618- in GEQ newMap (rhs - rhsAdjust)
622+ in GEQ newMap (rhs - rhsAdjust)
619623 EQ m rhs ->
620624 let (newMap, rhsAdjust) = substituteVarInMap origVar shiftedVar shiftBy m
621- in EQ newMap (rhs - rhsAdjust)
625+ in EQ newMap (rhs - rhsAdjust)
622626 where
623627 substituteVarInMap :: Var -> Var -> SimplexNum -> VarLitMapSum -> (VarLitMapSum , SimplexNum )
624628 substituteVarInMap oldVar newVar shift m =
@@ -669,24 +673,21 @@ unapplyTransformToVarMap transform valMap =
669673 case transform of
670674 -- AddLowerBound: No variable substitution was done, nothing to unapply
671675 AddLowerBound {} -> valMap
672-
673676 -- AddUpperBound: No variable substitution was done, nothing to unapply
674677 AddUpperBound {} -> valMap
675-
676678 -- Shift: originalVar = shiftedVar + shiftBy
677679 -- So originalVar's value = shiftedVar's value + shiftBy
678680 Shift origVar shiftedVar shiftBy ->
679681 let shiftedVal = M. findWithDefault 0 shiftedVar valMap
680682 origVal = shiftedVal + shiftBy
681- in M. insert origVar origVal (M. delete shiftedVar valMap)
682-
683+ in M. insert origVar origVal (M. delete shiftedVar valMap)
683684 -- Split: originalVar = posVar - negVar
684685 -- So originalVar's value = posVar's value - negVar's value
685686 Split origVar posVar negVar ->
686687 let posVal = M. findWithDefault 0 posVar valMap
687688 negVal = M. findWithDefault 0 negVar valMap
688689 origVal = posVal - negVal
689- in M. insert origVar origVal (M. delete posVar (M. delete negVar valMap))
690+ in M. insert origVar origVal (M. delete posVar (M. delete negVar valMap))
690691
691692-- | Perform the simplex pivot algorithm on a system with basic vars, assume that the first row is the 'ObjectiveFunction'.
692693simplexPivot :: (MonadIO m , MonadLogger m ) => PivotObjective -> Dict -> m (Maybe Dict )
@@ -794,13 +795,13 @@ simplexPivot objective@(PivotObjective {variable = objectiveVar, function = obje
794795 dictEntertingRow
795796 & # varMapSum
796797 %~ ( \ basicEquation ->
797- -- uncurry
798- M. insert
799- leavingVariable
800- (- 1 )
801- (filterOutEnteringVarTerm basicEquation)
802- & traverse
803- %~ divideByNegatedEnteringVariableCoeff
798+ -- uncurry
799+ M. insert
800+ leavingVariable
801+ (- 1 )
802+ (filterOutEnteringVarTerm basicEquation)
803+ & traverse
804+ %~ divideByNegatedEnteringVariableCoeff
804805 )
805806 & # constant
806807 %~ divideByNegatedEnteringVariableCoeff
0 commit comments