add MIP-glpk package

.travis.yml

@@ -11,7 +11,7 @@
 # Use new container infrastructure to enable caching
 sudo: false
 
-dist: xenial
+dist: bionic
 
 # Do not choose a language; we provide our own build tools.
 language: generic
@@ -22,8 +22,13 @@ addons:
     - libgmp-dev
     - coinor-cbc
     - glpk-utils
+    - libglpk-dev
     - lp-solve
     - c2hs
+  homebrew:
+    update: true
+    packages:
+    - glpk
 
 matrix:
   include:
@@ -66,6 +71,7 @@ install:
 - if [ -n "$HADDOCK" ]; then FLAGS="$FLAGS --haddock --no-haddock-deps"; fi
 - if [ -n "$NOZLIB" ]; then FLAGS="$FLAGS --flag MIP:-WithZLIB"; fi
 - if [ "$TRAVIS_OS_NAME" = "linux" ]; then FLAGS="$FLAGS --flag MIP:TestCBC --flag MIP:TestGlpsol --flag MIP:TestLPSolve"; fi
+- if [ "$TRAVIS_OS_NAME" = "osx" ]; then FLAGS="$FLAGS --flag MIP:TestGlpsol"; fi
 
 # Build dependencies
 - stack --jobs 2 --no-terminal --install-ghc build --test --bench --only-dependencies $FLAGS

MIP-glpk/ChangeLog.md

@@ -0,0 +1,3 @@
+# Changelog for MIP-glpk
+
+## Unreleased changes

MIP-glpk/LICENSE

@@ -0,0 +1,30 @@
+Copyright Masahiro Sakai (c) 2020
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Author name here nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

MIP-glpk/README.md

@@ -0,0 +1 @@
+# MIP-glpk

MIP-glpk/Setup.hs

@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain

MIP-glpk/package.yaml

@@ -0,0 +1,44 @@
+name:                MIP-glpk
+version:             0.1.1.0
+github:              "msakai/haskell-MIP"
+license:             BSD3
+author:              "Masahiro Sakai"
+maintainer:          "masahiro.sakai@gmail.com"
+copyright:           "2020 Masahiro Sakai"
+
+extra-source-files:
+- README.md
+- ChangeLog.md
+
+# Metadata used when publishing your package
+synopsis:            A GLPK backend to the MIP library.
+category:            Math, Algorithms, Optimisation, Optimization
+
+# To avoid duplicated efforts in documentation and dealing with the
+# complications of embedding Haddock markup inside cabal files, it is
+# common to point users to the README.md file.
+description:         Please see the README on GitHub at <https://github.com/msakai/haskell-MIP/tree/master/MIP#readme>
+
+dependencies:
+- base >= 4.7 && < 5
+- containers
+- scientific
+- extended-reals >=0.1 && <1.0
+- MIP
+- glpk-headers >=0.4.1
+
+library:
+  source-dirs: src
+  extra-libraries: glpk
+
+tests:
+  MIP-glpk-test:
+    main:                TestSuite.hs
+    source-dirs:         test
+    dependencies:
+    - MIP-glpk
+    - data-default-class
+    - tasty >=0.10.1
+    - tasty-hunit >=0.9 && <0.11
+    - tasty-quickcheck >=0.8 && <0.11
+    - tasty-th

MIP-glpk/samples/lp/infeasible.lp

@@ -0,0 +1,10 @@
+Minimize
+  0 x
+Subject To
+  3 x >= 1
+  3 x <= 2
+Bounds
+  -inf <= x <= +inf
+General
+  x
+End

MIP-glpk/samples/lp/test.lp

@@ -0,0 +1,13 @@
+\...
+Maximize
+ obj: x1 + 2 x2 + 3 x3 + x4
+Subject To
+ c1: - x1 + x2 + x3 + 10 x4 <= 20
+ c2: x1 - 3 x2 + x3 <= 30
+ c3: x2 - 3.5 x4 = 0
+Bounds
+ 0 <= x1 <= 40
+ 2 <= x4 <= 3
+General
+ x4
+End

MIP-glpk/samples/lp/unbounded-ip.lp

@@ -0,0 +1,10 @@
+Maximize
+  obj: x + y
+Subject To
+  c1: y - 2 x >= 0
+Bounds 
+  1 <= x
+  0 <= y
+General
+  x y
+End

MIP-glpk/src/Numeric/Optimization/MIP/Solver/GLPK.hs

@@ -0,0 +1,170 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# OPTIONS_HADDOCK show-extensions #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Optimization.MIP.Solver.GLPK
+-- Copyright   :  (c) Masahiro Sakai 2020
+-- License     :  BSD-style
+--
+-- Maintainer  :  masahiro.sakai@gmail.com
+-- Stability   :  provisional
+-- Portability :  non-portable
+--
+-----------------------------------------------------------------------------
+module Numeric.Optimization.MIP.Solver.GLPK
+  ( GLPK (..)
+  , glpk
+  ) where
+
+import Control.Exception
+import Control.Monad
+import qualified Data.Map.Strict as Map
+import Data.Scientific (Scientific, fromFloatDigits, toRealFloat)
+import qualified Data.Set as Set
+import Foreign
+import Foreign.C
+
+import Data.ExtendedReal
+
+import qualified Numeric.Optimization.MIP as MIP
+import Numeric.Optimization.MIP.Solver.Base
+import qualified Math.Programming.Glpk.Header as Raw
+
+data GLPK
+  = GLPK
+
+instance Default GLPK where
+  def = GLPK
+
+glpk :: GLPK
+glpk = GLPK
+
+instance IsSolver GLPK IO where
+  solve _solver opt prob = do
+    bracket Raw.glp_create_prob Raw.glp_delete_prob $ \prob' -> do
+      let vs = MIP.variables  prob
+          varToCol = Map.fromList $ zip (Set.toAscList vs) [1..]
+          exprToMap (MIP.Expr ts) = Map.fromListWith (+) $ do
+            t <- ts
+            case t of
+              MIP.Term c [] -> return (0, c)
+              MIP.Term c [v] -> return (varToCol Map.! v, c)
+              MIP.Term _ _ -> error "GLPK does not support non-linear term"
+
+      -- Variables
+      _ <- Raw.glp_add_cols prob' $ fromIntegral $ Map.size $ MIP.varType prob
+      forM_ (Map.toList varToCol) $ \(v, col) -> do
+        let (lb, ub) = MIP.getBounds prob v
+        Raw.glp_set_col_kind prob' col $
+          case MIP.getVarType prob v of
+            MIP.SemiContinuousVariable -> error "GLPK does not support semi-continuous variables"
+            MIP.SemiIntegerVariable -> error "GLPK does not support semi-integer variables"
+            MIP.ContinuousVariable -> Raw.glpkContinuous
+            MIP.IntegerVariable ->
+              case (lb, ub) of
+                (Finite 0, Finite 1) -> Raw.glpkBinary
+                _ -> Raw.glpkInteger
+        case fromBound lb ub of
+          (constrType, lb', ub') -> Raw.glp_set_col_bnds prob' col constrType lb' ub'
+
+      -- Objective Function
+      let obj = MIP.objectiveFunction prob
+      Raw.glp_set_obj_dir prob' $
+        case MIP.objDir obj of
+          MIP.OptMax -> Raw.glpkMax
+          MIP.OptMin -> Raw.glpkMin
+      forM_ (Map.toList (exprToMap (MIP.objExpr obj))) $ \(col, c) -> do
+        Raw.glp_set_obj_coef prob' col (toRealFloat c)
+
+      -- Constraints
+      let constrs = MIP.constraints prob
+      _ <- Raw.glp_add_rows prob' $ fromIntegral $ length constrs
+      forM_ (zip [1..] constrs) $ \(row, constr) -> do
+        case MIP.constrIndicator constr of
+          Nothing -> return ()
+          Just _ -> error "Indicator constraints are not supported"
+        when (MIP.constrIsLazy constr) $ do
+          error "GLPK does not support lazy constraints"
+        case fromBound (MIP.constrLB constr) (MIP.constrUB constr) of
+          (constrType, lb', ub') -> Raw.glp_set_row_bnds prob' row constrType lb' ub'
+        -- TODO: check constant terms
+        let m = exprToMap (MIP.constrExpr constr)
+            ts = Map.toList m
+            n = Map.size m
+        Raw.allocaGlpkArray (map fst ts) $ \ind -> do
+          Raw.allocaGlpkArray (map (toRealFloat . snd) ts) $ \val -> do
+            Raw.glp_set_mat_row prob' row (fromIntegral n) ind val
+
+      when (length (MIP.sosConstraints prob) > 0) $ do
+        error "GLPK does not support SOS constraints"
+
+      when (length (MIP.userCuts prob) > 0) $ do
+        error "GLPK does not support user cuts"
+
+      let loggingCallback :: CString -> IO Int
+          loggingCallback p = do
+            s <- peekCString p
+            solveLogger opt s
+            return 1
+
+      -- Solving
+      alloca $ \p -> do
+        Raw.glp_init_iocp p
+        iocp <- peek p
+        poke p $
+          iocp
+          { Raw.iocpPresolve = Raw.glpkPresolveOn
+          , Raw.iocpTimeLimitMillis =
+              case solveTimeLimit opt of
+                Nothing -> Raw.iocpTimeLimitMillis iocp -- maxBound :: CInt
+                Just sec -> round (sec * 1000)
+          }
+
+        status <-
+          bracket (newStablePtr loggingCallback) freeStablePtr $ \loggingCallbackPtr ->
+            bracket_ (Raw.glp_term_hook termHookFunPtr (castStablePtrToPtr loggingCallbackPtr)) (Raw.glp_term_hook nullFunPtr nullPtr) $
+              Raw.glp_intopt prob' p
+
+        objVal <- liftM fromFloatDigits $ Raw.glp_mip_obj_val prob'
+        varVal <- mapM (liftM fromFloatDigits . Raw.glp_mip_col_val prob') varToCol
+        let sol = MIP.Solution
+              { MIP.solStatus =
+                  if status == Raw.glpkMIPSuccess then MIP.StatusOptimal -- ???
+                  else if status == Raw.glpkMIPBadBound then MIP.StatusInfeasible
+                  else if status == Raw.glpkMIPNoBasis then MIP.StatusUnknown
+                  else if status == Raw.glpkMIPPrimalInfeasible then MIP.StatusInfeasible
+                  else if status == Raw.glpkMIPDualInfeasible then MIP.StatusInfeasibleOrUnbounded
+                  else if status == Raw.glpkMIPFailure then MIP.StatusUnknown
+                  else if status == Raw.glpkMIPRelativeGap then MIP.StatusUnknown -- ???
+                  else if status == Raw.glpkMIPTimeLimit then MIP.StatusUnknown -- ???
+                  else if status == Raw.glpkMIPStopped then MIP.StatusUnknown -- ???
+                  else error ("unknown mip status: " ++ show status)
+              , MIP.solObjectiveValue = Just objVal
+              , MIP.solVariables = varVal
+              }
+        return sol
+
+fromBound :: MIP.BoundExpr Scientific -> MIP.BoundExpr Scientific -> (Raw.GlpkConstraintType, CDouble, CDouble)
+fromBound NegInf PosInf = (Raw.glpkFree, 0, 0)
+fromBound (Finite lb') (Finite ub')
+  | lb' == ub' = (Raw.glpkFixed, toRealFloat lb', toRealFloat ub')
+  | otherwise  = (Raw.glpkBounded, toRealFloat lb', toRealFloat ub')
+fromBound (Finite lb') PosInf = (Raw.glpkGT, toRealFloat lb', 0)
+fromBound NegInf (Finite ub') = (Raw.glpkLT, 0, toRealFloat ub')
+fromBound _ NegInf = (Raw.glpkBounded, 1, 0)  -- inconsistent
+fromBound PosInf _ = (Raw.glpkBounded, 1, 0)  -- inconsistent
+
+
+termHook :: Ptr () -> CString -> IO CInt
+termHook p s = do
+  let sp :: StablePtr (CString -> IO Int)
+      sp = castPtrToStablePtr p
+  callback <- deRefStablePtr sp
+  liftM fromIntegral $ callback s
+
+foreign export ccall "haskell_mip_glpk_term_hook"
+  termHook :: Ptr () -> CString -> IO CInt
+
+foreign import ccall "&haskell_mip_glpk_term_hook"
+  termHookFunPtr :: FunPtr (Ptr () -> CString -> IO CInt)

MIP-glpk/test/TestSuite.hs

@@ -0,0 +1,53 @@
+-- {-# OPTIONS_GHC -Wall -Wno-unused-top-binds #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE TemplateHaskell #-}
+module Main where
+
+import Control.Monad
+import Data.Default.Class
+import qualified Data.Map as Map
+import Test.Tasty
+import Test.Tasty.HUnit
+import qualified Numeric.Optimization.MIP as MIP
+import Numeric.Optimization.MIP.Solver
+import Numeric.Optimization.MIP.Solver.GLPK
+
+-- ------------------------------------------------------------------------
+
+case_glpk :: Assertion
+case_glpk = do
+  prob <- MIP.readFile def "samples/lp/test.lp"
+  sol <- solve glpk def prob
+  sol @?=
+    MIP.Solution
+    { MIP.solStatus = MIP.StatusOptimal
+    , MIP.solObjectiveValue = Just 122.5
+    , MIP.solVariables = Map.fromList [("x1", 40), ("x2", 10.5), ("x3", 19.5), ("x4", 3)]
+    }
+
+case_glpk_unbounded :: Assertion
+case_glpk_unbounded = do
+  prob <- MIP.readFile def "samples/lp/unbounded-ip.lp"
+  sol <- solve glpk def prob
+  let status = MIP.solStatus sol
+  unless (status == MIP.StatusUnbounded || status == MIP.StatusFeasible || status == MIP.StatusInfeasibleOrUnbounded) $
+    assertFailure $ unlines $
+      [ "expected: StatusUnbounded, StatusFeasible or StatusInfeasibleOrUnbounded"
+      , " but got: " ++ show status
+      ]
+
+case_glpk_infeasible :: Assertion
+case_glpk_infeasible = do
+  prob <- MIP.readFile def "samples/lp/infeasible.lp"
+  sol <- solve glpk def prob
+  MIP.solStatus sol @?= MIP.StatusInfeasible
+
+-- ------------------------------------------------------------------------
+
+main :: IO ()
+main = defaultMain $ testGroup "MIP-glpk test suite"
+  [ testCase "glpk" case_glpk
+  , testCase "glpk unbounded" case_glpk_unbounded
+  , testCase "glpk infeasible" case_glpk_infeasible
+  ]

appveyor.yml

@@ -35,8 +35,12 @@ install:
 
 - curl -owinglpk-4.65.zip -L --insecure https://jaist.dl.sourceforge.net/project/winglpk/winglpk/GLPK-4.65/winglpk-4.65.zip
 - 7z x winglpk-4.65.zip
-- if %PLATFORM%==x86 set PATH=%APPVEYOR_BUILD_FOLDER%\glpk-4.65\w32;%PATH%
-- if %PLATFORM%==x64 set PATH=%APPVEYOR_BUILD_FOLDER%\glpk-4.65\w64;%PATH%
+- if %PLATFORM%==x86 set GLPK_PLATFORM=w32
+- if %PLATFORM%==x64 set GLPK_PLATFORM=w64
+- set PATH=%APPVEYOR_BUILD_FOLDER%\glpk-4.65\%GLPK_PLATFORM%;%PATH%
+- set GLPK_LIB_DIR=%APPVEYOR_BUILD_FOLDER%\glpk-4.65\%GLPK_PLATFORM%
+- set GLPK_INCLUDE_DIR=%APPVEYOR_BUILD_FOLDER%\glpk-4.65\src
+- copy glpk-4.65\%GLPK_PLATFORM%\glpk_4_65.lib glpk-4.65\%GLPK_PLATFORM%\glpk.lib
 
 - if %PLATFORM%==x86 (
     curl -olp_solve_5.5.2.5_exe_win32.zip -L --insecure https://jaist.dl.sourceforge.net/project/lpsolve/lpsolve/5.5.2.5/lp_solve_5.5.2.5_exe_win32.zip &&
@@ -54,4 +58,4 @@ build_script:
 test_script:
 # The ugly echo "" hack is to avoid complaints about 0 being an invalid file
 # descriptor
-- echo "" | stack --no-terminal test --bench --no-run-benchmarks --flag MIP:TestCBC --flag MIP:TestGlpsol --flag MIP:TestLPSolve
+- echo "" | stack --no-terminal test --bench --no-run-benchmarks --flag MIP:TestCBC --flag MIP:TestGlpsol --flag MIP:TestLPSolve --extra-include-dirs=%GLPK_INCLUDE_DIR% --extra-lib-dirs=%GLPK_LIB_DIR%