feat (wip): VarDomain supports optional upper and lower bounds

Author: rasheedja

ChangeLog.md

@@ -2,8 +2,23 @@
 
 ## Unreleased changes
 
+- **BREAKING CHANGE**: Restructured `VarDomain` type to support upper bounds
+  - Replaced `NonNegative`, `LowerBound SimplexNum`, and `Unbounded` constructors with 
+    a single `Bounded { lowerBound :: Maybe SimplexNum, upperBound :: Maybe SimplexNum }` record
+  - Added smart constructors for convenience: `unbounded`, `nonNegative`, `lowerBoundOnly`, 
+    `upperBoundOnly`, and `boundedRange`
+  - `Bounded Nothing Nothing` is equivalent to `Unbounded`
+  - `Bounded (Just 0) Nothing` is equivalent to `NonNegative`
+  - Upper bounds are now supported and automatically added as LEQ constraints
+- Added `AddUpperBound` constructor to `VarTransform` for upper bound constraint generation
+- Updated `getTransform` to return a list of transforms (can now generate both lower and upper bound transforms)
 - Use Hspec for tests
 - Add nix flake
+- twoPhaseSimplex now takes a VarDomainMap (as the first param)
+  - You can specify each Var's domain using smart constructors: `nonNegative`, `unbounded`, 
+    `lowerBoundOnly`, `upperBoundOnly`, or `boundedRange`
+  - If a VarDomain for a Var is undefined, it's assumed to be `unbounded`
+  - If you want to keep the same behaviour as before (all vars non-negative), use `nonNegative` for all Vars
 
 ## [v0.2.0.0](https://github.com/rasheedja/LPPaver/tree/v0.2.0.0)
 

src/Linear/Simplex/Solver/TwoPhase.hs

@@ -49,6 +49,7 @@ import qualified Data.Set as Set
 import GHC.Real (Ratio)
 import Linear.Simplex.Types
 import Linear.Simplex.Util
+import qualified Control.Applicative as LPPaver
 
 -- | 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.
@@ -402,7 +403,9 @@ optimizeFeasibleSystem objFunction fsys@(FeasibleSystem {dict = phase1Dict, ..})
 -- Variables not in the VarDomainMap are assumed to be Unbounded (no lower bound).
 -- This function applies necessary transformations before solving and unapplies them after.
 -- The returned Result contains variable values and objective value in the original space.
--- TODO: use this as twoPhaseSimplex, add instructions in CHANGELOG for old users
+-- TODO: we need to be able to support multiple objective functions for the LPPaver.
+-- one way to do this is to have a list of objective functions and optimize them one by one.
+-- think about cases where the opitmal result is infinity
 twoPhaseSimplex :: (MonadIO m, MonadLogger m) => VarDomainMap -> ObjectiveFunction -> [PolyConstraint] -> m (Maybe Result)
 twoPhaseSimplex domainMap objFunction constraints = do
   logMsg LevelInfo $
@@ -489,26 +492,40 @@ collectAllVars objFunction constraints =
 -- Returns updated (transforms, nextFreshVar).
 generateTransform :: M.Map Var VarDomain -> Var -> ([VarTransform], Var) -> ([VarTransform], Var)
 generateTransform domainMap var (transforms, nextFreshVar) =
-  let domain = M.findWithDefault Unbounded var domainMap
-  in case getTransform nextFreshVar var domain of
-       Nothing -> (transforms, nextFreshVar)
-       Just t@(AddLowerBound {}) -> (t : transforms, nextFreshVar)
-       Just t@(Shift {}) -> (t : transforms, nextFreshVar + 1)
-       Just t@(Split {}) -> (t : transforms, nextFreshVar + 2)
-
--- | Determine what transform (if any) is needed for a variable given its domain.
-getTransform :: Var -> Var -> VarDomain -> Maybe VarTransform
-getTransform nextFreshVar var domain =
-  case domain of
-    NonNegative -> Nothing
-    
-    LowerBound l 
-      | l == 0 -> Nothing
-      | l > 0  -> Just $ AddLowerBound var l
-      | otherwise -> Just $ Shift var nextFreshVar l  -- l < 0, need to shift
-    
-    Unbounded ->
-      Just $ Split var nextFreshVar (nextFreshVar + 1)
+  let domain = M.findWithDefault unbounded var domainMap
+      (newTransforms, varOffset) = getTransform nextFreshVar var domain
+  in (newTransforms ++ transforms, nextFreshVar + varOffset)
+
+-- | Determine what transforms are needed for a variable given its domain.
+-- Returns a list of transforms and the number of fresh variables consumed.
+getTransform :: Var -> Var -> VarDomain -> ([VarTransform], Var)
+getTransform nextFreshVar var (Bounded mLower mUpper) =
+  let -- Handle lower bound
+      (lowerTransforms, varOffset) = case mLower of
+        Nothing -> ([], 0)  -- No lower bound: will need Split
+        Just l
+          | l == 0    -> ([], 0)  -- NonNegative: no transform needed
+          | l > 0     -> ([AddLowerBound var l], 0)  -- Positive lower bound: add constraint
+          | otherwise -> ([Shift var nextFreshVar l], 1)  -- Negative lower bound: shift
+      
+      -- Handle upper bound (if present)
+      upperTransforms = case mUpper of
+        Nothing -> []
+        Just u  -> [AddUpperBound var u]
+      
+      -- If no lower bound (Nothing), we need Split transformation
+      -- Split replaces the variable, so upper bound would apply to the original var
+      -- which gets expressed as posVar - negVar
+      (finalTransforms, finalOffset) = case mLower of
+        Nothing -> 
+          -- Unbounded: split the variable
+          -- Note: upperTransforms will still be added and will apply to the original variable
+          -- expression (posVar - negVar) via the constraint system
+          (Split var nextFreshVar (nextFreshVar + 1) : upperTransforms, 2)
+        Just _ ->
+          (lowerTransforms ++ upperTransforms, varOffset)
+  
+  in (finalTransforms, finalOffset)
 
 -- | Apply all transforms to the objective function and constraints.
 applyTransforms :: [VarTransform] -> ObjectiveFunction -> [PolyConstraint] -> (ObjectiveFunction, [PolyConstraint])
@@ -523,6 +540,10 @@ applyTransform transform (objFunction, constraints) =
     AddLowerBound v bound ->
       (objFunction, GEQ (M.singleton v 1) bound : constraints)
 
+    -- AddUpperBound: Add a LEQ constraint for the variable
+    AddUpperBound v bound ->
+      (objFunction, LEQ (M.singleton v 1) bound : constraints)
+    
     -- Shift: originalVar = shiftedVar + shiftBy (where shiftBy < 0)
     -- Substitute: wherever we see originalVar, replace with shiftedVar
     -- and adjust the RHS by -coeff * shiftBy
@@ -624,6 +645,9 @@ unapplyTransform transform result@(Result {varValMap = valMap, ..}) =
     -- AddLowerBound: No variable substitution was done, nothing to unapply
     AddLowerBound {} -> result
 
+    -- AddUpperBound: No variable substitution was done, nothing to unapply
+    AddUpperBound {} -> result
+    
     -- Shift: originalVar = shiftedVar + shiftBy
     -- So originalVar's value = shiftedVar's value + shiftBy
     Shift origVar shiftedVar shiftBy ->

src/Linear/Simplex/Types.hs

@@ -122,16 +122,45 @@ data PivotObjective = PivotObjective
   }
   deriving (Show, Read, Eq, Generic)
 
--- | Domain specification for a variable's lower bound.
--- Note: This only concerns lower bounds. Upper bounds are handled via constraints.
--- Variables not in the VarDomainMap are assumed to be Unbounded.
-data VarDomain 
-  = NonNegative           -- ^ var >= 0 (standard simplex assumption, no transformation needed)
-  | LowerBound SimplexNum -- ^ var >= L for some L (if L < 0: shift, if L > 0: add constraint)
-  | Unbounded             -- ^ No lower bound (split into difference of two non-negative vars)
-  -- TODO: Upperbound can still be useful, can negate it to get a loewr bound, can add it to the constraints
+-- | Domain specification for a variable's bounds.
+-- Variables not in the VarDomainMap are assumed to be Unbounded (both bounds Nothing).
+-- 
+-- Bounds semantics:
+--   * @lowerBound = Just L@ means var >= L
+--   * @lowerBound = Nothing@ means no lower bound (var can be arbitrarily negative)
+--   * @upperBound = Just U@ means var <= U  
+--   * @upperBound = Nothing@ means no upper bound (var can be arbitrarily positive)
+--
+-- Note: @Bounded Nothing Nothing@ is equivalent to unbounded. Use the smart constructors
+-- ('unbounded', 'nonNegative', etc.) for clarity.
+data VarDomain = Bounded 
+  { lowerBound :: Maybe SimplexNum  -- ^ Lower bound (Nothing = -∞)
+  , upperBound :: Maybe SimplexNum  -- ^ Upper bound (Nothing = +∞)
+  }
   deriving stock (Show, Read, Eq, Generic)
 
+-- | Smart constructor for an unbounded variable (no lower or upper bound).
+-- The variable can take any real value.
+unbounded :: VarDomain
+unbounded = Bounded Nothing Nothing
+
+-- | Smart constructor for a non-negative variable (var >= 0).
+-- This is the standard simplex assumption.
+nonNegative :: VarDomain
+nonNegative = Bounded (Just 0) Nothing
+
+-- | Smart constructor for a variable with only a lower bound (var >= L).
+lowerBoundOnly :: SimplexNum -> VarDomain
+lowerBoundOnly l = Bounded (Just l) Nothing
+
+-- | Smart constructor for a variable with only an upper bound (var <= U).
+upperBoundOnly :: SimplexNum -> VarDomain
+upperBoundOnly u = Bounded Nothing (Just u)
+
+-- | Smart constructor for a variable with both lower and upper bounds (L <= var <= U).
+boundedRange :: SimplexNum -> SimplexNum -> VarDomain
+boundedRange l u = Bounded (Just l) (Just u)
+
 -- | Map from variables to their domain specifications.
 -- Variables not in this map are assumed to be Unbounded.
 newtype VarDomainMap = VarDomainMap { unVarDomainMap :: M.Map Var VarDomain }
@@ -143,6 +172,10 @@ data VarTransform
       { var :: !Var
       , bound :: !SimplexNum 
       } -- ^ var >= bound where bound > 0. Adds GEQ constraint to system.
+  | AddUpperBound
+      { var :: !Var
+      , bound :: !SimplexNum
+      } -- ^ var <= bound. Adds LEQ constraint to system.
   | Shift 
       { originalVar :: !Var
       , shiftedVar :: !Var