diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,4 +1,12 @@ -# Changelog for math-programming-glpk +# Changelog for `math-programming-glpk` + +## [0.5.0] -- 18 January 2023 + +### Changed + +- Make the `Glpk` monad work with threaded runtimes. The underlying + library requires that GLPK environments be accessed from a single + thread, so we run the monad in bound threads. ## [0.4.1] -- 5 July 2020 diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright Patrick Steele (c) 2018 +Copyright (c) 2018-2023, Patrick Steele All rights reserved. diff --git a/README.md b/README.md deleted file mode 100644 --- a/README.md +++ /dev/null @@ -1,5 +0,0 @@ -# math-programming-glpk: A GLPK backend to the math-programming library - -This package provides an interface to the -[GLPK](https://www.gnu.org/software/glpk/) LP and MIP solver for the -[`math-programming`](https://github.com/prsteele/math-programming) package. diff --git a/math-programming-glpk.cabal b/math-programming-glpk.cabal --- a/math-programming-glpk.cabal +++ b/math-programming-glpk.cabal @@ -1,65 +1,47 @@ -cabal-version: 1.12 -name: math-programming-glpk -version: 0.4.1 -license: BSD3 -license-file: LICENSE -copyright: 2018 Patrick Steele -maintainer: steele.pat@gmail.com -author: Patrick Steele -homepage: https://github.com/prsteele/math-programming-glpk#readme -bug-reports: https://github.com/prsteele/math-programming-glpk/issues -synopsis: A GLPK backend to the math-programming library. -description: - Please see the README on GitHub at - -category: Math -build-type: Simple -extra-source-files: - README.md - ChangeLog.md +cabal-version: 2.4 +name: math-programming-glpk +version: 0.5.0 +synopsis: A GLPK backend to the math-programming library. +description: Please see the . +bug-reports: https://github.com/prsteele/math-programming/issues +license: BSD-3-Clause +license-file: LICENSE +author: Patrick Steele +maintainer: steele.pat@gmail.com +copyright: 2018-2023, Patrick Steele +category: Math +build-type: Simple +extra-source-files: ChangeLog.md source-repository head - type: git - location: https://github.com/prsteele/math-programming-glpk + type: git + location: https://github.com/prsteele/math-programming library - exposed-modules: - Math.Programming.Glpk - Math.Programming.Glpk.Internal - - hs-source-dirs: src - other-modules: Paths_math_programming_glpk - default-language: Haskell2010 - ghc-options: -Wall - build-depends: - base >=4.7 && <5, - containers >=0.6.0.1 && <0.7, - glpk-headers >=0.4.0 && <0.5, - math-programming >=0.4.0 && <0.5, - mtl >=2.2.2 && <2.3, - text >=1.2.3.1 && <1.3 - -test-suite math-programming-glpk-test - type: exitcode-stdio-1.0 - main-is: Driver.hs - hs-source-dirs: test - other-modules: - ApiTests - RegressionTests - Paths_math_programming_glpk + default-language: Haskell2010 + hs-source-dirs: src + exposed-modules: Math.Programming.Glpk + other-modules: Math.Programming.Glpk.Internal + ghc-options: -Wall + build-depends: base <5 + , glpk-headers ^>=0.5 + , math-programming ^>=0.5 + , mtl ^>=2.2 + , text ^>=1.2 + , transformers ^>=0.5 + , unliftio ^>=0.2 - default-language: Haskell2010 - extra-libraries: glpk - ghc-options: -Wall - build-depends: - base >=4.7 && <5, - containers >=0.6.0.1 && <0.7, - glpk-headers >=0.4.0 && <0.5, - math-programming >=0.4.0 && <0.5, - math-programming-glpk -any, - math-programming-tests >=0.4.0 && <0.5, - mtl >=2.2.2 && <2.3, - tasty >=1.2.3 && <1.3, - tasty-discover >=4.2.1 && <4.3, - tasty-hunit >=0.10.0.2 && <0.11, - text >=1.2.3.1 && <1.3 +test-suite math-programming-glpk-test-suite + default-language: Haskell2010 + type: exitcode-stdio-1.0 + main-is: Spec.hs + other-modules: Math.Programming.GlpkSpec + hs-source-dirs: test + ghc-options: -threaded -with-rtsopts=-N -Wall + build-depends: math-programming-glpk + , base <5 + , math-programming ^>=0.5 + , math-programming-tests ^>=0.5 + , unliftio ^>=0.2 + , hspec ^>=2.8 + , hspec-core ^>=2.8 diff --git a/src/Math/Programming/Glpk.hs b/src/Math/Programming/Glpk.hs --- a/src/Math/Programming/Glpk.hs +++ b/src/Math/Programming/Glpk.hs @@ -2,17 +2,23 @@ -- -- This package allows both linear and mixed-integer programs to be solved. module Math.Programming.Glpk - ( Glpk - , runGlpk - , writeFormulation - -- ** Controlling GLPK behavior - -- $settings - , GlpkEnv - , GlpkError + ( Glpk, + GlpkVariable, + GlpkConstraint, + GlpkObjective, + runGlpk, + writeFormulation, + GlpkT, + MonadGlpk, + + -- ** Controlling GLPK behavior + -- $settings + GlpkEnv, + GlpkException, ) - where +where -import Math.Programming.Glpk.Internal +import Math.Programming.Glpk.Internal -- $settings -- diff --git a/src/Math/Programming/Glpk/Internal.hs b/src/Math/Programming/Glpk/Internal.hs --- a/src/Math/Programming/Glpk/Internal.hs +++ b/src/Math/Programming/Glpk/Internal.hs @@ -1,466 +1,543 @@ -{-| This module contains the full definitions backing the simplified API -exposed in 'Math.Programming.Glpk'. --} -{-# LANGUAGE DeriveFunctor #-} -{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE OverloadedStrings #-} + +-- | This module contains the full definitions backing the simplified API +-- exposed in 'Math.Programming.Glpk'. module Math.Programming.Glpk.Internal where -import Control.Exception -import Control.Monad -import Control.Monad.Except -import Control.Monad.Reader -import Data.IORef -import Data.List -import Data.Typeable -import Data.Void -import Foreign.C.String -import Foreign.C.Types -import Foreign.Marshal.Alloc -import Foreign.Ptr -import Foreign.Storable +import Control.Monad +import Control.Monad.Trans +import Control.Monad.Trans.Reader +import Control.Monad.Trans.State +import Data.Functor +import qualified Data.Text as T +import Data.Typeable +import Data.Void +import Foreign.C.String +import Foreign.C.Types +import Foreign.Marshal.Alloc +import Foreign.Ptr +import Foreign.Storable +import Math.Programming +import Math.Programming.Glpk.Header +import UnliftIO +import UnliftIO.Concurrent -import Math.Programming -import Math.Programming.Glpk.Header +-- | A reference to a GLPK variable. +type GlpkVariable = GlpkPtr Column --- | An environment to solve math programs using GLPK. -newtype Glpk a = Glpk { _runGlpk :: ExceptT GlpkError (ReaderT GlpkEnv IO) a } +-- | A reference to a GLPK constraint. +type GlpkConstraint = GlpkPtr Row + +-- | A placeholder for an objective. +-- +-- GLPK supports only single-objective problems, and so no indices +-- need to be stored. +newtype GlpkObjective = GlpkObjective () + +class (MonadLP GlpkVariable GlpkConstraint GlpkObjective m, MonadIP GlpkVariable GlpkConstraint GlpkObjective m) => MonadGlpk m where + writeFormulation :: FilePath -> m () + +-- | An interface to the low-level GLPK API. +-- +-- High-level solver settings can be modified by altering the +-- 'SimplexMethodControlParameters' and 'MIPControlParameters' values +-- for LP and IP solves, respectively. +data GlpkEnv = GlpkEnv + { -- | A pointer to the Problem object. Most GLPK routines take this + -- as the first argument. + _glpkEnvProblem :: Ptr Problem, + -- | The variables in the model + _glpkVariables :: IORef [GlpkVariable], + -- | The next unique ID to assign to a variable. + _glpkNextVariableId :: IORef Integer, + -- | The constraints in the model + _glpkConstraints :: IORef [GlpkConstraint], + -- | The next unique ID to assign to a variable. + _glpkNextConstraintId :: IORef Integer, + -- | The control parameters for the simplex method + _glpkSimplexControl :: IORef SimplexMethodControlParameters, + -- | The control parameters for the MIP solver + _glpkMIPControl :: IORef (MIPControlParameters Void), + -- | The type of the last solve. This is needed to know whether to + -- retrieve simplex, interior point, or MIP solutions. + _glpkLastSolveType :: IORef (Maybe SolveType) + } + +-- | A pointer to a GLPK row or column. +-- +-- We assign an immutable unique value to each 'GlpkPtr' we create. +-- +-- Internally, GLPK refers to variables and constraints by their +-- column and row indices, respectively. These indices can change when +-- rows and columns are deleted, so we update this value as necessary. +data GlpkPtr a = GlpkPtr + { -- | An immutable, unique value associated with this pointer. + _glpkPtrId :: Integer, + -- | The referenced object. + _glpkPtrRef :: IORef a, + -- | Whether this reference has been deleted from the problem. + _glpkPtrDeleted :: IORef Bool + } + +instance Eq (GlpkPtr a) where + (GlpkPtr x _ _) == (GlpkPtr y _ _) = x == y + +instance Ord (GlpkPtr a) where + compare (GlpkPtr x _ _) (GlpkPtr y _ _) = compare x y + +-- | An error that GLPK can encounter. +data GlpkException + = UnknownVariable + | UnknownCode T.Text CInt + | GlpkFailure T.Text deriving - ( Functor - , Applicative - , Monad - , MonadIO - , MonadReader GlpkEnv - , MonadError GlpkError + ( Show, + Typeable ) -instance LPMonad Glpk where - type Numeric Glpk = Double - - data Variable Glpk - = Variable { fromVariable :: GlpkVariable } - deriving - ( Eq - , Ord - , Show - ) +instance Exception GlpkException - data Constraint Glpk - = Constraint { fromConstraint :: GlpkConstraint } - deriving - ( Eq - , Ord - , Show - ) +-- | An environment to solve math programs using GLPK. +newtype GlpkT m a = GlpkT {_runGlpk :: ReaderT GlpkEnv m a} + deriving + ( Functor, + Applicative, + Monad, + MonadIO, + MonadUnliftIO, + MonadTrans + ) - data Objective Glpk = Objective +type Glpk = GlpkT IO +instance MonadLP GlpkVariable GlpkConstraint GlpkObjective Glpk where addVariable = addVariable' - removeVariable = removeVariable' + deleteVariable = deleteVariable' getVariableName = getVariableName' setVariableName = setVariableName' + getVariableValue = getVariableValue' getVariableBounds = getVariableBounds' setVariableBounds = setVariableBounds' - getVariableValue = getVariableValue' + addConstraint = addConstraint' - removeConstraint = removeConstraint' + deleteConstraint = deleteConstraint' getConstraintName = getConstraintName' setConstraintName = setConstraintName' - getDualValue = getDualValue' + getConstraintValue = getDualValue + addObjective = addObjective' + deleteObjective = deleteObjective' getObjectiveName = getObjectiveName' setObjectiveName = setObjectiveName' - getObjectiveSense = getObjectiveSense' - setObjectiveSense = setObjectiveSense' getObjectiveValue = getObjectiveValue' + getObjectiveSense = getSense' + setObjectiveSense = setSense' + getTimeout = getTimeout' setTimeout = setTimeout' optimizeLP = optimizeLP' -instance IPMonad Glpk where - optimizeIP = optimizeIP' +instance MonadIP GlpkVariable GlpkConstraint GlpkObjective Glpk where getVariableDomain = getVariableDomain' setVariableDomain = setVariableDomain' getRelativeMIPGap = getRelativeMIPGap' setRelativeMIPGap = setRelativeMIPGap' + optimizeIP = optimizeIP' -runGlpk :: Glpk a -> IO (Either GlpkError a) -runGlpk glpk = do - -- Turn off terminal output. If we don't, users won't be able to - -- inhibit terminal output generated from our setup. - _ <- glp_term_out glpkOff +instance MonadGlpk Glpk where + writeFormulation = writeFormulation' - bracket glp_create_prob glp_delete_prob $ \problem -> do - -- Load the default simplex control parameters - defaultSimplexControl <- alloca $ \simplexControlPtr -> do - glp_init_smcp simplexControlPtr - peek simplexControlPtr +instance MonadGlpk m => MonadGlpk (ReaderT r m) where + writeFormulation = lift . writeFormulation - -- Load the default MIP control parameters - defaultMipControl <- alloca $ \mipControlPtr -> do - glp_init_iocp mipControlPtr - peek mipControlPtr +instance MonadGlpk m => MonadGlpk (StateT s m) where + writeFormulation = lift . writeFormulation - -- Turn on presolve, because it seems insane not to. - -- - -- In particular, this ensures that a naked call to optimizeIP - -- doesn't fail because of the lack of a basis. Sophisticated users - -- can control this parameter as they see fit befor the first - -- optimization call. - let simplexControl = defaultSimplexControl { smcpPresolve = glpkPresolveOn } - mipControl = defaultMipControl { iocpPresolve = glpkPresolveOn } +withGlpkErrorHook :: (Ptr a -> IO CInt) -> Ptr a -> IO b -> IO b +withGlpkErrorHook hook ptr actions = + bracket (mkHaskellErrorHook hook) freeHaskellFunPtr $ \hookPtr -> do + glp_error_hook hookPtr ptr + actions - env <- GlpkEnv problem - <$> newIORef [] - <*> newIORef [] - <*> newIORef simplexControl - <*> newIORef mipControl - <*> newIORef Nothing +removeGlpkErrorHook :: IO () +removeGlpkErrorHook = glp_error_hook nullFunPtr nullPtr - runReaderT (runExceptT (_runGlpk glpk)) env +runGlpk :: Glpk a -> IO a +runGlpk program = + let withGlpkEnv actions = + bracket glp_init_env (const glp_free_env) $ \case + 0 -> actions + 1 -> throwIO (GlpkFailure "GLPK already initialized") + 2 -> throwIO (GlpkFailure "GLPK failed to initialize; not enough memory") + 3 -> throwIO (GlpkFailure "GLPK failed to initialize; unsupported programming model") + r -> throwIO (GlpkFailure ("GLPK failed to initialize; unknown status code " <> T.pack (show r))) + in runInBoundThread $ + withGlpkEnv $ + flip finally removeGlpkErrorHook $ + withGlpkErrorHook (const glp_free_env) nullPtr $ + runGlpk' program +getDefaultSimplexControlParameters :: IO SimplexMethodControlParameters +getDefaultSimplexControlParameters = do + params <- alloca $ \simplexControlPtr -> do + glp_init_smcp simplexControlPtr + peek simplexControlPtr -data SolveType = LP | MIP | InteriorPoint + -- Turn on presolve. Users can simply turn this off before the first + -- optimization call if they desire. + pure (params {smcpPresolve = glpkPresolveOn}) --- | An interface to the low-level GLPK API. --- --- High-level solver settings can be modified by altering the --- 'SimplexMethodControlParameters' and 'MIPControlParameters' values --- for LP and IP solves, respectively. -data GlpkEnv - = GlpkEnv - { _glpkEnvProblem :: Ptr Problem - -- ^ A pointer to the Problem object. Most GLPK routines take this - -- as the first argument. - , _glpkVariables :: IORef [GlpkVariable] - -- ^ The variables in the model - , _glpkConstraints :: IORef [GlpkConstraint] - -- ^ The constraints in the model - , _glpkSimplexControl :: IORef SimplexMethodControlParameters - -- ^ The control parameters for the simplex method - , _glpkMIPControl :: IORef (MIPControlParameters Void) - -- ^ The control parameters for the MIP solver - , _glpkLastSolveType :: IORef (Maybe SolveType) - -- ^ The type of the last solve. This is needed to know whether to - -- retrieve simplex, interior point, or MIP solutions. - } +getDefaultMIPControlParameters :: IO (MIPControlParameters Void) +getDefaultMIPControlParameters = do + params <- alloca $ \mipControlPtr -> do + glp_init_iocp mipControlPtr + peek mipControlPtr -data NamedRef a - = NamedRef - { namedRefId :: Int - , namedRefRef :: IORef a - } + -- Turn on presolve. Users can simply turn this off before the first + -- optimization call if they desire. + pure (params {iocpPresolve = glpkPresolveOn}) -instance Eq (NamedRef a) where - x == y = namedRefId x == namedRefId y +runGlpk' :: Glpk a -> IO a +runGlpk' glpk = do + -- Turn off terminal output. If we don't, users won't be able to + -- inhibit terminal output generated from our setup. + _ <- glp_term_out glpkOff -instance Ord (NamedRef a) where - x <= y = namedRefId x <= namedRefId y + bracket glp_create_prob glp_delete_prob $ \problem -> do + env <- + GlpkEnv problem + <$> newIORef [] + <*> newIORef 0 + <*> newIORef [] + <*> newIORef 0 + <*> (getDefaultSimplexControlParameters >>= newIORef) + <*> (getDefaultMIPControlParameters >>= newIORef) + <*> newIORef Nothing -instance Show (NamedRef a) where - show = show . namedRefId + runReaderT (_runGlpk glpk) env -type GlpkConstraint = NamedRef Row +data SolveType = LP | MIP | InteriorPoint -type GlpkVariable = NamedRef Column +-- | Retrieve a component of the Glpk context +askGlpk :: Monad m => (GlpkEnv -> a) -> GlpkT m a +askGlpk = flip fmap (GlpkT ask) -askProblem :: Glpk (Ptr Problem) -askProblem = asks _glpkEnvProblem +-- | The underlying Glpk problem pointer +askProblem :: Monad m => GlpkT m (Ptr Problem) +askProblem = askGlpk _glpkEnvProblem +-- | All registered variables askVariablesRef :: Glpk (IORef [GlpkVariable]) -askVariablesRef = asks _glpkVariables +askVariablesRef = askGlpk _glpkVariables +-- | All registered constraints askConstraintsRef :: Glpk (IORef [GlpkConstraint]) -askConstraintsRef = asks _glpkConstraints - -register :: Glpk (IORef [NamedRef a]) -> NamedRef a -> Glpk () -register askRef x = do - ref <- askRef - liftIO $ modifyIORef' ref (x :) +askConstraintsRef = askGlpk _glpkConstraints -unregister :: (Enum a) => Glpk (IORef [NamedRef a]) -> NamedRef a -> Glpk () -unregister askRef x = - let - decrement (NamedRef _ ref) = modifyIORef' ref pred +-- | Note that a new row or column has been added to the to the problem. +register :: GlpkPtr a -> IORef [GlpkPtr a] -> Glpk () +register newPtr ptrRefs = do + liftIO $ modifyIORef' ptrRefs (newPtr :) - mogrify [] = return () - mogrify (z: zs) | z <= x = return () - | otherwise = decrement z >> mogrify zs - in do - ref <- askRef - liftIO $ do - -- Remove the element to be unregistered - modifyIORef' ref (delete x) +-- | Note that a row or column has been deleted from the problem, and +-- update row or column indices accordingly. +unregister :: Integral a => GlpkPtr a -> IORef [GlpkPtr a] -> Glpk () +unregister deletedPtr ptrsRef = + let update removed (GlpkPtr _ ptr _) = do + z <- readIORef ptr + when (z > removed) $ + modifyIORef' ptr pred + in liftIO $ do + -- If the reference was already deleted, do nothing + deleted <- readIORef (_glpkPtrDeleted deletedPtr) + unless deleted $ do + -- Mark deletion + writeIORef (_glpkPtrDeleted deletedPtr) True - -- Modify the referenced values that were greater than the - -- referenced element - readIORef ref >>= mogrify + -- Remove the element to be unregistered + modifyIORef' ptrsRef (filter ((/= _glpkPtrId deletedPtr) . _glpkPtrId)) -data GlpkError - = UnknownVariable GlpkVariable - | UnknownCode String CInt - deriving - ( Show - ) + -- Modify the referenced values that were greater than the + -- referenced element. + deletedId <- readIORef (_glpkPtrRef deletedPtr) + ptrs <- readIORef ptrsRef + mapM_ (update deletedId) ptrs -readColumn :: Variable Glpk -> Glpk Column -readColumn = liftIO . readIORef . namedRefRef . fromVariable +readColumn :: GlpkVariable -> Glpk Column +readColumn = liftIO . readIORef . _glpkPtrRef -readRow :: Constraint Glpk -> Glpk Row -readRow = liftIO . readIORef . namedRefRef . fromConstraint +readRow :: GlpkConstraint -> Glpk Row +readRow = liftIO . readIORef . _glpkPtrRef -addVariable' :: Glpk (Variable Glpk) +addVariable' :: Glpk GlpkVariable addVariable' = do problem <- askProblem variable <- liftIO $ do column <- glp_add_cols problem 1 - glp_set_col_bnds problem column glpkFree 0 0 + GlpkPtr (fromIntegral column) + <$> newIORef column + <*> newIORef False - columnRef <- newIORef column - return $ NamedRef (fromIntegral column) columnRef - register askVariablesRef variable - return (Variable variable) + askVariablesRef >>= register variable + setVariableName' variable (defaultVariableName variable) + pure variable -setVariableName' :: Variable Glpk -> String -> Glpk () +defaultVariableName :: GlpkVariable -> T.Text +defaultVariableName (GlpkPtr x _ _) = "x" <> T.pack (show x) + +defaultConstraintName :: GlpkConstraint -> T.Text +defaultConstraintName (GlpkPtr x _ _) = "c" <> T.pack (show x) + +setVariableName' :: GlpkVariable -> T.Text -> Glpk () setVariableName' variable name = do - problem <- askProblem - column <- readColumn variable - liftIO $ withCString name (glp_set_col_name problem column) + problem <- askProblem + column <- readColumn variable + liftIO $ withCText name (glp_set_col_name problem column) -getVariableName' :: Variable Glpk -> Glpk String +getVariableName' :: GlpkVariable -> Glpk T.Text getVariableName' variable = do problem <- askProblem column <- readColumn variable - liftIO $ glp_get_col_name problem column >>= peekCString + name <- liftIO $ glp_get_col_name problem column >>= peekCString + pure (T.pack name) -removeVariable' :: Variable Glpk -> Glpk () -removeVariable' variable = do +deleteVariable' :: GlpkVariable -> Glpk () +deleteVariable' variable = do problem <- askProblem column <- readColumn variable liftIO $ allocaGlpkArray [column] (glp_del_cols problem 1) - unregister askVariablesRef (fromVariable variable) + askVariablesRef >>= unregister variable -addConstraint' :: Inequality (LinearExpression Double (Variable Glpk)) -> Glpk (Constraint Glpk) +addConstraint' :: Inequality (Expr GlpkVariable) -> Glpk GlpkConstraint addConstraint' (Inequality ordering lhs rhs) = - let - LinearExpression terms constant = (simplify (lhs .-. rhs)) :: LinearExpression Double (Variable Glpk) + let LinExpr terms constant = simplify (lhs .-. rhs) :: Expr GlpkVariable - constraintType :: GlpkConstraintType - constraintType = case ordering of - LT -> glpkLT - GT -> glpkGT - EQ -> glpkFixed + constraintType :: GlpkConstraintType + constraintType = case ordering of + LT -> glpkLT + GT -> glpkGT + EQ -> glpkFixed - constraintRhs :: CDouble - constraintRhs = realToFrac (negate constant) + constraintRhs :: CDouble + constraintRhs = realToFrac (negate constant) - numVars :: CInt - numVars = fromIntegral (length terms) + numVars :: CInt + numVars = fromIntegral (length terms) - variables :: [Variable Glpk] - variables = map snd terms + variables :: [GlpkVariable] + variables = map snd terms - coefficients :: [CDouble] - coefficients = map (realToFrac . fst) terms - in do - problem <- askProblem - columns <- mapM readColumn variables - constraintId <- liftIO $ do - row <- glp_add_rows problem 1 - rowRef <- newIORef row - allocaGlpkArray columns $ \columnArr -> - allocaGlpkArray coefficients $ \coefficientArr -> do - glp_set_row_bnds problem row constraintType constraintRhs constraintRhs - glp_set_mat_row problem row numVars columnArr coefficientArr - return $ NamedRef (fromIntegral row) rowRef + coefficients :: [CDouble] + coefficients = map (realToFrac . fst) terms + in do + problem <- askProblem + columns <- mapM readColumn variables + constraintPtr <- liftIO $ do + row <- glp_add_rows problem 1 + allocaGlpkArray columns $ \columnArr -> + allocaGlpkArray coefficients $ \coefficientArr -> do + glp_set_row_bnds problem row constraintType constraintRhs constraintRhs + glp_set_mat_row problem row numVars columnArr coefficientArr - register askConstraintsRef constraintId - return (Constraint constraintId) + GlpkPtr (fromIntegral row) + <$> newIORef row + <*> newIORef False -setConstraintName' :: Constraint Glpk -> String -> Glpk () + askConstraintsRef >>= register constraintPtr + setConstraintName' constraintPtr (defaultConstraintName constraintPtr) + pure constraintPtr + +setConstraintName' :: GlpkConstraint -> T.Text -> Glpk () setConstraintName' constraintId name = do problem <- askProblem row <- readRow constraintId - liftIO $ withCString name (glp_set_row_name problem row) + liftIO $ withCText name (glp_set_row_name problem row) -getConstraintName' :: Constraint Glpk -> Glpk String +getConstraintName' :: GlpkConstraint -> Glpk T.Text getConstraintName' constraint = do problem <- askProblem row <- readRow constraint - liftIO $ glp_get_row_name problem row >>= peekCString + name <- liftIO $ glp_get_row_name problem row >>= peekCString + pure (T.pack name) -getDualValue' :: Constraint Glpk -> Glpk Double -getDualValue' constraint = do +getDualValue :: GlpkConstraint -> Glpk Double +getDualValue constraint = do problem <- askProblem row <- readRow constraint fmap realToFrac . liftIO $ glp_get_row_dual problem row -removeConstraint' :: Constraint Glpk -> Glpk () -removeConstraint' constraintId = do +deleteConstraint' :: GlpkConstraint -> Glpk () +deleteConstraint' constraint = do problem <- askProblem - row <- readRow constraintId + row <- readRow constraint liftIO $ allocaGlpkArray [row] (glp_del_rows problem 1) - unregister askConstraintsRef (fromConstraint constraintId) + askConstraintsRef >>= unregister constraint -addObjective' :: LinearExpression Double (Variable Glpk) -> Glpk (Objective Glpk) +addObjective' :: Expr GlpkVariable -> Glpk GlpkObjective addObjective' expr = - let - LinearExpression terms constant = simplify expr - in do - problem <- askProblem + let LinExpr terms constant = simplify expr + in do + problem <- askProblem - -- Set the constant term - liftIO $ glp_set_obj_coef problem (GlpkInt 0) (realToFrac constant) + -- Set the constant term + liftIO $ glp_set_obj_coef problem (GlpkInt 0) (realToFrac constant) - -- Set the variable terms - forM_ terms $ \(coef, variable) -> do - column <- readColumn variable - liftIO $ glp_set_obj_coef problem column (realToFrac coef) + -- Set the variable terms + forM_ terms $ \(coef, variable) -> do + column <- readColumn variable + liftIO $ glp_set_obj_coef problem column (realToFrac coef) - pure Objective + pure (GlpkObjective ()) -getObjectiveName' :: Objective Glpk -> Glpk String +-- | Delete an objective +-- +-- There is nothing to actually delete, so we just set a zero objective +deleteObjective' :: GlpkObjective -> Glpk () +deleteObjective' _ = void (addObjective' mempty) + +getObjectiveName' :: GlpkObjective -> Glpk T.Text getObjectiveName' _ = do problem <- askProblem - liftIO $ glp_get_obj_name problem >>= peekCString + name <- liftIO $ glp_get_obj_name problem >>= peekCString + pure (T.pack name) -setObjectiveName' :: Objective Glpk -> String -> Glpk () +setObjectiveName' :: GlpkObjective -> T.Text -> Glpk () setObjectiveName' _ name = do problem <- askProblem - liftIO $ withCString name (glp_set_obj_name problem) + liftIO $ withCText name (glp_set_obj_name problem) -getObjectiveSense' :: Objective Glpk -> Glpk Sense -getObjectiveSense' _ = do +getSense' :: GlpkObjective -> Glpk Sense +getSense' _ = do problem <- askProblem direction <- liftIO $ glp_get_obj_dir problem if direction == glpkMin then pure Minimization else pure Maximization -setObjectiveSense' :: Objective Glpk -> Sense -> Glpk () -setObjectiveSense' _ sense = - let - direction = case sense of - Minimization -> glpkMin - Maximization -> glpkMax - in do - problem <- askProblem - liftIO $ glp_set_obj_dir problem direction +setSense' :: GlpkObjective -> Sense -> Glpk () +setSense' _ sense = + let direction = case sense of + Minimization -> glpkMin + Maximization -> glpkMax + in do + problem <- askProblem + liftIO $ glp_set_obj_dir problem direction -getObjectiveValue' :: Objective Glpk -> Glpk Double +getObjectiveValue' :: GlpkObjective -> Glpk Double getObjectiveValue' _ = do problem <- askProblem - lastSolveRef <- asks _glpkLastSolveType + lastSolveRef <- askGlpk _glpkLastSolveType lastSolve <- (liftIO . readIORef) lastSolveRef fmap realToFrac . liftIO $ case lastSolve of - Just MIP -> glp_mip_obj_val problem - Just LP -> glp_get_obj_val problem + Just MIP -> glp_mip_obj_val problem + Just LP -> glp_get_obj_val problem Just InteriorPoint -> glp_ipt_obj_val problem - Nothing -> glp_get_obj_val problem -- There's been no solve, so who cares + Nothing -> glp_get_obj_val problem -- There's been no solve, so who cares optimizeLP' :: Glpk SolutionStatus optimizeLP' = do - -- Note that we've run an LP solve - solveTypeRef <- asks _glpkLastSolveType - liftIO $ writeIORef solveTypeRef (Just LP) + -- Note that we've run an LP solve + solveTypeRef <- askGlpk _glpkLastSolveType + liftIO $ writeIORef solveTypeRef (Just LP) - -- Run Simplex - problem <- askProblem - controlRef <- asks _glpkSimplexControl - liftIO $ do - control <- readIORef controlRef - alloca $ \controlPtr -> do - poke controlPtr control - _<- glp_simplex problem controlPtr - glp_get_status problem >>= pure . solutionStatus + -- Run Simplex + problem <- askProblem + controlRef <- askGlpk _glpkSimplexControl + liftIO $ do + control <- readIORef controlRef + alloca $ \controlPtr -> do + poke controlPtr control + _ <- glp_simplex problem controlPtr + glp_get_status problem Data.Functor.<&> solutionStatus optimizeIP' :: Glpk SolutionStatus optimizeIP' = do - -- Note that we've run a MIP solve - solveTypeRef <- asks _glpkLastSolveType - liftIO $ writeIORef solveTypeRef (Just MIP) + -- Note that we've run a MIP solve + solveTypeRef <- askGlpk _glpkLastSolveType + liftIO $ writeIORef solveTypeRef (Just MIP) - problem <- askProblem - controlRef <- asks _glpkMIPControl - liftIO $ do - control <- readIORef controlRef - alloca $ \controlPtr -> do - poke controlPtr control - _ <- glp_intopt problem controlPtr - glp_mip_status problem >>= pure . solutionStatus + problem <- askProblem + controlRef <- askGlpk _glpkMIPControl + liftIO $ do + control <- readIORef controlRef + alloca $ \controlPtr -> do + poke controlPtr control + _ <- glp_intopt problem controlPtr + glp_mip_status problem Data.Functor.<&> solutionStatus -setVariableBounds' :: Variable Glpk -> Bounds Double -> Glpk () +setVariableBounds' :: GlpkVariable -> Bounds -> Glpk () setVariableBounds' variable bounds = - let - (boundType, cLow, cHigh) = case bounds of - Free -> (glpkFree, 0, 0) - NonNegativeReals -> (glpkGT, 0, 0) - NonPositiveReals -> (glpkLT, 0, 0) - Interval low high -> (glpkBounded, realToFrac low, realToFrac high) - in do - problem <- askProblem - column <- readColumn variable - liftIO $ glp_set_col_bnds problem column boundType cLow cHigh + let (boundType, cLow, cHigh) = case bounds of + Free -> (glpkFree, 0, 0) + NonNegativeReals -> (glpkGT, 0, 0) + NonPositiveReals -> (glpkLT, 0, 0) + Interval low high -> (glpkBounded, realToFrac low, realToFrac high) + in do + problem <- askProblem + column <- readColumn variable + liftIO $ glp_set_col_bnds problem column boundType cLow cHigh -getVariableBounds' :: Variable Glpk -> Glpk (Bounds Double) +getVariableBounds' :: GlpkVariable -> Glpk Bounds getVariableBounds' variable = - let - boundsFor lb ub | lb == -maxCDouble && ub == maxCDouble = Free - | lb == -maxCDouble && ub == 0.0 = NonPositiveReals - | lb == 0.0 && ub == maxCDouble = NonNegativeReals - | otherwise = Interval lb' ub' - where - lb' = realToFrac lb - ub' = realToFrac ub - in do - problem <- askProblem - column <- readColumn variable - lb <- liftIO (glp_get_col_lb problem column) - ub <- liftIO (glp_get_col_ub problem column) - return (boundsFor lb ub) + let boundsFor lb ub + | lb == - maxCDouble && ub == maxCDouble = Free + | lb == - maxCDouble && ub == 0.0 = NonPositiveReals + | lb == 0.0 && ub == maxCDouble = NonNegativeReals + | otherwise = Interval lb' ub' + where + lb' = realToFrac lb + ub' = realToFrac ub + in do + problem <- askProblem + column <- readColumn variable + lb <- liftIO (glp_get_col_lb problem column) + ub <- liftIO (glp_get_col_ub problem column) + return (boundsFor lb ub) -setVariableDomain' :: Variable Glpk -> Domain -> Glpk () +setVariableDomain' :: GlpkVariable -> Domain -> Glpk () setVariableDomain' variable domain = - let - vType = case domain of - Continuous -> glpkContinuous - Integer -> glpkInteger - Binary -> glpkBinary - in do - problem <- askProblem - column <- readColumn variable - liftIO $ glp_set_col_kind problem column vType + let vType = case domain of + Continuous -> glpkContinuous + Integer -> glpkInteger + Binary -> glpkBinary + in do + problem <- askProblem + column <- readColumn variable + liftIO $ glp_set_col_kind problem column vType -getVariableDomain' :: Variable Glpk -> Glpk Domain +getVariableDomain' :: GlpkVariable -> Glpk Domain getVariableDomain' variable = - let - getDomain :: GlpkVariableType -> Glpk Domain - getDomain vType | vType == glpkContinuous = return Continuous - getDomain vType | vType == glpkInteger = return Integer - getDomain vType | vType == glpkBinary = return Binary - | otherwise = throwError unknownCode - where - typeName = show (typeOf vType) - GlpkVariableType code = vType - unknownCode = UnknownCode typeName code - in do - problem <- askProblem - column <- readColumn variable - getDomain =<< liftIO (glp_get_col_kind problem column) + let getDomain' :: GlpkVariableType -> Glpk Domain + getDomain' vType | vType == glpkContinuous = return Continuous + getDomain' vType | vType == glpkInteger = return Integer + getDomain' vType + | vType == glpkBinary = return Binary + | otherwise = throwIO unknownCode + where + typeName = T.pack . show . typeOf $ vType + GlpkVariableType code = vType + unknownCode = UnknownCode typeName code + in do + problem <- askProblem + column <- readColumn variable + getDomain' =<< liftIO (glp_get_col_kind problem column) -getVariableValue' :: Variable Glpk -> Glpk Double +getVariableValue' :: GlpkVariable -> Glpk Double getVariableValue' variable = do - lastSolveRef <- asks _glpkLastSolveType + lastSolveRef <- askGlpk _glpkLastSolveType lastSolve <- (liftIO . readIORef) lastSolveRef let method = case lastSolve of - Nothing -> glp_get_col_prim - Just LP -> glp_get_col_prim - Just MIP -> glp_mip_col_val + Nothing -> glp_get_col_prim + Just LP -> glp_get_col_prim + Just MIP -> glp_mip_col_val Just InteriorPoint -> glp_ipt_col_prim problem <- askProblem @@ -469,51 +546,49 @@ getTimeout' :: RealFrac a => Glpk a getTimeout' = - let - fromMillis :: RealFrac a => CInt -> a - fromMillis millis = realToFrac millis / 1000 - in do - controlRef <- asks _glpkSimplexControl - control <- liftIO (readIORef controlRef) - return $ fromMillis (smcpTimeLimitMillis control) + let fromMillis :: RealFrac a => CInt -> a + fromMillis millis = realToFrac millis / 1000 + in do + controlRef <- askGlpk _glpkSimplexControl + control <- liftIO (readIORef controlRef) + return $ fromMillis (smcpTimeLimitMillis control) setTimeout' :: RealFrac a => a -> Glpk () setTimeout' seconds = - let - millis :: Integer - millis = round (seconds * 1000) - in do - controlRef <- asks _glpkSimplexControl - control <- liftIO (readIORef controlRef) - let control' = control { smcpTimeLimitMillis = fromIntegral millis } - liftIO (writeIORef controlRef control') + let millis :: Integer + millis = round (seconds * 1000) + in do + controlRef <- askGlpk _glpkSimplexControl + control <- liftIO (readIORef controlRef) + let control' = control {smcpTimeLimitMillis = fromIntegral millis} + liftIO (writeIORef controlRef control') setRelativeMIPGap' :: RealFrac a => a -> Glpk () setRelativeMIPGap' gap = do - controlRef <- asks _glpkMIPControl + controlRef <- askGlpk _glpkMIPControl control <- liftIO (readIORef controlRef) - let control' = control { iocpRelativeMIPGap = realToFrac gap } + let control' = control {iocpRelativeMIPGap = realToFrac gap} liftIO (writeIORef controlRef control') getRelativeMIPGap' :: RealFrac a => Glpk a getRelativeMIPGap' = do - controlRef <- asks _glpkMIPControl + controlRef <- askGlpk _glpkMIPControl control <- liftIO (readIORef controlRef) return $ realToFrac (iocpRelativeMIPGap control) solutionStatus :: GlpkSolutionStatus -> SolutionStatus solutionStatus status - | status == glpkOptimal = Optimal - | status == glpkFeasible = Feasible + | status == glpkOptimal = Optimal + | status == glpkFeasible = Feasible | status == glpkInfeasible = Infeasible | status == glpkNoFeasible = Infeasible - | status == glpkUnbounded = Unbounded - | status == glpkUndefined = Infeasible - | otherwise = Error + | status == glpkUnbounded = Unbounded + | status == glpkUndefined = Infeasible + | otherwise = Error -- | Write out the current formulation to a file. -writeFormulation :: FilePath -> Glpk () -writeFormulation fileName = do +writeFormulation' :: FilePath -> Glpk () +writeFormulation' fileName = do problem <- askProblem _ <- liftIO $ withCString fileName (glp_write_lp problem nullPtr) return () @@ -526,3 +601,6 @@ (_, maxExponent) = floatRange (undefined :: CDouble) significand' = base ^ precision - 1 exponent' = maxExponent - precision + +withCText :: T.Text -> (CString -> IO a) -> IO a +withCText = withCString . T.unpack diff --git a/test/ApiTests.hs b/test/ApiTests.hs deleted file mode 100644 --- a/test/ApiTests.hs +++ /dev/null @@ -1,18 +0,0 @@ -module ApiTests where - -import Test.Tasty -import Test.Tasty.HUnit - -import Math.Programming.Glpk -import Math.Programming.Tests - -test_tree :: TestTree -test_tree = makeAllTests "GLPK" glpkRunner - -glpkRunner :: Glpk () -> IO () -glpkRunner program = do - result <- runGlpk program - case result of - Left errorMsg -> assertFailure (show errorMsg) - Right () -> return () - return () diff --git a/test/Driver.hs b/test/Driver.hs deleted file mode 100644 --- a/test/Driver.hs +++ /dev/null @@ -1,10 +0,0 @@ -{-# OPTIONS_GHC -F -pgmF tasty-discover -optF --ignores="*~" #-} -{-| -Module : Main - -The entry point for tests. We use tasty-discover to automatically -discover test functions in this directory, so simply preface tasty -tests with test_ for tasty TestTrees, or prop_ for QuickCheck -properties. --} -module Main where diff --git a/test/Math/Programming/GlpkSpec.hs b/test/Math/Programming/GlpkSpec.hs new file mode 100644 --- /dev/null +++ b/test/Math/Programming/GlpkSpec.hs @@ -0,0 +1,120 @@ +module Math.Programming.GlpkSpec where + +import Control.Monad +import Control.Monad.IO.Class +import Data.IORef +import Math.Programming +import Math.Programming.Glpk +import Math.Programming.Tests +import Math.Programming.Tests.Fuzz +import Math.Programming.Tests.IP (simpleMIPTest) +import Math.Programming.Tests.LP (dietProblemTest) +import Test.Hspec +import UnliftIO.Async + +logFormulation :: IORef Int -> Glpk () +logFormulation iRef = do + i <- liftIO $ readIORef iRef + writeFormulation ("instance" <> show i <> ".mip") + liftIO $ modifyIORef' iRef succ + +spec :: Spec +spec = do + makeAllTests "GLPK" runGlpk + + -- Toggle this to log formulations generated during fuzzing + let writeFuzzedFormulations = False + when writeFuzzedFormulations $ do + iRef <- runIO $ newIORef 0 + makeFuzzTests (\action -> runGlpk (action >> logFormulation iRef)) + + describe "Regression tests" $ do + it "solves an LP with free variables" testFreeVariablesLP + it "solves an IP with free variables" testFreeVariablesIP + it "finds an infeasible LP to be infeasible" testInfeasibleLP + it "finds an infeasible IP to be infeasible" testInfeasibleIP + + describe "Threaded runtime tests" $ do + it "can solve problems in parallel" testParallelSolves + +assertFeasible :: SolutionStatus -> Glpk () +assertFeasible result = + liftIO $ case result of + Error -> expectationFailure "Failed to solve program" + Unbounded -> expectationFailure "Unbounded program" + Infeasible -> expectationFailure "Infeasible program" + _ -> pure () + +testFreeVariablesLP :: IO () +testFreeVariablesLP = runGlpk $ do + x <- free + y <- free + z <- free + + _ <- var x .== 0 + _ <- var y .== 3.1 + _ <- var z .== -3.1 + + optimizeLP >>= assertFeasible + + vx <- getVariableValue x + vy <- getVariableValue y + vz <- getVariableValue z + + liftIO $ 0 `shouldBe` vx + liftIO $ 3.1 `shouldBe` vy + liftIO $ -3.1 `shouldBe` vz + +testFreeVariablesIP :: IO () +testFreeVariablesIP = runGlpk $ do + x <- integer + y <- integer + z <- integer + + _ <- var x .== 0 + _ <- var y .== 3 + _ <- var z .== -3 + + optimizeIP >>= assertFeasible + + vx <- getVariableValue x + vy <- getVariableValue y + vz <- getVariableValue z + + liftIO $ 0 `shouldBe` vx + liftIO $ 3 `shouldBe` vy + liftIO $ -3 `shouldBe` vz + +testInfeasibleLP :: IO () +testInfeasibleLP = runGlpk $ do + x <- free + _ <- var x .>= 2 + _ <- var x .<= 1 + + status <- optimizeLP + + liftIO $ Infeasible `shouldBe` status + +testInfeasibleIP :: IO () +testInfeasibleIP = runGlpk $ do + x <- integer + _ <- var x .>= 2 + _ <- var x .<= 1 + + status <- optimizeIP + + liftIO $ Infeasible `shouldBe` status + +testParallelSolves :: IO () +testParallelSolves = + let problems = + [ dietProblemTest, + simpleMIPTest, + dietProblemTest, + simpleMIPTest, + dietProblemTest, + simpleMIPTest, + dietProblemTest, + simpleMIPTest + ] + in pooledMapConcurrently_ runGlpk problems diff --git a/test/RegressionTests.hs b/test/RegressionTests.hs deleted file mode 100644 --- a/test/RegressionTests.hs +++ /dev/null @@ -1,91 +0,0 @@ -module RegressionTests where - -import Control.Monad.IO.Class -import Test.Tasty -import Test.Tasty.HUnit - -import Math.Programming -import Math.Programming.Glpk - -test_tree :: TestTree -test_tree = testGroup "Regression tests" - [ testCase "Free variables (LP)" testFreeVariablesLP - , testCase "Free variables (IP)" testFreeVariablesIP - , testCase "Infeasible (LP)" testInfeasibleLP - , testCase "Infeasible (IP)" testInfeasibleIP - ] - -assertFeasible :: SolutionStatus -> Glpk () -assertFeasible result - = liftIO $ case result of - Error -> assertFailure "Failed to solve program" - Unbounded -> assertFailure "Unbounded program" - Infeasible -> assertFailure "Infeasible program" - _ -> pure () - -assertRunGlpk :: Glpk a -> Assertion -assertRunGlpk program = do - eResult <- runGlpk program - case eResult of - Left err -> assertFailure ("Error solving program: " <> show err) - Right _ -> pure () - -testFreeVariablesLP :: Assertion -testFreeVariablesLP = assertRunGlpk $ do - x <- free - y <- free - z <- free - - _ <- x @==# 0 - _ <- y @==# 3.1 - _ <- z @==# -3.1 - - optimizeLP >>= assertFeasible - - vx <- getVariableValue x - vy <- getVariableValue y - vz <- getVariableValue z - - liftIO $ 0 @=? vx - liftIO $ 3.1 @=? vy - liftIO $ -3.1 @=? vz - -testFreeVariablesIP :: Assertion -testFreeVariablesIP = assertRunGlpk $ do - x <- integer - y <- integer - z <- integer - - _ <- x @==# 0 - _ <- y @==# 3 - _ <- z @==# -3 - - optimizeIP >>= assertFeasible - - vx <- getVariableValue x - vy <- getVariableValue y - vz <- getVariableValue z - - liftIO $ 0 @=? vx - liftIO $ 3 @=? vy - liftIO $ -3 @=? vz - -testInfeasibleLP :: Assertion -testInfeasibleLP = assertRunGlpk $ do - x <- free - _ <- x @>=# 2 - _ <- x @<=# 1 - - status <- optimizeLP - - liftIO $ Infeasible @=? status - -testInfeasibleIP :: Assertion -testInfeasibleIP = assertRunGlpk $ do - x <- integer - _ <- x @>=# 2 - _ <- x @<=# 1 - - status <- optimizeIP - - liftIO $ Infeasible @=? status diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -0,0 +1,3 @@ +{-# OPTIONS_GHC -F -pgmF hspec-discover #-} + +module Spec where