diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,8 @@ +# Changelog for MIP + +## Unreleased changes + +## 0.1.0.0 + +* Separated from toysolver package + https://github.com/msakai/toysolver/tree/7096038ece0a5f860a951567689ef9a03ac0355d diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/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. diff --git a/MIP.cabal b/MIP.cabal new file mode 100644 --- /dev/null +++ b/MIP.cabal @@ -0,0 +1,168 @@ +cabal-version: 1.18 + +-- This file has been generated from package.yaml by hpack version 0.33.0. +-- +-- see: https://github.com/sol/hpack +-- +-- hash: ac803e98baaa0544b947d71e759475b741c54bd5c39105056cfd3d1fbc25b4e0 + +name: MIP +version: 0.1.0.0 +synopsis: Library for using Mixed Integer Programming (MIP) +description: Please see the README on GitHub at +category: Math, Algorithms, Optimisation, Optimization +homepage: https://github.com/msakai/haskell-MIP#readme +bug-reports: https://github.com/msakai/haskell-MIP/issues +author: Masahiro Sakai +maintainer: masahiro.sakai@gmail.com +copyright: 2020 Masahiro Sakai +license: BSD3 +license-file: LICENSE +build-type: Simple +extra-source-files: + README.md + ChangeLog.md + doc-images/MIP-Status-diagram.tex +extra-doc-files: + doc-images/MIP-Status-diagram.png + +source-repository head + type: git + location: https://github.com/msakai/haskell-MIP + +flag TestCBC + description: run test cases that depends on cbc command + manual: True + default: False + +flag TestCPLEX + description: run test cases that depends on cplex command + manual: True + default: False + +flag TestGlpsol + description: run test cases that depends on glpsol command + manual: True + default: False + +flag TestGurobiCl + description: run test cases that depends on gurobi_cl command + manual: True + default: False + +flag TestLPSolve + description: run test cases that depends on lp_solve command + manual: True + default: False + +flag TestSCIP + description: run test cases that depends on scip command + manual: True + default: False + +flag WithZlib + description: Use zlib package to support gzipped files + manual: True + default: True + +library + exposed-modules: + Numeric.Optimization.MIP + Numeric.Optimization.MIP.Base + Numeric.Optimization.MIP.FileUtils + Numeric.Optimization.MIP.LPFile + Numeric.Optimization.MIP.MPSFile + Numeric.Optimization.MIP.Solution.CBC + Numeric.Optimization.MIP.Solution.CPLEX + Numeric.Optimization.MIP.Solution.GLPK + Numeric.Optimization.MIP.Solution.Gurobi + Numeric.Optimization.MIP.Solution.SCIP + Numeric.Optimization.MIP.Solver + Numeric.Optimization.MIP.Solver.Base + Numeric.Optimization.MIP.Solver.CBC + Numeric.Optimization.MIP.Solver.CPLEX + Numeric.Optimization.MIP.Solver.Glpsol + Numeric.Optimization.MIP.Solver.GurobiCl + Numeric.Optimization.MIP.Solver.LPSolve + Numeric.Optimization.MIP.Solver.SCIP + other-modules: + Paths_MIP + Numeric.Optimization.MIP.Internal.ProcessUtil + Numeric.Optimization.MIP.Internal.Util + hs-source-dirs: + src + other-extensions: CPP + build-depends: + OptDir + , base >=4.7 && <5 + , containers >=0.5.0 + , data-default-class + , extended-reals >=0.1 && <1.0 + , filepath + , intern >=0.9.1.2 && <1.0.0.0 + , lattices + , megaparsec >=5 && <9 + , mtl >=2.1.2 + , process >=1.1.0.2 + , scientific + , stm >=2.3 + , temporary >=1.2 + , text >=1.1.0.0 + , xml-conduit + if flag(WithZlib) + cpp-options: -DWITH_ZLIB + build-depends: + bytestring + , bytestring-encoding + , case-insensitive + , zlib + else + cpp-options: + default-language: Haskell2010 + +test-suite MIP-test + type: exitcode-stdio-1.0 + main-is: TestSuite.hs + other-modules: + Test.LPFile + Test.MIP + Test.MIPSolver + Test.MPSFile + Paths_MIP + hs-source-dirs: + test + build-depends: + MIP + , base >=4.7 && <5 + , containers >=0.5.0 + , data-default-class + , lattices + , tasty >=0.10.1 + , tasty-hunit >=0.9 && <0.11 + , tasty-quickcheck >=0.8 && <0.11 + , tasty-th + if flag(TestCBC) + cpp-options: -DTEST_CBC + else + cpp-options: + if flag(TestCPLEX) + cpp-options: -DTEST_CPLEX + else + cpp-options: + if flag(TestGlpsol) + cpp-options: -DTEST_GLPSOL + else + cpp-options: + if flag(TestGurobiCl) + cpp-options: -DTEST_GUROBI_CL + else + cpp-options: + if flag(TestLPSolve) + cpp-options: -DTEST_LP_SOLVE + else + cpp-options: + if flag(TestSCIP) + cpp-options: -DTEST_SCIP + else + cpp-options: + default-language: Haskell2010 diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,8 @@ +# MIP + +Library for using Mixed Integer Programming (MIP) in Haskell. +This library contains functions like: + +* Reading / Writing MIP problem files (e.g. `LP` file or `MPS` file), +* Invokling MIP solvers like Gurobi, CPLEX, CBC, GLPK, lp_solve, +* Reading solution files of those solvers. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/doc-images/MIP-Status-diagram.png b/doc-images/MIP-Status-diagram.png new file mode 100644 Binary files /dev/null and b/doc-images/MIP-Status-diagram.png differ diff --git a/doc-images/MIP-Status-diagram.tex b/doc-images/MIP-Status-diagram.tex new file mode 100644 --- /dev/null +++ b/doc-images/MIP-Status-diagram.tex @@ -0,0 +1,18 @@ +\documentclass{article} +\usepackage{tikz} +\begin{document} + +\begin{tikzpicture}[scale=.7] + \node (Optimal) at (-3,2) {Optimal}; + \node (Feasible) at (-3,0) {Feasible}; + \node (InfeasibleOrUnbounded) at (2,0) {InfeasibleOrUnbounded}; + \node (Unbounded) at (2,2) {Unbounded}; + \node (Infeasible) at (6,2) {Infeasible}; + \node (Unknown) at (0,-2) {Unknown}; + \draw (Unknown) -- (Feasible) -- (Optimal); + \draw (Feasible) -- (Unbounded); + \draw (Unknown) -- (InfeasibleOrUnbounded) -- (Unbounded); + \draw (InfeasibleOrUnbounded) -- (Infeasible); +\end{tikzpicture} + +\end{document} diff --git a/src/Numeric/Optimization/MIP.hs b/src/Numeric/Optimization/MIP.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP.hs @@ -0,0 +1,166 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP +-- Copyright : (c) Masahiro Sakai 2011-2014 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +-- Mixed-Integer Programming Problems with some commmonly used extensions +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP + ( module Numeric.Optimization.MIP.Base + , readFile + , readLPFile + , readMPSFile + , parseLPString + , parseMPSString + , writeFile + , writeLPFile + , writeMPSFile + , toLPString + , toMPSString + , ParseError + ) where + +import Prelude hiding (readFile, writeFile) +import Control.Exception +import Data.Char +import Data.Scientific (Scientific) +import qualified Data.Text.Lazy as TL +import qualified Data.Text.Lazy.IO as TLIO +import System.FilePath (takeExtension, splitExtension) +import System.IO hiding (readFile, writeFile) + +import Numeric.Optimization.MIP.Base +import Numeric.Optimization.MIP.FileUtils (ParseError) +import qualified Numeric.Optimization.MIP.LPFile as LPFile +import qualified Numeric.Optimization.MIP.MPSFile as MPSFile + +#ifdef WITH_ZLIB +import qualified Codec.Compression.GZip as GZip +import qualified Data.ByteString.Lazy as BL +import Data.ByteString.Lazy.Encoding (encode, decode) +import qualified Data.CaseInsensitive as CI +import GHC.IO.Encoding (getLocaleEncoding) +#endif + +-- | Parse .lp or .mps file based on file extension +readFile :: FileOptions -> FilePath -> IO (Problem Scientific) +readFile opt fname = + case getExt fname of + ".lp" -> readLPFile opt fname + ".mps" -> readMPSFile opt fname + ext -> ioError $ userError $ "unknown extension: " ++ ext + +-- | Parse a file containing LP file data. +readLPFile :: FileOptions -> FilePath -> IO (Problem Scientific) +#ifndef WITH_ZLIB +readLPFile = LPFile.parseFile +#else +readLPFile opt fname = do + s <- readTextFile opt fname + let ret = LPFile.parseString opt fname s + case ret of + Left e -> throw e + Right a -> return a +#endif + +-- | Parse a file containing MPS file data. +readMPSFile :: FileOptions -> FilePath -> IO (Problem Scientific) +#ifndef WITH_ZLIB +readMPSFile = MPSFile.parseFile +#else +readMPSFile opt fname = do + s <- readTextFile opt fname + let ret = MPSFile.parseString opt fname s + case ret of + Left e -> throw e + Right a -> return a +#endif + +readTextFile :: FileOptions -> FilePath -> IO TL.Text +#ifndef WITH_ZLIB +readTextFile opt fname = do + h <- openFile fname ReadMode + case optFileEncoding opt of + Nothing -> return () + Just enc -> hSetEncoding h enc + TLIO.hGetContents h +#else +readTextFile opt fname = do + enc <- case optFileEncoding opt of + Nothing -> getLocaleEncoding + Just enc -> return enc + let f = if CI.mk (takeExtension fname) == ".gz" then GZip.decompress else id + s <- BL.readFile fname + return $ decode enc $ f s +#endif + +-- | Parse a string containing LP file data. +parseLPString :: FileOptions -> String -> String -> Either (ParseError String) (Problem Scientific) +parseLPString = LPFile.parseString + +-- | Parse a string containing MPS file data. +parseMPSString :: FileOptions -> String -> String -> Either (ParseError String) (Problem Scientific) +parseMPSString = MPSFile.parseString + +writeFile :: FileOptions -> FilePath -> Problem Scientific -> IO () +writeFile opt fname prob = + case getExt fname of + ".lp" -> writeLPFile opt fname prob + ".mps" -> writeMPSFile opt fname prob + ext -> ioError $ userError $ "unknown extension: " ++ ext + +getExt :: String -> String +getExt fname | (base, ext) <- splitExtension fname = + case map toLower ext of +#ifdef WITH_ZLIB + ".gz" -> getExt base +#endif + s -> s + +writeLPFile :: FileOptions -> FilePath -> Problem Scientific -> IO () +writeLPFile opt fname prob = + case LPFile.render opt prob of + Left err -> ioError $ userError err + Right s -> writeTextFile opt fname s + +writeMPSFile :: FileOptions -> FilePath -> Problem Scientific -> IO () +writeMPSFile opt fname prob = + case MPSFile.render opt prob of + Left err -> ioError $ userError err + Right s -> writeTextFile opt fname s + +writeTextFile :: FileOptions -> FilePath -> TL.Text -> IO () +writeTextFile opt fname s = do + let writeSimple = do + withFile fname WriteMode $ \h -> do + case optFileEncoding opt of + Nothing -> return () + Just enc -> hSetEncoding h enc + TLIO.hPutStr h s +#ifdef WITH_ZLIB + if CI.mk (takeExtension fname) /= ".gz" then do + writeSimple + else do + enc <- case optFileEncoding opt of + Nothing -> getLocaleEncoding + Just enc -> return enc + BL.writeFile fname $ GZip.compress $ encode enc s +#else + writeSimple +#endif + +toLPString :: FileOptions -> Problem Scientific -> Either String TL.Text +toLPString = LPFile.render + +toMPSString :: FileOptions -> Problem Scientific -> Either String TL.Text +toMPSString = MPSFile.render diff --git a/src/Numeric/Optimization/MIP/Base.hs b/src/Numeric/Optimization/MIP/Base.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Base.hs @@ -0,0 +1,465 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE CPP #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Base +-- Copyright : (c) Masahiro Sakai 2011-2019 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +-- Mixed-Integer Programming Problems with some commmonly used extensions +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Base + ( + -- * The MIP Problem type + Problem (..) + , Label + + -- * Variables + , Var + , toVar + , fromVar + + -- ** Variable types + , VarType (..) + , getVarType + + -- ** Variable bounds + , BoundExpr + , Extended (..) + , Bounds + , defaultBounds + , defaultLB + , defaultUB + , getBounds + + -- ** Variable getters + , variables + , integerVariables + , semiContinuousVariables + , semiIntegerVariables + + -- * Expressions + , Expr (..) + , varExpr + , constExpr + , terms + , Term (..) + + -- * Objective function + , OptDir (..) + , ObjectiveFunction (..) + + -- * Constraints + + -- ** Linear (or Quadratic or Polynomial) constraints + , Constraint (..) + , (.==.) + , (.<=.) + , (.>=.) + , RelOp (..) + + -- ** SOS constraints + , SOSType (..) + , SOSConstraint (..) + + -- * Solutions + , Solution (..) + , Status (..) + , meetStatus + + -- * File I/O options + , FileOptions (..) + + -- * Utilities + , Variables (..) + , intersectBounds + ) where + +#if !MIN_VERSION_lattices(2,0,0) +import Algebra.Lattice +#endif +import Algebra.PartialOrd +import Control.Arrow ((***)) +import Data.Default.Class +import Data.Map (Map) +import qualified Data.Map as Map +import Data.Set (Set) +import qualified Data.Set as Set +import Data.Interned (unintern) +import Data.Interned.Text +import Data.ExtendedReal +import Data.OptDir +import Data.String +import qualified Data.Text as T +import System.IO (TextEncoding) + +infix 4 .<=., .>=., .==. + +-- --------------------------------------------------------------------------- + +-- | Problem +data Problem c + = Problem + { name :: Maybe T.Text + , objectiveFunction :: ObjectiveFunction c + , constraints :: [Constraint c] + , sosConstraints :: [SOSConstraint c] + , userCuts :: [Constraint c] + , varType :: Map Var VarType + , varBounds :: Map Var (Bounds c) + } + deriving (Show, Eq, Ord) + +instance Default (Problem c) where + def = Problem + { name = Nothing + , objectiveFunction = def + , constraints = [] + , sosConstraints = [] + , userCuts = [] + , varType = Map.empty + , varBounds = Map.empty + } + +instance Functor Problem where + fmap f prob = + prob + { objectiveFunction = fmap f (objectiveFunction prob) + , constraints = map (fmap f) (constraints prob) + , sosConstraints = map (fmap f) (sosConstraints prob) + , userCuts = map (fmap f) (userCuts prob) + , varBounds = fmap (fmap f *** fmap f) (varBounds prob) + } + +-- | label +type Label = T.Text + +-- --------------------------------------------------------------------------- + +-- | variable +type Var = InternedText + +-- | convert a string into a variable +toVar :: String -> Var +toVar = fromString + +-- | convert a variable into a string +fromVar :: Var -> String +fromVar = T.unpack . unintern + +data VarType + = ContinuousVariable + | IntegerVariable + | SemiContinuousVariable + | SemiIntegerVariable + deriving (Eq, Ord, Show) + +instance Default VarType where + def = ContinuousVariable + +-- | looking up bounds for a variable +getVarType :: Problem c -> Var -> VarType +getVarType mip v = Map.findWithDefault def v (varType mip) + +-- | type for representing lower/upper bound of variables +type BoundExpr c = Extended c + +-- | type for representing lower/upper bound of variables +type Bounds c = (BoundExpr c, BoundExpr c) + +-- | default bounds +defaultBounds :: Num c => Bounds c +defaultBounds = (defaultLB, defaultUB) + +-- | default lower bound (0) +defaultLB :: Num c => BoundExpr c +defaultLB = Finite 0 + +-- | default upper bound (+∞) +defaultUB :: BoundExpr c +defaultUB = PosInf + +-- | looking up bounds for a variable +getBounds :: Num c => Problem c -> Var -> Bounds c +getBounds mip v = Map.findWithDefault defaultBounds v (varBounds mip) + +intersectBounds :: Ord c => Bounds c -> Bounds c -> Bounds c +intersectBounds (lb1,ub1) (lb2,ub2) = (max lb1 lb2, min ub1 ub2) + +-- --------------------------------------------------------------------------- + +-- | expressions +newtype Expr c = Expr [Term c] + deriving (Eq, Ord, Show) + +varExpr :: Num c => Var -> Expr c +varExpr v = Expr [Term 1 [v]] + +constExpr :: (Eq c, Num c) => c -> Expr c +constExpr 0 = Expr [] +constExpr c = Expr [Term c []] + +terms :: Expr c -> [Term c] +terms (Expr ts) = ts + +instance Num c => Num (Expr c) where + Expr e1 + Expr e2 = Expr (e1 ++ e2) + Expr e1 * Expr e2 = Expr [Term (c1*c2) (vs1 ++ vs2) | Term c1 vs1 <- e1, Term c2 vs2 <- e2] + negate (Expr e) = Expr [Term (-c) vs | Term c vs <- e] + abs = id + signum _ = 1 + fromInteger 0 = Expr [] + fromInteger c = Expr [Term (fromInteger c) []] + +instance Functor Expr where + fmap f (Expr ts) = Expr $ map (fmap f) ts + +splitConst :: Num c => Expr c -> (Expr c, c) +splitConst e = (e2, c2) + where + e2 = Expr [t | t@(Term _ (_:_)) <- terms e] + c2 = sum [c | Term c [] <- terms e] + +-- | terms +data Term c = Term c [Var] + deriving (Eq, Ord, Show) + +instance Functor Term where + fmap f (Term c vs) = Term (f c) vs + +-- --------------------------------------------------------------------------- + +-- | objective function +data ObjectiveFunction c + = ObjectiveFunction + { objLabel :: Maybe Label + , objDir :: OptDir + , objExpr :: Expr c + } + deriving (Eq, Ord, Show) + +instance Default (ObjectiveFunction c) where + def = + ObjectiveFunction + { objLabel = Nothing + , objDir = OptMin + , objExpr = Expr [] + } + +instance Functor ObjectiveFunction where + fmap f obj = obj{ objExpr = fmap f (objExpr obj) } + +-- --------------------------------------------------------------------------- + +-- | constraint +data Constraint c + = Constraint + { constrLabel :: Maybe Label + , constrIndicator :: Maybe (Var, c) + , constrExpr :: Expr c + , constrLB :: BoundExpr c + , constrUB :: BoundExpr c + , constrIsLazy :: Bool + } + deriving (Eq, Ord, Show) + +(.==.) :: Num c => Expr c -> Expr c -> Constraint c +lhs .==. rhs = + case splitConst (lhs - rhs) of + (e, c) -> def{ constrExpr = e, constrLB = Finite (- c), constrUB = Finite (- c) } + +(.<=.) :: Num c => Expr c -> Expr c -> Constraint c +lhs .<=. rhs = + case splitConst (lhs - rhs) of + (e, c) -> def{ constrExpr = e, constrUB = Finite (- c) } + +(.>=.) :: Num c => Expr c -> Expr c -> Constraint c +lhs .>=. rhs = + case splitConst (lhs - rhs) of + (e, c) -> def{ constrExpr = e, constrLB = Finite (- c) } + +instance Default (Constraint c) where + def = Constraint + { constrLabel = Nothing + , constrIndicator = Nothing + , constrExpr = Expr [] + , constrLB = NegInf + , constrUB = PosInf + , constrIsLazy = False + } + +instance Functor Constraint where + fmap f c = + c + { constrIndicator = fmap (id *** f) (constrIndicator c) + , constrExpr = fmap f (constrExpr c) + , constrLB = fmap f (constrLB c) + , constrUB = fmap f (constrUB c) + } + +-- | relational operators +data RelOp = Le | Ge | Eql + deriving (Eq, Ord, Enum, Show) + +-- --------------------------------------------------------------------------- + +-- | types of SOS (special ordered sets) constraints +data SOSType + = S1 -- ^ Type 1 SOS constraint + | S2 -- ^ Type 2 SOS constraint + deriving (Eq, Ord, Enum, Show, Read) + +-- | SOS (special ordered sets) constraints +data SOSConstraint c + = SOSConstraint + { sosLabel :: Maybe Label + , sosType :: SOSType + , sosBody :: [(Var, c)] + } + deriving (Eq, Ord, Show) + +instance Functor SOSConstraint where + fmap f c = c{ sosBody = map (id *** f) (sosBody c) } + +-- --------------------------------------------------------------------------- + +-- | MIP status with the following partial order: +-- +-- <> +data Status + = StatusUnknown + | StatusFeasible + | StatusOptimal + | StatusInfeasibleOrUnbounded + | StatusInfeasible + | StatusUnbounded + deriving (Eq, Ord, Enum, Bounded, Show) + +instance PartialOrd Status where + leq a b = (a,b) `Set.member` rel + where + rel = unsafeLfpFrom rel0 $ \r -> + Set.union r (Set.fromList [(x,z) | (x,y) <- Set.toList r, (y',z) <- Set.toList r, y == y']) + rel0 = Set.fromList $ + [(x,x) | x <- [minBound .. maxBound]] ++ + [ (StatusUnknown, StatusFeasible) + , (StatusUnknown, StatusInfeasibleOrUnbounded) + , (StatusFeasible, StatusOptimal) + , (StatusFeasible, StatusUnbounded) + , (StatusInfeasibleOrUnbounded, StatusUnbounded) + , (StatusInfeasibleOrUnbounded, StatusInfeasible) + ] + + +meetStatus :: Status -> Status -> Status +StatusUnknown `meetStatus` _b = StatusUnknown +StatusFeasible `meetStatus` b + | StatusFeasible `leq` b = StatusFeasible + | otherwise = StatusUnknown +StatusOptimal `meetStatus` StatusOptimal = StatusOptimal +StatusOptimal `meetStatus` b + | StatusFeasible `leq` b = StatusFeasible + | otherwise = StatusUnknown +StatusInfeasibleOrUnbounded `meetStatus` b + | StatusInfeasibleOrUnbounded `leq` b = StatusInfeasibleOrUnbounded + | otherwise = StatusUnknown +StatusInfeasible `meetStatus` StatusInfeasible = StatusInfeasible +StatusInfeasible `meetStatus` b + | StatusInfeasibleOrUnbounded `leq` b = StatusInfeasibleOrUnbounded + | otherwise = StatusUnknown +StatusUnbounded `meetStatus` StatusUnbounded = StatusUnbounded +StatusUnbounded `meetStatus` b + | StatusFeasible `leq` b = StatusFeasible + | StatusInfeasibleOrUnbounded `leq` b = StatusInfeasibleOrUnbounded + | otherwise = StatusUnknown + +#if !MIN_VERSION_lattices(2,0,0) + +instance MeetSemiLattice Status where + meet = meetStatus + +#endif + + +data Solution r + = Solution + { solStatus :: Status + , solObjectiveValue :: Maybe r + , solVariables :: Map Var r + } + deriving (Eq, Ord, Show) + +instance Functor Solution where + fmap f (Solution status obj vs) = Solution status (fmap f obj) (fmap f vs) + +instance Default (Solution r) where + def = Solution + { solStatus = StatusUnknown + , solObjectiveValue = Nothing + , solVariables = Map.empty + } + +-- --------------------------------------------------------------------------- + +class Variables a where + vars :: a -> Set Var + +instance Variables a => Variables [a] where + vars = Set.unions . map vars + +instance (Variables a, Variables b) => Variables (Either a b) where + vars (Left a) = vars a + vars (Right b) = vars b + +instance Variables (Problem c) where + vars = variables + +instance Variables (Expr c) where + vars (Expr e) = vars e + +instance Variables (Term c) where + vars (Term _ xs) = Set.fromList xs + +instance Variables (ObjectiveFunction c) where + vars ObjectiveFunction{ objExpr = e } = vars e + +instance Variables (Constraint c) where + vars Constraint{ constrIndicator = ind, constrExpr = e } = Set.union (vars e) vs2 + where + vs2 = maybe Set.empty (Set.singleton . fst) ind + +instance Variables (SOSConstraint c) where + vars SOSConstraint{ sosBody = xs } = Set.fromList (map fst xs) + +-- --------------------------------------------------------------------------- + +variables :: Problem c -> Set Var +variables mip = Map.keysSet $ varType mip + +integerVariables :: Problem c -> Set Var +integerVariables mip = Map.keysSet $ Map.filter (IntegerVariable ==) (varType mip) + +semiContinuousVariables :: Problem c -> Set Var +semiContinuousVariables mip = Map.keysSet $ Map.filter (SemiContinuousVariable ==) (varType mip) + +semiIntegerVariables :: Problem c -> Set Var +semiIntegerVariables mip = Map.keysSet $ Map.filter (SemiIntegerVariable ==) (varType mip) + +-- --------------------------------------------------------------------------- + +data FileOptions + = FileOptions + { optFileEncoding :: Maybe TextEncoding + } deriving (Show) + +instance Default FileOptions where + def = + FileOptions + { optFileEncoding = Nothing + } diff --git a/src/Numeric/Optimization/MIP/FileUtils.hs b/src/Numeric/Optimization/MIP/FileUtils.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/FileUtils.hs @@ -0,0 +1,30 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE CPP #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.FileUtils +-- Copyright : (c) Masahiro Sakai 2018 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.FileUtils + ( ParseError + ) where + +#if MIN_VERSION_megaparsec(6,0,0) +import Data.Void +#endif +import qualified Text.Megaparsec as MP + +#if MIN_VERSION_megaparsec(7,0,0) +type ParseError s = MP.ParseErrorBundle s Void +#elif MIN_VERSION_megaparsec(6,0,0) +type ParseError s = MP.ParseError (MP.Token s) Void +#else +type ParseError s = MP.ParseError (MP.Token s) MP.Dec +#endif diff --git a/src/Numeric/Optimization/MIP/Internal/ProcessUtil.hs b/src/Numeric/Optimization/MIP/Internal/ProcessUtil.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Internal/ProcessUtil.hs @@ -0,0 +1,94 @@ +{-# LANGUAGE CPP #-} +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Internal.ProcessUtil +-- Copyright : (c) Masahiro Sakai 2014 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Internal.ProcessUtil + ( runProcessWithOutputCallback + ) where + +import Control.Concurrent +import Control.Concurrent.STM +import Control.Exception (SomeException, try, mask, throwIO) +import qualified Control.Exception as C +import Control.Monad +import Foreign.C +import System.Exit +import System.IO +import System.IO.Error +import System.Process + +#ifdef __GLASGOW_HASKELL__ +import GHC.IO.Exception ( IOErrorType(..), IOException(..) ) +#endif + +runProcessWithOutputCallback + :: FilePath -- ^ Filename of the executable (see 'proc' for details) + -> [String] -- ^ any arguments + -> String -- ^ standard input + -> (String -> IO ()) -- ^ callback function which is called when a line is read from stdout + -> (String -> IO ()) -- ^ callback function which is called when a line is read from stderr + -> IO ExitCode +runProcessWithOutputCallback cmd args input putMsg putErr = do + (Just inh, Just outh, Just errh, processh) <- createProcess + (proc cmd args) + { std_in = CreatePipe + , std_out = CreatePipe + , std_err = CreatePipe + } + req <- newEmptyTMVarIO + let f act = atomically (putTMVar req act) + m1 = forever (hGetLine outh >>= \s -> f (putMsg s)) + `catchIOError` (\e -> if isEOFError e then return () else ioError e) + m2 = forever (hGetLine errh >>= \s -> f (putErr s)) + `catchIOError` (\e -> if isEOFError e then return () else ioError e) + withForkWait m1 $ \waitOut -> do + withForkWait m2 $ \waitErr -> do + -- now write any input + unless (null input) $ + ignoreSigPipe $ hPutStr inh input + -- hClose performs implicit hFlush, and thus may trigger a SIGPIPE + ignoreSigPipe $ hClose inh + + hSetBuffering outh LineBuffering + hSetBuffering errh LineBuffering + let loop = join $ atomically $ msum $ + [ do act <- takeTMVar req + return $ act >> loop + , do guard =<< isEmptyTMVar req + waitOut + waitErr + return $ return () + ] + loop + hClose outh + hClose errh + waitForProcess processh + +withForkWait :: IO () -> (STM () -> IO a) -> IO a +withForkWait async body = do + waitVar <- newEmptyTMVarIO :: IO (TMVar (Either SomeException ())) + mask $ \restore -> do + tid <- forkIO $ try (restore async) >>= \v -> atomically (putTMVar waitVar v) + let wait = takeTMVar waitVar >>= either throwSTM return + restore (body wait) `C.onException` killThread tid + +ignoreSigPipe :: IO () -> IO () +#if defined(__GLASGOW_HASKELL__) +ignoreSigPipe = C.handle $ \e -> case e of + IOError { ioe_type = ResourceVanished + , ioe_errno = Just ioe } + | Errno ioe == ePIPE -> return () + _ -> throwIO e +#else +ignoreSigPipe = id +#endif diff --git a/src/Numeric/Optimization/MIP/Internal/Util.hs b/src/Numeric/Optimization/MIP/Internal/Util.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Internal/Util.hs @@ -0,0 +1,95 @@ +{-# OPTIONS_GHC -Wall #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Internal.Util +-- Copyright : (c) Masahiro Sakai 2011-2012 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : portable +-- +-- Some utility functions. +-- +----------------------------------------------------------------------------- + +module Numeric.Optimization.MIP.Internal.Util where + +import Control.Monad +import Data.Ratio +import Data.Set (Set) +import qualified Data.Set as Set +import System.IO +import GHC.IO.Encoding + +-- | Combining two @Maybe@ values using given function. +combineMaybe :: (a -> a -> a) -> Maybe a -> Maybe a -> Maybe a +combineMaybe _ Nothing y = y +combineMaybe _ x Nothing = x +combineMaybe f (Just x) (Just y) = Just (f x y) + +-- | is the number integral? +-- +-- @ +-- isInteger x = fromInteger (round x) == x +-- @ +isInteger :: RealFrac a => a -> Bool +isInteger x = fromInteger (round x) == x + +-- | fractional part +-- +-- @ +-- fracPart x = x - fromInteger (floor x) +-- @ +fracPart :: RealFrac a => a -> a +fracPart x = x - fromInteger (floor x) + +showRational :: Bool -> Rational -> String +showRational asRatio v + | denominator v == 1 = show (numerator v) + | asRatio = show (numerator v) ++ "/" ++ show (denominator v) + | otherwise = show (fromRational v :: Double) + +showRationalAsFiniteDecimal :: Rational -> Maybe String +showRationalAsFiniteDecimal x = do + let a :: Integer + (a,b) = properFraction (abs x) + s1 = if x < 0 then "-" else "" + s2 = show a + s3 <- if b == 0 + then return ".0" + else liftM ("." ++ ) $ loop Set.empty b + return $ s1 ++ s2 ++ s3 + where + loop :: Set Rational -> Rational -> Maybe String + loop _ 0 = return "" + loop rs r + | r `Set.member` rs = mzero + | otherwise = do + let a :: Integer + (a,b) = properFraction (r * 10) + s <- loop (Set.insert r rs) b + return $ show a ++ s + +{-# INLINE revSequence #-} +revSequence :: Monad m => [m a] -> m [a] +revSequence = go [] + where + go xs [] = return xs + go xs (m:ms) = do + x <- m + go (x:xs) ms + +{-# INLINE revMapM #-} +revMapM :: Monad m => (a -> m b) -> ([a] -> m [b]) +revMapM f = revSequence . map f + +{-# INLINE revForM #-} +revForM :: Monad m => [a] -> (a -> m b) -> m [b] +revForM = flip revMapM + +setEncodingChar8 :: IO () +setEncodingChar8 = do + setLocaleEncoding char8 + setForeignEncoding char8 + setFileSystemEncoding char8 diff --git a/src/Numeric/Optimization/MIP/LPFile.hs b/src/Numeric/Optimization/MIP/LPFile.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/LPFile.hs @@ -0,0 +1,743 @@ +{-# OPTIONS_GHC -Wall -fno-warn-unused-do-bind #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.LPFile +-- Copyright : (c) Masahiro Sakai 2011-2014 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +-- A CPLEX .lp format parser library. +-- +-- References: +-- +-- * +-- +-- * +-- +-- * +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.LPFile + ( parseString + , parseFile + , ParseError + , parser + , render + ) where + +import Control.Applicative hiding (many) +import Control.Exception (throwIO) +import Control.Monad +import Control.Monad.Writer +import Control.Monad.ST +import Data.Char +import Data.Default.Class +import Data.Interned +import Data.List +import Data.Maybe +#if !MIN_VERSION_base(4,9,0) +import Data.Monoid +#endif +import Data.Scientific (Scientific) +import Data.Map (Map) +import qualified Data.Map as Map +import Data.Set (Set) +import qualified Data.Set as Set +import Data.STRef +import Data.String +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +import Data.Text.Lazy.Builder (Builder) +import qualified Data.Text.Lazy.Builder as B +import qualified Data.Text.Lazy.Builder.Scientific as B +import qualified Data.Text.Lazy.IO as TLIO +import Data.OptDir +import System.IO +#if MIN_VERSION_megaparsec(6,0,0) +import Text.Megaparsec hiding (label, skipManyTill, ParseError) +import Text.Megaparsec.Char hiding (string', char') +import qualified Text.Megaparsec.Char.Lexer as P +#else +import Text.Megaparsec hiding (label, string', char', ParseError) +import qualified Text.Megaparsec.Lexer as P +import Text.Megaparsec.Prim (MonadParsec ()) +#endif + +import qualified Numeric.Optimization.MIP.Base as MIP +import Numeric.Optimization.MIP.FileUtils (ParseError) +import Numeric.Optimization.MIP.Internal.Util (combineMaybe) + +-- --------------------------------------------------------------------------- + +#if MIN_VERSION_megaparsec(6,0,0) +type C e s m = (MonadParsec e s m, Token s ~ Char, IsString (Tokens s)) +#else +type C e s m = (MonadParsec e s m, Token s ~ Char) +#endif + +-- | Parse a string containing LP file data. +-- The source name is only | used in error messages and may be the empty string. +#if MIN_VERSION_megaparsec(6,0,0) +parseString :: (Stream s, Token s ~ Char, IsString (Tokens s)) => MIP.FileOptions -> String -> s -> Either (ParseError s) (MIP.Problem Scientific) +#else +parseString :: (Stream s, Token s ~ Char) => MIP.FileOptions -> String -> s -> Either (ParseError s) (MIP.Problem Scientific) +#endif +parseString _ = parse (parser <* eof) + +-- | Parse a file containing LP file data. +parseFile :: MIP.FileOptions -> FilePath -> IO (MIP.Problem Scientific) +parseFile opt fname = do + h <- openFile fname ReadMode + case MIP.optFileEncoding opt of + Nothing -> return () + Just enc -> hSetEncoding h enc + ret <- parse (parser <* eof) fname <$> TLIO.hGetContents h + case ret of + Left e -> throwIO (e :: ParseError TL.Text) + Right a -> return a + +-- --------------------------------------------------------------------------- + +#if MIN_VERSION_megaparsec(7,0,0) +anyChar :: C e s m => m Char +anyChar = anySingle +#endif + +char' :: C e s m => Char -> m Char +char' c = (char c <|> char (toUpper c)) show c + +string' :: C e s m => String -> m () +string' s = mapM_ char' s show s + +sep :: C e s m => m () +sep = skipMany ((comment >> return ()) <|> (spaceChar >> return ())) + +comment :: C e s m => m () +comment = do + char '\\' + skipManyTill anyChar (try newline) + +tok :: C e s m => m a -> m a +tok p = do + x <- p + sep + return x + +ident :: C e s m => m String +ident = tok $ do + x <- letterChar <|> oneOf syms1 + xs <- many (alphaNumChar <|> oneOf syms2) + let s = x:xs + guard $ map toLower s `Set.notMember` reserved + return s + where + syms1 = "!\"#$%&()/,;?@_`'{}|~" + syms2 = '.' : syms1 + +variable :: C e s m => m MIP.Var +variable = liftM MIP.toVar ident + +label :: C e s m => m MIP.Label +label = do + name <- ident + tok $ char ':' + return $! T.pack name + +reserved :: Set String +reserved = Set.fromList + [ "bound", "bounds" + , "gen", "general", "generals" + , "bin", "binary", "binaries" + , "semi", "semi-continuous", "semis" + , "sos" + , "end" + , "subject" + ] + +-- --------------------------------------------------------------------------- + +-- | LP file parser +#if MIN_VERSION_megaparsec(6,0,0) +parser :: (MonadParsec e s m, Token s ~ Char, IsString (Tokens s)) => m (MIP.Problem Scientific) +#else +parser :: (MonadParsec e s m, Token s ~ Char) => m (MIP.Problem Scientific) +#endif +parser = do + name <- optional $ try $ do + space + string' "\\* Problem: " + liftM fromString $ manyTill anyChar (try (string " *\\\n")) + sep + obj <- problem + + cs <- liftM concat $ many $ msum $ + [ liftM (map Left) constraintSection + , liftM (map Left) lazyConstraintsSection + , liftM (map Right) userCutsSection + ] + + bnds <- option Map.empty (try boundsSection) + exvs <- many (liftM Left generalSection <|> liftM Right binarySection) + let ints = Set.fromList $ concat [x | Left x <- exvs] + bins = Set.fromList $ concat [x | Right x <- exvs] + bnds2 <- return $ Map.unionWith MIP.intersectBounds + bnds (Map.fromAscList [(v, (MIP.Finite 0, MIP.Finite 1)) | v <- Set.toAscList bins]) + scs <- liftM Set.fromList $ option [] (try semiSection) + + ss <- option [] (try sosSection) + end + let vs = Set.unions $ map MIP.vars cs ++ + [ Map.keysSet bnds2 + , ints + , bins + , scs + , MIP.vars obj + , MIP.vars ss + ] + isInt v = v `Set.member` ints || v `Set.member` bins + isSemi v = v `Set.member` scs + return $ + MIP.Problem + { MIP.name = name + , MIP.objectiveFunction = obj + , MIP.constraints = [c | Left c <- cs] + , MIP.userCuts = [c | Right c <- cs] + , MIP.sosConstraints = ss + , MIP.varType = Map.fromAscList + [ ( v + , if isInt v then + if isSemi v then MIP.SemiIntegerVariable + else MIP.IntegerVariable + else + if isSemi v then MIP.SemiContinuousVariable + else MIP.ContinuousVariable + ) + | v <- Set.toAscList vs ] + , MIP.varBounds = Map.fromAscList [ (v, Map.findWithDefault MIP.defaultBounds v bnds2) | v <- Set.toAscList vs] + } + +problem :: C e s m => m (MIP.ObjectiveFunction Scientific) +problem = do + flag <- (try minimize >> return OptMin) + <|> (try maximize >> return OptMax) + name <- optional (try label) + obj <- expr + return def{ MIP.objLabel = name, MIP.objDir = flag, MIP.objExpr = obj } + +minimize, maximize :: C e s m => m () +minimize = tok $ string' "min" >> optional (string' "imize") >> return () +maximize = tok $ string' "max" >> optional (string' "imize") >> return () + +end :: C e s m => m () +end = tok $ string' "end" + +-- --------------------------------------------------------------------------- + +constraintSection :: C e s m => m [MIP.Constraint Scientific] +constraintSection = subjectTo >> many (try (constraint False)) + +subjectTo :: C e s m => m () +subjectTo = msum + [ try $ tok (string' "subject") >> tok (string' "to") + , try $ tok (string' "such") >> tok (string' "that") + , try $ tok (string' "st") + , try $ tok (string' "s") >> optional (tok (char '.')) >> tok (string' "t") + >> tok (char '.') >> return () + ] + +constraint :: C e s m => Bool -> m (MIP.Constraint Scientific) +constraint isLazy = do + name <- optional (try label) + g <- optional $ try indicator + + -- It seems that CPLEX allows empty lhs, but GLPK rejects it. + e <- expr + op <- relOp + s <- option 1 sign + rhs <- liftM (s*) number + + let (lb,ub) = + case op of + MIP.Le -> (MIP.NegInf, MIP.Finite rhs) + MIP.Ge -> (MIP.Finite rhs, MIP.PosInf) + MIP.Eql -> (MIP.Finite rhs, MIP.Finite rhs) + + return $ MIP.Constraint + { MIP.constrLabel = name + , MIP.constrIndicator = g + , MIP.constrExpr = e + , MIP.constrLB = lb + , MIP.constrUB = ub + , MIP.constrIsLazy = isLazy + } + +relOp :: C e s m => m MIP.RelOp +relOp = tok $ msum + [ char '<' >> optional (char '=') >> return MIP.Le + , char '>' >> optional (char '=') >> return MIP.Ge + , char '=' >> msum [ char '<' >> return MIP.Le + , char '>' >> return MIP.Ge + , return MIP.Eql + ] + ] + +indicator :: C e s m => m (MIP.Var, Scientific) +indicator = do + var <- variable + tok (char '=') + val <- tok ((char '0' >> return 0) <|> (char '1' >> return 1)) + tok $ string "->" + return (var, val) + +lazyConstraintsSection :: C e s m => m [MIP.Constraint Scientific] +lazyConstraintsSection = do + tok $ string' "lazy" + tok $ string' "constraints" + many $ try $ constraint True + +userCutsSection :: C e s m => m [MIP.Constraint Scientific] +userCutsSection = do + tok $ string' "user" + tok $ string' "cuts" + many $ try $ constraint False + +type Bounds2 c = (Maybe (MIP.BoundExpr c), Maybe (MIP.BoundExpr c)) + +boundsSection :: C e s m => m (Map MIP.Var (MIP.Bounds Scientific)) +boundsSection = do + tok $ string' "bound" >> optional (char' 's') + liftM (Map.map g . Map.fromListWith f) $ many (try bound) + where + f (lb1,ub1) (lb2,ub2) = (combineMaybe max lb1 lb2, combineMaybe min ub1 ub2) + g (lb, ub) = ( fromMaybe MIP.defaultLB lb + , fromMaybe MIP.defaultUB ub + ) + +bound :: C e s m => m (MIP.Var, Bounds2 Scientific) +bound = msum + [ try $ do + v <- try variable + msum + [ do + op <- relOp + b <- boundExpr + return + ( v + , case op of + MIP.Le -> (Nothing, Just b) + MIP.Ge -> (Just b, Nothing) + MIP.Eql -> (Just b, Just b) + ) + , do + tok $ string' "free" + return (v, (Just MIP.NegInf, Just MIP.PosInf)) + ] + , do + b1 <- liftM Just boundExpr + op1 <- relOp + guard $ op1 == MIP.Le + v <- variable + b2 <- option Nothing $ do + op2 <- relOp + guard $ op2 == MIP.Le + liftM Just boundExpr + return (v, (b1, b2)) + ] + +boundExpr :: C e s m => m (MIP.BoundExpr Scientific) +boundExpr = msum + [ try (tok (char '+') >> inf >> return MIP.PosInf) + , try (tok (char '-') >> inf >> return MIP.NegInf) + , do + s <- option 1 sign + x <- number + return $ MIP.Finite (s*x) + ] + +inf :: C e s m => m () +inf = tok (string "inf" >> optional (string "inity")) >> return () + +-- --------------------------------------------------------------------------- + +generalSection :: C e s m => m [MIP.Var] +generalSection = do + tok $ string' "gen" >> optional (string' "eral" >> optional (string' "s")) + many (try variable) + +binarySection :: C e s m => m [MIP.Var] +binarySection = do + tok $ string' "bin" >> optional (string' "ar" >> (string' "y" <|> string' "ies")) + many (try variable) + +semiSection :: C e s m => m [MIP.Var] +semiSection = do + tok $ string' "semi" >> optional (string' "-continuous" <|> string' "s") + many (try variable) + +sosSection :: C e s m => m [MIP.SOSConstraint Scientific] +sosSection = do + tok $ string' "sos" + many $ try $ do + (l,t) <- try (do{ l <- label; t <- typ; return (Just l, t) }) + <|> (do{ t <- typ; return (Nothing, t) }) + xs <- many $ try $ do + v <- variable + tok $ char ':' + w <- number + return (v,w) + return $ MIP.SOSConstraint l t xs + where + typ = do + t <- tok $ (char' 's' >> ((char '1' >> return MIP.S1) <|> (char '2' >> return MIP.S2))) + tok (string "::") + return t + +-- --------------------------------------------------------------------------- + +expr :: forall e s m. C e s m => m (MIP.Expr Scientific) +expr = try expr1 <|> return 0 + where + expr1 :: m (MIP.Expr Scientific) + expr1 = do + t <- term True + ts <- many (term False) + return $ foldr (+) 0 (t : ts) + +sign :: (C e s m, Num a) => m a +sign = tok ((char '+' >> return 1) <|> (char '-' >> return (-1))) + +term :: C e s m => Bool -> m (MIP.Expr Scientific) +term flag = do + s <- if flag then optional sign else liftM Just sign + c <- optional number + e <- liftM MIP.varExpr variable <|> qexpr + return $ case combineMaybe (*) s c of + Nothing -> e + Just d -> MIP.constExpr d * e + +qexpr :: C e s m => m (MIP.Expr Scientific) +qexpr = do + tok (char '[') + t <- qterm True + ts <- many (qterm False) + let e = MIP.Expr (t:ts) + tok (char ']') + -- Gurobi allows ommiting "/2" + (do mapM_ (tok . char) ("/2" :: String) -- Explicit type signature is necessary because the type of mapM_ in GHC-7.10 is generalized for arbitrary Foldable + return $ MIP.constExpr (1/2) * e) + <|> return e + +qterm :: C e s m => Bool -> m (MIP.Term Scientific) +qterm flag = do + s <- if flag then optional sign else liftM Just sign + c <- optional number + es <- do + e <- qfactor + es <- many (tok (char '*') >> qfactor) + return $ e ++ concat es + return $ case combineMaybe (*) s c of + Nothing -> MIP.Term 1 es + Just d -> MIP.Term d es + +qfactor :: C e s m => m [MIP.Var] +qfactor = do + v <- variable + msum [ tok (char '^') >> tok (char '2') >> return [v,v] + , return [v] + ] + +number :: forall e s m. C e s m => m Scientific +#if MIN_VERSION_megaparsec(6,0,0) +number = tok $ P.signed sep P.scientific +#else +number = tok $ P.signed sep P.number +#endif + +skipManyTill :: Alternative m => m a -> m end -> m () +skipManyTill p end' = scan + where + scan = (end' *> pure ()) <|> (p *> scan) + +-- --------------------------------------------------------------------------- + +type M a = Writer Builder a + +execM :: M a -> TL.Text +execM m = B.toLazyText $ execWriter m + +writeString :: T.Text -> M () +writeString s = tell $ B.fromText s + +writeChar :: Char -> M () +writeChar c = tell $ B.singleton c + +-- --------------------------------------------------------------------------- + +-- | Render a problem into a string. +render :: MIP.FileOptions -> MIP.Problem Scientific -> Either String TL.Text +render _ mip = Right $ execM $ render' $ normalize mip + +writeVar :: MIP.Var -> M () +writeVar v = writeString $ unintern v + +render' :: MIP.Problem Scientific -> M () +render' mip = do + case MIP.name mip of + Just name -> writeString $ "\\* Problem: " <> name <> " *\\\n" + Nothing -> return () + + let obj = MIP.objectiveFunction mip + + writeString $ + case MIP.objDir obj of + OptMin -> "MINIMIZE" + OptMax -> "MAXIMIZE" + writeChar '\n' + + renderLabel (MIP.objLabel obj) + renderExpr True (MIP.objExpr obj) + writeChar '\n' + + writeString "SUBJECT TO\n" + forM_ (MIP.constraints mip) $ \c -> do + unless (MIP.constrIsLazy c) $ do + renderConstraint c + writeChar '\n' + + let lcs = [c | c <- MIP.constraints mip, MIP.constrIsLazy c] + unless (null lcs) $ do + writeString "LAZY CONSTRAINTS\n" + forM_ lcs $ \c -> do + renderConstraint c + writeChar '\n' + + let cuts = [c | c <- MIP.userCuts mip] + unless (null cuts) $ do + writeString "USER CUTS\n" + forM_ cuts $ \c -> do + renderConstraint c + writeChar '\n' + + let ivs = MIP.integerVariables mip `Set.union` MIP.semiIntegerVariables mip + (bins,gens) = Set.partition (\v -> MIP.getBounds mip v == (MIP.Finite 0, MIP.Finite 1)) ivs + scs = MIP.semiContinuousVariables mip `Set.union` MIP.semiIntegerVariables mip + + writeString "BOUNDS\n" + forM_ (Map.toAscList (MIP.varBounds mip)) $ \(v, (lb,ub)) -> do + unless (v `Set.member` bins) $ do + renderBoundExpr lb + writeString " <= " + writeVar v + writeString " <= " + renderBoundExpr ub + writeChar '\n' + + unless (Set.null gens) $ do + writeString "GENERALS\n" + renderVariableList $ Set.toList gens + + unless (Set.null bins) $ do + writeString "BINARIES\n" + renderVariableList $ Set.toList bins + + unless (Set.null scs) $ do + writeString "SEMI-CONTINUOUS\n" + renderVariableList $ Set.toList scs + + unless (null (MIP.sosConstraints mip)) $ do + writeString "SOS\n" + forM_ (MIP.sosConstraints mip) $ \(MIP.SOSConstraint l typ xs) -> do + renderLabel l + writeString $ fromString $ show typ + writeString " ::" + forM_ xs $ \(v, r) -> do + writeString " " + writeVar v + writeString " : " + tell $ B.scientificBuilder r + writeChar '\n' + + writeString "END\n" + +-- FIXME: Gurobi は quadratic term が最後に一つある形式でないとダメっぽい +renderExpr :: Bool -> MIP.Expr Scientific -> M () +renderExpr isObj e = fill 80 (ts1 ++ ts2) + where + (ts,qts) = partition isLin (MIP.terms e) + isLin (MIP.Term _ []) = True + isLin (MIP.Term _ [_]) = True + isLin _ = False + + ts1 = map f ts + ts2 + | null qts = [] + | otherwise = + -- マイナスで始めるとSCIP 2.1.1 は「cannot have '-' in front of quadratic part ('[')」というエラーを出す + -- SCIP-3.1.0 does not allow spaces between '/' and '2'. + ["+ ["] ++ map g qts ++ [if isObj then "] /2" else "]"] + + f :: MIP.Term Scientific -> T.Text + f (MIP.Term c []) = showConstTerm c + f (MIP.Term c [v]) = showCoeff c <> fromString (MIP.fromVar v) + f _ = error "should not happen" + + g :: MIP.Term Scientific -> T.Text + g (MIP.Term c vs) = + (if isObj then showCoeff (2*c) else showCoeff c) <> + mconcat (intersperse " * " (map (fromString . MIP.fromVar) vs)) + +showValue :: Scientific -> T.Text +showValue = fromString . show + +showCoeff :: Scientific -> T.Text +showCoeff c = + if c' == 1 + then s + else s <> showValue c' <> " " + where + c' = abs c + s = if c >= 0 then "+ " else "- " + +showConstTerm :: Scientific -> T.Text +showConstTerm c = s <> showValue (abs c) + where + s = if c >= 0 then "+ " else "- " + +renderLabel :: Maybe MIP.Label -> M () +renderLabel l = + case l of + Nothing -> return () + Just s -> writeString s >> writeString ": " + +renderOp :: MIP.RelOp -> M () +renderOp MIP.Le = writeString "<=" +renderOp MIP.Ge = writeString ">=" +renderOp MIP.Eql = writeString "=" + +renderConstraint :: MIP.Constraint Scientific -> M () +renderConstraint c@MIP.Constraint{ MIP.constrExpr = e, MIP.constrLB = lb, MIP.constrUB = ub } = do + renderLabel (MIP.constrLabel c) + case MIP.constrIndicator c of + Nothing -> return () + Just (v,vval) -> do + writeVar v + writeString " = " + tell $ B.scientificBuilder vval + writeString " -> " + + renderExpr False e + writeChar ' ' + let (op, val) = + case (lb, ub) of + (MIP.NegInf, MIP.Finite x) -> (MIP.Le, x) + (MIP.Finite x, MIP.PosInf) -> (MIP.Ge, x) + (MIP.Finite x1, MIP.Finite x2) | x1==x2 -> (MIP.Eql, x1) + _ -> error "Numeric.Optimization.MIP.LPFile.renderConstraint: should not happen" + renderOp op + writeChar ' ' + tell $ B.scientificBuilder val + +renderBoundExpr :: MIP.BoundExpr Scientific -> M () +renderBoundExpr (MIP.Finite r) = tell $ B.scientificBuilder r +renderBoundExpr MIP.NegInf = writeString "-inf" +renderBoundExpr MIP.PosInf = writeString "+inf" + +renderVariableList :: [MIP.Var] -> M () +renderVariableList vs = fill 80 (map unintern vs) >> writeChar '\n' + +fill :: Int -> [T.Text] -> M () +fill width str = go str 0 + where + go [] _ = return () + go (x:xs) 0 = writeString x >> go xs (T.length x) + go (x:xs) w = + if w + 1 + T.length x <= width + then writeChar ' ' >> writeString x >> go xs (w + 1 + T.length x) + else writeChar '\n' >> go (x:xs) 0 + +-- --------------------------------------------------------------------------- + +{- +compileExpr :: Expr -> Maybe (Map Var Scientific) +compileExpr e = do + xs <- forM e $ \(Term c vs) -> + case vs of + [v] -> return (v, c) + _ -> mzero + return (Map.fromList xs) +-} + +-- --------------------------------------------------------------------------- + +normalize :: (Eq r, Num r) => MIP.Problem r -> MIP.Problem r +normalize = removeEmptyExpr . removeRangeConstraints + +removeRangeConstraints :: (Eq r, Num r) => MIP.Problem r -> MIP.Problem r +removeRangeConstraints prob = runST $ do + vsRef <- newSTRef $ MIP.variables prob + cntRef <- newSTRef (0::Int) + newvsRef <- newSTRef [] + + let gensym = do + vs <- readSTRef vsRef + let loop !c = do + let v = MIP.toVar ("~r_" ++ show c) + if v `Set.member` vs then + loop (c+1) + else do + writeSTRef cntRef $! c+1 + modifySTRef vsRef (Set.insert v) + return v + loop =<< readSTRef cntRef + + cs2 <- forM (MIP.constraints prob) $ \c -> do + case (MIP.constrLB c, MIP.constrUB c) of + (MIP.NegInf, MIP.Finite _) -> return c + (MIP.Finite _, MIP.PosInf) -> return c + (MIP.Finite x1, MIP.Finite x2) | x1 == x2 -> return c + (lb, ub) -> do + v <- gensym + modifySTRef newvsRef ((v, (lb,ub)) :) + return $ + c + { MIP.constrExpr = MIP.constrExpr c - MIP.varExpr v + , MIP.constrLB = MIP.Finite 0 + , MIP.constrUB = MIP.Finite 0 + } + + newvs <- liftM reverse $ readSTRef newvsRef + return $ + prob + { MIP.constraints = cs2 + , MIP.varType = MIP.varType prob `Map.union` Map.fromList [(v, MIP.ContinuousVariable) | (v,_) <- newvs] + , MIP.varBounds = MIP.varBounds prob `Map.union` (Map.fromList newvs) + } + +removeEmptyExpr :: Num r => MIP.Problem r -> MIP.Problem r +removeEmptyExpr prob = + prob + { MIP.objectiveFunction = obj{ MIP.objExpr = convertExpr (MIP.objExpr obj) } + , MIP.constraints = map convertConstr $ MIP.constraints prob + , MIP.userCuts = map convertConstr $ MIP.userCuts prob + } + where + obj = MIP.objectiveFunction prob + + convertExpr (MIP.Expr []) = MIP.Expr [MIP.Term 0 [MIP.toVar "x0"]] + convertExpr e = e + + convertConstr constr = + constr + { MIP.constrExpr = convertExpr $ MIP.constrExpr constr + } diff --git a/src/Numeric/Optimization/MIP/MPSFile.hs b/src/Numeric/Optimization/MIP/MPSFile.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/MPSFile.hs @@ -0,0 +1,922 @@ +{-# OPTIONS_GHC -Wall -fno-warn-unused-do-bind #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.MPSFile +-- Copyright : (c) Masahiro Sakai 2012-2014 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +-- A .mps format parser library. +-- +-- References: +-- +-- * +-- +-- * +-- +-- * +-- +-- * +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.MPSFile + ( parseString + , parseFile + , ParseError + , parser + , render + ) where + +#if !MIN_VERSION_base(4,8,0) +import Control.Applicative ((<$>), (<*)) +#endif +import Control.Exception (throwIO) +import Control.Monad +import Control.Monad.Writer +import Data.Default.Class +import Data.Maybe +#if !MIN_VERSION_base(4,9,0) +import Data.Monoid +#endif +import Data.Set (Set) +import qualified Data.Set as Set +import Data.Map (Map) +import qualified Data.Map as Map +import Data.Scientific +import Data.Interned +import Data.Interned.Text +import Data.String +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +import Data.Text.Lazy.Builder (Builder) +import qualified Data.Text.Lazy.Builder as B +import qualified Data.Text.Lazy.IO as TLIO +import System.IO +#if MIN_VERSION_megaparsec(6,0,0) +import Text.Megaparsec hiding (ParseError) +import Text.Megaparsec.Char hiding (string', newline) +import qualified Text.Megaparsec.Char as P +import qualified Text.Megaparsec.Char.Lexer as Lexer +#else +import qualified Text.Megaparsec as P +import Text.Megaparsec hiding (string', newline, ParseError) +import qualified Text.Megaparsec.Lexer as Lexer +import Text.Megaparsec.Prim (MonadParsec ()) +#endif + +import Data.OptDir +import qualified Numeric.Optimization.MIP.Base as MIP +import Numeric.Optimization.MIP.FileUtils (ParseError) + +type Column = MIP.Var +type Row = InternedText + +data BoundType + = LO -- lower bound + | UP -- upper bound + | FX -- variable is fixed at the specified value + | FR -- free variable (no lower or upper bound) + | MI -- infinite lower bound + | PL -- infinite upper bound + | BV -- variable is binary (equal 0 or 1) + | LI -- lower bound for integer variable + | UI -- upper bound for integer variable + | SC -- upper bound for semi-continuous variable + | SI -- upper bound for semi-integer variable + deriving (Eq, Ord, Show, Read, Enum, Bounded) + +-- --------------------------------------------------------------------------- + +#if MIN_VERSION_megaparsec(6,0,0) +type C e s m = (MonadParsec e s m, Token s ~ Char, IsString (Tokens s)) +#else +type C e s m = (MonadParsec e s m, Token s ~ Char) +#endif + +-- | Parse a string containing MPS file data. +-- The source name is only | used in error messages and may be the empty string. +#if MIN_VERSION_megaparsec(6,0,0) +parseString :: (Stream s, Token s ~ Char, IsString (Tokens s)) => MIP.FileOptions -> String -> s -> Either (ParseError s) (MIP.Problem Scientific) +#else +parseString :: (Stream s, Token s ~ Char) => MIP.FileOptions -> String -> s -> Either (ParseError s) (MIP.Problem Scientific) +#endif +parseString _ = parse (parser <* eof) + +-- | Parse a file containing MPS file data. +parseFile :: MIP.FileOptions -> FilePath -> IO (MIP.Problem Scientific) +parseFile opt fname = do + h <- openFile fname ReadMode + case MIP.optFileEncoding opt of + Nothing -> return () + Just enc -> hSetEncoding h enc + ret <- parse (parser <* eof) fname <$> TLIO.hGetContents h + case ret of + Left e -> throwIO (e :: ParseError TL.Text) + Right a -> return a + +-- --------------------------------------------------------------------------- + + +#if MIN_VERSION_megaparsec(7,0,0) +anyChar :: C e s m => m Char +anyChar = anySingle +#endif + +space' :: C e s m => m Char +space' = oneOf [' ', '\t'] + +spaces' :: C e s m => m () +spaces' = skipMany space' + +spaces1' :: C e s m => m () +spaces1' = skipSome space' + +commentline :: C e s m => m () +commentline = do + _ <- char '*' + _ <- manyTill anyChar P.newline + return () + +newline' :: C e s m => m () +newline' = do + spaces' + _ <- P.newline + skipMany commentline + return () + +tok :: C e s m => m a -> m a +tok p = do + x <- p + msum [eof, lookAhead (char '\n' >> return ()), spaces1'] + return x + +row :: C e s m => m Row +row = liftM intern ident + +column :: C e s m => m Column +column = liftM intern $ ident + +ident :: C e s m => m T.Text +ident = liftM fromString $ tok $ some $ noneOf [' ', '\t', '\n'] + +stringLn :: C e s m => String -> m () +stringLn s = string (fromString s) >> newline' + +number :: forall e s m. C e s m => m Scientific +#if MIN_VERSION_megaparsec(6,0,0) +number = tok $ Lexer.signed (return ()) Lexer.scientific +#else +number = tok $ Lexer.signed (return ()) Lexer.number +#endif + +-- --------------------------------------------------------------------------- + +-- | MPS file parser +#if MIN_VERSION_megaparsec(6,0,0) +parser :: (MonadParsec e s m, Token s ~ Char, IsString (Tokens s)) => m (MIP.Problem Scientific) +#else +parser :: (MonadParsec e s m, Token s ~ Char) => m (MIP.Problem Scientific) +#endif +parser = do + many commentline + + name <- nameSection + + -- http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_objsen.html + -- CPLEX extends the MPS standard by allowing two additional sections: OBJSEN and OBJNAME. + -- If these options are used, they must appear in order and as the first and second sections after the NAME section. + objsense <- optional $ objSenseSection + objname <- optional $ objNameSection + + rows <- rowsSection + + -- http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_usercuts.html + -- The order of sections must be ROWS USERCUTS. + usercuts <- option [] userCutsSection + + -- http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_lazycons.html + -- The order of sections must be ROWS USERCUTS LAZYCONS. + lazycons <- option [] lazyConsSection + + (cols, intvs1) <- colsSection + rhss <- rhsSection + rngs <- option Map.empty rangesSection + bnds <- option [] boundsSection + + -- http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_quadobj.html + -- Following the BOUNDS section, a QMATRIX section may be specified. + qobj <- msum [quadObjSection, qMatrixSection, return []] + + -- http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_sos.html + -- Note that in an MPS file, the SOS section must follow the BOUNDS section. + sos <- option [] sosSection + + -- http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_qcmatrix.html + -- QCMATRIX sections appear after the optional SOS section. + qterms <- liftM Map.fromList $ many qcMatrixSection + + -- http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_indicators.html + -- The INDICATORS section follows any quadratic constraint section and any quadratic objective section. + inds <- option Map.empty indicatorsSection + + string "ENDATA" + P.space + + let objrow = + case objname of + Nothing -> head [r | (Nothing, r) <- rows] -- XXX + Just r -> intern r + objdir = + case objsense of + Nothing -> OptMin + Just d -> d + vs = Map.keysSet cols `Set.union` Set.fromList [col | (_,col,_) <- bnds] + intvs2 = Set.fromList [col | (t,col,_) <- bnds, t `elem` [BV,LI,UI]] + scvs = Set.fromList [col | (SC,col,_) <- bnds] + sivs = Set.fromList [col | (SI,col,_) <- bnds] + + let explicitBounds = Map.fromListWith f + [ case typ of + LO -> (col, (Just (MIP.Finite val), Nothing)) + UP -> (col, (Nothing, Just (MIP.Finite val))) + FX -> (col, (Just (MIP.Finite val), Just (MIP.Finite val))) + FR -> (col, (Just MIP.NegInf, Just MIP.PosInf)) + MI -> (col, (Just MIP.NegInf, Nothing)) + PL -> (col, (Nothing, Just MIP.PosInf)) + BV -> (col, (Just (MIP.Finite 0), Just (MIP.Finite 1))) + LI -> (col, (Just (MIP.Finite val), Nothing)) + UI -> (col, (Nothing, Just (MIP.Finite val))) + SC -> (col, (Nothing, Just (MIP.Finite val))) + SI -> (col, (Nothing, Just (MIP.Finite val))) + | (typ,col,val) <- bnds ] + where + f (a1,b1) (a2,b2) = (g a1 a2, g b1 b2) + g _ (Just x) = Just x + g x Nothing = x + + let bounds = Map.fromList + [ case Map.lookup v explicitBounds of + Nothing -> + if v `Set.member` intvs1 + then + -- http://eaton.math.rpi.edu/cplex90html/reffileformatscplex/reffileformatscplex9.html + -- If no bounds are specified for the variables within markers, bounds of 0 (zero) and 1 (one) are assumed. + (v, (MIP.Finite 0, MIP.Finite 1)) + else + (v, (MIP.Finite 0, MIP.PosInf)) + Just (Nothing, Just (MIP.Finite ub)) | ub < 0 -> + {- + http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_records.html + If no bounds are specified, CPLEX assumes a lower + bound of 0 (zero) and an upper bound of +∞. If only a + single bound is specified, the unspecified bound + remains at 0 or +∞, whichever applies, with one + exception. If an upper bound of less than 0 is + specified and no other bound is specified, the lower + bound is automatically set to -∞. CPLEX deviates + slightly from a convention used by some MPS readers + when it encounters an upper bound of 0 (zero). Rather + than automatically set this variable’s lower bound to + -∞, CPLEX accepts both a lower and upper bound of 0, + effectively fixing that variable at 0. CPLEX resets + the lower bound to -∞ only if the upper bound is less + than 0. A warning message is issued when this + exception is encountered. + -} + (v, (MIP.NegInf, MIP.Finite ub)) + {- + lp_solve uses 1 as default lower bound for semi-continuous variable. + + But Gurobi Optimizer uses 0 as default lower bound for semi-continuous variable. + Here we adopt Gurobi's way. + -} +{- + Just (Nothing, ub) | v `Set.member` scvs -> + (v, (MIP.Finite 1, fromMaybe MIP.PosInf ub)) +-} + Just (lb,ub) -> + (v, (fromMaybe (MIP.Finite 0) lb, fromMaybe MIP.PosInf ub)) + | v <- Set.toList vs ] + + let rowCoeffs :: Map Row (Map Column Scientific) + rowCoeffs = Map.fromListWith Map.union [(r, Map.singleton col coeff) | (col,m) <- Map.toList cols, (r,coeff) <- Map.toList m] + + let f :: Bool -> (Maybe MIP.RelOp, Row) -> [MIP.Constraint Scientific] + f _isLazy (Nothing, _row) = [] + f isLazy (Just op, r) = do + let lhs = [MIP.Term c [col] | (col,c) <- Map.toList (Map.findWithDefault Map.empty r rowCoeffs)] + ++ Map.findWithDefault [] r qterms + let rhs = Map.findWithDefault 0 r rhss + (lb,ub) = + case Map.lookup r rngs of + Nothing -> + case op of + MIP.Ge -> (MIP.Finite rhs, MIP.PosInf) + MIP.Le -> (MIP.NegInf, MIP.Finite rhs) + MIP.Eql -> (MIP.Finite rhs, MIP.Finite rhs) + Just rng -> + case op of + MIP.Ge -> (MIP.Finite rhs, MIP.Finite (rhs + abs rng)) + MIP.Le -> (MIP.Finite (rhs - abs rng), MIP.Finite rhs) + MIP.Eql -> + if rng < 0 + then (MIP.Finite (rhs + rng), MIP.Finite rhs) + else (MIP.Finite rhs, MIP.Finite (rhs + rng)) + return $ + MIP.Constraint + { MIP.constrLabel = Just $ unintern r + , MIP.constrIndicator = Map.lookup r inds + , MIP.constrIsLazy = isLazy + , MIP.constrExpr = MIP.Expr lhs + , MIP.constrLB = lb + , MIP.constrUB = ub + } + + let mip = + MIP.Problem + { MIP.name = name + , MIP.objectiveFunction = def + { MIP.objDir = objdir + , MIP.objLabel = Just (unintern objrow) + , MIP.objExpr = MIP.Expr $ [MIP.Term c [col] | (col,m) <- Map.toList cols, c <- maybeToList (Map.lookup objrow m)] ++ qobj + } + , MIP.constraints = concatMap (f False) rows ++ concatMap (f True) lazycons + , MIP.sosConstraints = sos + , MIP.userCuts = concatMap (f False) usercuts + , MIP.varType = Map.fromAscList + [ ( v + , if v `Set.member` sivs then + MIP.SemiIntegerVariable + else if v `Set.member` intvs1 && v `Set.member` scvs then + MIP.SemiIntegerVariable + else if v `Set.member` intvs1 || v `Set.member` intvs2 then + MIP.IntegerVariable + else if v `Set.member` scvs then + MIP.SemiContinuousVariable + else + MIP.ContinuousVariable + ) + | v <- Set.toAscList vs ] + , MIP.varBounds = Map.fromAscList [(v, Map.findWithDefault MIP.defaultBounds v bounds) | v <- Set.toAscList vs] + } + + return mip + +nameSection :: C e s m => m (Maybe T.Text) +nameSection = do + string "NAME" + n <- optional $ try $ do + spaces1' + ident + newline' + return n + +objSenseSection :: C e s m => m OptDir +objSenseSection = do + try $ stringLn "OBJSENSE" + spaces1' + d <- (try (stringLn "MAX") >> return OptMax) + <|> (stringLn "MIN" >> return OptMin) + return d + +objNameSection :: C e s m => m T.Text +objNameSection = do + try $ stringLn "OBJNAME" + spaces1' + name <- ident + newline' + return name + +rowsSection :: C e s m => m [(Maybe MIP.RelOp, Row)] +rowsSection = do + try $ stringLn "ROWS" + rowsBody + +userCutsSection :: C e s m => m [(Maybe MIP.RelOp, Row)] +userCutsSection = do + try $ stringLn "USERCUTS" + rowsBody + +lazyConsSection :: C e s m => m [(Maybe MIP.RelOp, Row)] +lazyConsSection = do + try $ stringLn "LAZYCONS" + rowsBody + +rowsBody :: C e s m => m [(Maybe MIP.RelOp, Row)] +rowsBody = many $ do + spaces1' + op <- msum + [ char 'N' >> return Nothing + , char 'G' >> return (Just MIP.Ge) + , char 'L' >> return (Just MIP.Le) + , char 'E' >> return (Just MIP.Eql) + ] + spaces1' + name <- row + newline' + return (op, name) + +colsSection :: forall e s m. C e s m => m (Map Column (Map Row Scientific), Set Column) +colsSection = do + try $ stringLn "COLUMNS" + body False Map.empty Set.empty + where + body :: Bool -> Map Column (Map Row Scientific) -> Set Column -> m (Map Column (Map Row Scientific), Set Column) + body isInt rs ivs = msum + [ do _ <- spaces1' + x <- ident + msum + [ do isInt' <- try intMarker + body isInt' rs ivs + , do (k,v) <- entry x + let rs' = Map.insertWith Map.union k v rs + ivs' = if isInt then Set.insert k ivs else ivs + seq rs' $ seq ivs' $ body isInt rs' ivs' + ] + , return (rs, ivs) + ] + + intMarker :: m Bool + intMarker = do + string "'MARKER'" + spaces1' + b <- (try (string "'INTORG'") >> return True) + <|> (string "'INTEND'" >> return False) + newline' + return b + + entry :: T.Text -> m (Column, Map Row Scientific) + entry x = do + let col = intern x + rv1 <- rowAndVal + opt <- optional rowAndVal + newline' + case opt of + Nothing -> return (col, rv1) + Just rv2 -> return (col, Map.union rv1 rv2) + +rowAndVal :: C e s m => m (Map Row Scientific) +rowAndVal = do + r <- row + val <- number + return $ Map.singleton r val + +rhsSection :: C e s m => m (Map Row Scientific) +rhsSection = do + try $ stringLn "RHS" + liftM Map.unions $ many entry + where + entry = do + spaces1' + _name <- ident + rv1 <- rowAndVal + opt <- optional rowAndVal + newline' + case opt of + Nothing -> return rv1 + Just rv2 -> return $ Map.union rv1 rv2 + +rangesSection :: C e s m => m (Map Row Scientific) +rangesSection = do + try $ stringLn "RANGES" + liftM Map.unions $ many entry + where + entry = do + spaces1' + _name <- ident + rv1 <- rowAndVal + opt <- optional rowAndVal + newline' + case opt of + Nothing -> return rv1 + Just rv2 -> return $ Map.union rv1 rv2 + +boundsSection :: C e s m => m [(BoundType, Column, Scientific)] +boundsSection = do + try $ stringLn "BOUNDS" + many entry + where + entry = do + spaces1' + typ <- boundType + _name <- ident + col <- column + val <- if typ `elem` [FR, BV, MI, PL] + then return 0 + else number + newline' + return (typ, col, val) + +boundType :: C e s m => m BoundType +boundType = tok $ do + msum [try (string (fromString (show k))) >> return k | k <- [minBound..maxBound]] + +sosSection :: forall e s m. C e s m => m [MIP.SOSConstraint Scientific] +sosSection = do + try $ stringLn "SOS" + many entry + where + entry = do + spaces1' + typ <- (try (string "S1") >> return MIP.S1) + <|> (string "S2" >> return MIP.S2) + spaces1' + name <- ident + newline' + xs <- many (try identAndVal) + return $ MIP.SOSConstraint{ MIP.sosLabel = Just name, MIP.sosType = typ, MIP.sosBody = xs } + + identAndVal :: m (Column, Scientific) + identAndVal = do + spaces1' + col <- column + val <- number + newline' + return (col, val) + +quadObjSection :: C e s m => m [MIP.Term Scientific] +quadObjSection = do + try $ stringLn "QUADOBJ" + many entry + where + entry = do + spaces1' + col1 <- column + col2 <- column + val <- number + newline' + return $ MIP.Term (if col1 /= col2 then val else val / 2) [col1, col2] + +qMatrixSection :: C e s m => m [MIP.Term Scientific] +qMatrixSection = do + try $ stringLn "QMATRIX" + many entry + where + entry = do + spaces1' + col1 <- column + col2 <- column + val <- number + newline' + return $ MIP.Term (val / 2) [col1, col2] + +qcMatrixSection :: C e s m => m (Row, [MIP.Term Scientific]) +qcMatrixSection = do + try $ string "QCMATRIX" + spaces1' + r <- row + newline' + xs <- many entry + return (r, xs) + where + entry = do + spaces1' + col1 <- column + col2 <- column + val <- number + newline' + return $ MIP.Term val [col1, col2] + +indicatorsSection :: C e s m => m (Map Row (Column, Scientific)) +indicatorsSection = do + try $ stringLn "INDICATORS" + liftM Map.fromList $ many entry + where + entry = do + spaces1' + string "IF" + spaces1' + r <- row + var <- column + val <- number + newline' + return (r, (var, val)) + +-- --------------------------------------------------------------------------- + +type M a = Writer Builder a + +execM :: M a -> TL.Text +execM m = B.toLazyText $ execWriter m + +writeText :: T.Text -> M () +writeText s = tell $ B.fromText s + +writeChar :: Char -> M () +writeChar c = tell $ B.singleton c + +-- --------------------------------------------------------------------------- + +render :: MIP.FileOptions -> MIP.Problem Scientific -> Either String TL.Text +render _ mip | not (checkAtMostQuadratic mip) = Left "Expression must be atmost quadratic" +render _ mip = Right $ execM $ render' $ nameRows mip + +render' :: MIP.Problem Scientific -> M () +render' mip = do + let probName = fromMaybe "" (MIP.name mip) + + -- NAME section + -- The name starts in column 15 in fixed formats. + writeSectionHeader $ "NAME" <> T.replicate 10 " " <> probName + + let MIP.ObjectiveFunction + { MIP.objLabel = Just objName + , MIP.objDir = dir + , MIP.objExpr = obj + } = MIP.objectiveFunction mip + + -- OBJSENSE section + -- Note: GLPK-4.48 does not support this section. + writeSectionHeader "OBJSENSE" + case dir of + OptMin -> writeFields ["MIN"] + OptMax -> writeFields ["MAX"] + +{- + -- OBJNAME section + -- Note: GLPK-4.48 does not support this section. + writeSectionHeader "OBJNAME" + writeFields [objName] +-} + + let splitRange c = + case (MIP.constrLB c, MIP.constrUB c) of + (MIP.Finite x, MIP.PosInf) -> ((MIP.Ge, x), Nothing) + (MIP.NegInf, MIP.Finite x) -> ((MIP.Le, x), Nothing) + (MIP.Finite x1, MIP.Finite x2) + | x1 == x2 -> ((MIP.Eql, x1), Nothing) + | x1 < x2 -> ((MIP.Eql, x1), Just (x2 - x1)) + _ -> error "invalid constraint bound" + + let renderRows cs = do + forM_ cs $ \c -> do + let ((op,_), _) = splitRange c + let s = case op of + MIP.Le -> "L" + MIP.Ge -> "G" + MIP.Eql -> "E" + writeFields [s, fromJust $ MIP.constrLabel c] + + -- ROWS section + writeSectionHeader "ROWS" + writeFields ["N", objName] + renderRows [c | c <- MIP.constraints mip, not (MIP.constrIsLazy c)] + + -- USERCUTS section + unless (null (MIP.userCuts mip)) $ do + writeSectionHeader "USERCUTS" + renderRows (MIP.userCuts mip) + + -- LAZYCONS section + let lcs = [c | c <- MIP.constraints mip, MIP.constrIsLazy c] + unless (null lcs) $ do + writeSectionHeader "LAZYCONS" + renderRows lcs + + -- COLUMNS section + writeSectionHeader "COLUMNS" + let cols :: Map Column (Map T.Text Scientific) + cols = Map.fromListWith Map.union + [ (v, Map.singleton l d) + | (Just l, xs) <- + (Just objName, obj) : + [(MIP.constrLabel c, lhs) | c <- MIP.constraints mip ++ MIP.userCuts mip, let lhs = MIP.constrExpr c] + , MIP.Term d [v] <- MIP.terms xs + ] + f col xs = + forM_ (Map.toList xs) $ \(r, d) -> do + writeFields ["", unintern col, r, showValue d] + ivs = MIP.integerVariables mip `Set.union` MIP.semiIntegerVariables mip + forM_ (Map.toList (Map.filterWithKey (\col _ -> col `Set.notMember` ivs) cols)) $ \(col, xs) -> f col xs + unless (Set.null ivs) $ do + writeFields ["", "MARK0000", "'MARKER'", "", "'INTORG'"] + forM_ (Map.toList (Map.filterWithKey (\col _ -> col `Set.member` ivs) cols)) $ \(col, xs) -> f col xs + writeFields ["", "MARK0001", "'MARKER'", "", "'INTEND'"] + + -- RHS section + let rs = [(fromJust $ MIP.constrLabel c, rhs) | c <- MIP.constraints mip ++ MIP.userCuts mip, let ((_,rhs),_) = splitRange c, rhs /= 0] + writeSectionHeader "RHS" + forM_ rs $ \(name, val) -> do + writeFields ["", "rhs", name, showValue val] + + -- RANGES section + let rngs = [(fromJust $ MIP.constrLabel c, fromJust rng) | c <- MIP.constraints mip ++ MIP.userCuts mip, let ((_,_), rng) = splitRange c, isJust rng] + unless (null rngs) $ do + writeSectionHeader "RANGES" + forM_ rngs $ \(name, val) -> do + writeFields ["", "rhs", name, showValue val] + + -- BOUNDS section + writeSectionHeader "BOUNDS" + forM_ (Map.toList (MIP.varType mip)) $ \(col, vt) -> do + let (lb,ub) = MIP.getBounds mip col + case (lb,ub) of + (MIP.NegInf, MIP.PosInf) -> do + -- free variable (no lower or upper bound) + writeFields ["FR", "bound", unintern col] + + (MIP.Finite 0, MIP.Finite 1) | vt == MIP.IntegerVariable -> do + -- variable is binary (equal 0 or 1) + writeFields ["BV", "bound", unintern col] + + (MIP.Finite a, MIP.Finite b) | a == b -> do + -- variable is fixed at the specified value + writeFields ["FX", "bound", unintern col, showValue a] + + _ -> do + case lb of + MIP.PosInf -> error "should not happen" + MIP.NegInf -> do + -- Minus infinity + writeFields ["MI", "bound", unintern col] + MIP.Finite 0 | vt == MIP.ContinuousVariable -> return () + MIP.Finite a -> do + let t = case vt of + MIP.IntegerVariable -> "LI" -- lower bound for integer variable + _ -> "LO" -- Lower bound + writeFields [t, "bound", unintern col, showValue a] + + case ub of + MIP.NegInf -> error "should not happen" + MIP.PosInf | vt == MIP.ContinuousVariable -> return () + MIP.PosInf -> do + when (vt == MIP.SemiContinuousVariable || vt == MIP.SemiIntegerVariable) $ + error "cannot express +inf upper bound of semi-continuous or semi-integer variable" + writeFields ["PL", "bound", unintern col] -- Plus infinity + MIP.Finite a -> do + let t = case vt of + MIP.SemiContinuousVariable -> "SC" -- Upper bound for semi-continuous variable + MIP.SemiIntegerVariable -> + -- Gurobi uses "SC" while lpsolve uses "SI" for upper bound of semi-integer variable + "SC" + MIP.IntegerVariable -> "UI" -- Upper bound for integer variable + _ -> "UP" -- Upper bound + writeFields [t, "bound", unintern col, showValue a] + + -- QMATRIX section + -- Gurobiは対称行列になっていないと "qmatrix isn't symmetric" というエラーを発生させる + do let qm = Map.map (2*) $ quadMatrix obj + unless (Map.null qm) $ do + writeSectionHeader "QMATRIX" + forM_ (Map.toList qm) $ \(((v1,v2), val)) -> do + writeFields ["", unintern v1, unintern v2, showValue val] + + -- SOS section + unless (null (MIP.sosConstraints mip)) $ do + writeSectionHeader "SOS" + forM_ (MIP.sosConstraints mip) $ \sos -> do + let t = case MIP.sosType sos of + MIP.S1 -> "S1" + MIP.S2 -> "S2" + writeFields $ t : maybeToList (MIP.sosLabel sos) + forM_ (MIP.sosBody sos) $ \(var,val) -> do + writeFields ["", unintern var, showValue val] + + -- QCMATRIX section + let xs = [ (fromJust $ MIP.constrLabel c, qm) + | c <- MIP.constraints mip ++ MIP.userCuts mip + , let lhs = MIP.constrExpr c + , let qm = quadMatrix lhs + , not (Map.null qm) ] + unless (null xs) $ do + forM_ xs $ \(r, qm) -> do + -- The name starts in column 12 in fixed formats. + writeSectionHeader $ "QCMATRIX" <> T.replicate 3 " " <> r + forM_ (Map.toList qm) $ \((v1,v2), val) -> do + writeFields ["", unintern v1, unintern v2, showValue val] + + -- INDICATORS section + -- Note: Gurobi-5.6.3 does not support this section. + let ics = [c | c <- MIP.constraints mip, isJust $ MIP.constrIndicator c] + unless (null ics) $ do + writeSectionHeader "INDICATORS" + forM_ ics $ \c -> do + let Just (var,val) = MIP.constrIndicator c + writeFields ["IF", fromJust (MIP.constrLabel c), unintern var, showValue val] + + -- ENDATA section + writeSectionHeader "ENDATA" + +writeSectionHeader :: T.Text -> M () +writeSectionHeader s = writeText s >> writeChar '\n' + +-- Fields start in column 2, 5, 15, 25, 40 and 50 +writeFields :: [T.Text] -> M () +writeFields xs0 = f1 xs0 >> writeChar '\n' + where + -- columns 1-4 + f1 [] = return () + f1 [x] = writeChar ' ' >> writeText x + f1 (x:xs) = do + writeChar ' ' + writeText x + let len = T.length x + when (len < 2) $ writeText $ T.replicate (2 - len) " " + writeChar ' ' + f2 xs + + -- columns 5-14 + f2 [] = return () + f2 [x] = writeText x + f2 (x:xs) = do + writeText x + let len = T.length x + when (len < 9) $ writeText $ T.replicate (9 - len) " " + writeChar ' ' + f3 xs + + -- columns 15-24 + f3 [] = return () + f3 [x] = writeText x + f3 (x:xs) = do + writeText x + let len = T.length x + when (len < 9) $ writeText $ T.replicate (9 - len) " " + writeChar ' ' + f4 xs + + -- columns 25-39 + f4 [] = return () + f4 [x] = writeText x + f4 (x:xs) = do + writeText x + let len = T.length x + when (len < 14) $ writeText $ T.replicate (14 - len) " " + writeChar ' ' + f5 xs + + -- columns 40-49 + f5 [] = return () + f5 [x] = writeText x + f5 (x:xs) = do + writeText x + let len = T.length x + when (len < 19) $ writeText $ T.replicate (19 - len) " " + writeChar ' ' + f6 xs + + -- columns 50- + f6 [] = return () + f6 [x] = writeText x + f6 _ = error "MPSFile: >6 fields (this should not happen)" + +showValue :: Scientific -> T.Text +showValue = fromString . show + +nameRows :: MIP.Problem r -> MIP.Problem r +nameRows mip + = mip + { MIP.objectiveFunction = (MIP.objectiveFunction mip){ MIP.objLabel = Just objName' } + , MIP.constraints = f (MIP.constraints mip) [T.pack $ "row" ++ show n | n <- [(1::Int)..]] + , MIP.userCuts = f (MIP.userCuts mip) [T.pack $ "usercut" ++ show n | n <- [(1::Int)..]] + , MIP.sosConstraints = g (MIP.sosConstraints mip) [T.pack $ "sos" ++ show n | n <- [(1::Int)..]] + } + where + objName = MIP.objLabel $ MIP.objectiveFunction mip + used = Set.fromList $ catMaybes $ objName : [MIP.constrLabel c | c <- MIP.constraints mip ++ MIP.userCuts mip] ++ [MIP.sosLabel c | c <- MIP.sosConstraints mip] + objName' = fromMaybe (head [name | n <- [(1::Int)..], let name = T.pack ("obj" ++ show n), name `Set.notMember` used]) objName + + f [] _ = [] + f (c:cs) (name:names) + | isJust (MIP.constrLabel c) = c : f cs (name:names) + | name `Set.notMember` used = c{ MIP.constrLabel = Just name } : f cs names + | otherwise = f (c:cs) names + f _ [] = error "should not happen" + + g [] _ = [] + g (c:cs) (name:names) + | isJust (MIP.sosLabel c) = c : g cs (name:names) + | name `Set.notMember` used = c{ MIP.sosLabel = Just name } : g cs names + | otherwise = g (c:cs) names + g _ [] = error "should not happen" + +quadMatrix :: Fractional r => MIP.Expr r -> Map (MIP.Var, MIP.Var) r +quadMatrix e = Map.fromList $ do + let m = Map.fromListWith (+) [(if v1<=v2 then (v1,v2) else (v2,v1), c) | MIP.Term c [v1,v2] <- MIP.terms e] + ((v1,v2),c) <- Map.toList m + if v1==v2 then + [((v1,v2), c)] + else + [((v1,v2), c/2), ((v2,v1), c/2)] + +checkAtMostQuadratic :: forall r. MIP.Problem r -> Bool +checkAtMostQuadratic mip = all (all f . MIP.terms) es + where + es = MIP.objExpr (MIP.objectiveFunction mip) : + [lhs | c <- MIP.constraints mip ++ MIP.userCuts mip, let lhs = MIP.constrExpr c] + f :: MIP.Term r -> Bool + f (MIP.Term _ [_]) = True + f (MIP.Term _ [_,_]) = True + f _ = False + +-- --------------------------------------------------------------------------- diff --git a/src/Numeric/Optimization/MIP/Solution/CBC.hs b/src/Numeric/Optimization/MIP/Solution/CBC.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solution/CBC.hs @@ -0,0 +1,87 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solution.CBC +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solution.CBC + ( Solution (..) + , parse + , readFile + ) where + +import Prelude hiding (readFile, writeFile) +#if !MIN_VERSION_base(4,8,0) +import Control.Applicative +#endif + +import Control.Monad.Except +import Data.Interned (intern) +import Data.Map (Map) +import qualified Data.Map as Map +import Data.Scientific (Scientific) +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +import qualified Data.Text.Lazy.IO as TLIO +import Numeric.Optimization.MIP (Solution) +import qualified Numeric.Optimization.MIP as MIP + +parse :: TL.Text -> MIP.Solution Scientific +parse t = + case parse' $ TL.lines t of + Left e -> error e + Right x -> x + +parse' :: [TL.Text] -> Either String (MIP.Solution Scientific) +parse' (l1:ls) = do + (status, obj) <- + case TL.break ('-'==) l1 of + (s1,s2) -> + case TL.stripPrefix "- objective value " s2 of + Nothing -> throwError "fail to parse header" + Just s3 -> do + let s1' = TL.toStrict (TL.strip s1) + return + ( case Map.lookup s1' statusTable of + Just st -> st + Nothing -> + if T.isPrefixOf "Stopped on " s1' + then MIP.StatusUnknown + else MIP.StatusUnknown + , read (TL.unpack s3) + ) + let f :: [(MIP.Var, Scientific)] -> TL.Text -> Either String [(MIP.Var, Scientific)] + f vs t = + case TL.words t of + ("**":_no:var:val:_) -> return $ (intern (TL.toStrict var), read (TL.unpack val)) : vs + (_no:var:val:_) -> return $ (intern (TL.toStrict var), read (TL.unpack val)) : vs + [] -> return $ vs + _ -> throwError ("Numeric.Optimization.MIP.Solution.CBC: invalid line " ++ show t) + vs <- foldM f [] ls + return $ + MIP.Solution + { MIP.solStatus = status + , MIP.solObjectiveValue = Just obj + , MIP.solVariables = Map.fromList vs + } +parse' _ = throwError "must have >=1 lines" + +statusTable :: Map T.Text MIP.Status +statusTable = Map.fromList + [ ("Optimal", MIP.StatusOptimal) + , ("Unbounded", MIP.StatusInfeasibleOrUnbounded) + , ("Integer infeasible", MIP.StatusInfeasible) + , ("Infeasible", MIP.StatusInfeasible) + ] + +readFile :: FilePath -> IO (MIP.Solution Scientific) +readFile fname = parse <$> TLIO.readFile fname diff --git a/src/Numeric/Optimization/MIP/Solution/CPLEX.hs b/src/Numeric/Optimization/MIP/Solution/CPLEX.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solution/CPLEX.hs @@ -0,0 +1,163 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE OverloadedStrings #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solution.CPLEX +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solution.CPLEX + ( Solution (..) + , parse + , readFile + ) where + +import Prelude hiding (readFile) +#if !MIN_VERSION_base(4,8,0) +import Control.Applicative +#endif +import Data.Default.Class +import Data.Interned +import Data.IntMap (IntMap) +import qualified Data.IntMap as IntMap +import qualified Data.Map as Map +import Data.Maybe +import Data.Scientific (Scientific) +import Data.String +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +import qualified Text.XML as XML +import Text.XML.Cursor +import Numeric.Optimization.MIP (Solution) +import qualified Numeric.Optimization.MIP.Base as MIP + +parseDoc :: XML.Document -> MIP.Solution Scientific +parseDoc doc = + MIP.Solution + { MIP.solStatus = status + , MIP.solObjectiveValue = obj + , MIP.solVariables = Map.fromList vs + } + where + obj :: Maybe Scientific + obj = listToMaybe + $ fromDocument doc + $| element "CPLEXSolution" + &/ element "header" + >=> attribute "objectiveValue" + &| (read . T.unpack) + + status :: MIP.Status + status = head + $ fromDocument doc + $| element "CPLEXSolution" + &/ element "header" + >=> attribute "solutionStatusValue" + &| (flip (IntMap.findWithDefault MIP.StatusUnknown) table . read . T.unpack) + + f :: Cursor -> [(MIP.Var, Scientific)] + f x + | XML.NodeElement e <- node x = maybeToList $ do + let m = XML.elementAttributes e + name <- Map.lookup "name" m + value <- read . T.unpack <$> Map.lookup "value" m + return (intern name, value) + | otherwise = [] + vs = fromDocument doc + $| element "CPLEXSolution" + &/ element "variables" + &/ element "variable" + >=> f + +-- https://www.ibm.com/support/knowledgecenter/en/SSSA5P_12.7.0/ilog.odms.cplex.help/refcallablelibrary/macros/Solution_status_codes.html +table :: IntMap MIP.Status +table = IntMap.fromList + [ (1, MIP.StatusOptimal) -- CPX_STAT_OPTIMAL + , (2, MIP.StatusUnbounded) -- CPX_STAT_UNBOUNDED + , (3, MIP.StatusInfeasible) -- CPX_STAT_INFEASIBLE + , (4, MIP.StatusInfeasibleOrUnbounded) -- CPX_STAT_INForUNBD + , (5, MIP.StatusOptimal) -- CPX_STAT_OPTIMAL_INFEAS +{- + , (6, ) -- CPX_STAT_NUM_BEST + , (10, ) -- CPX_STAT_ABORT_IT_LIM + , (11, ) -- CPX_STAT_ABORT_TIME_LIM + , (12, ) -- CPX_STAT_ABORT_OBJ_LIM + , (13, ) -- CPX_STAT_ABORT_USER + , (14, ) -- CPX_STAT_FEASIBLE_RELAXED_SUM + , (15, ) -- CPX_STAT_OPTIMAL_RELAXED_SUM + , (16, ) -- CPX_STAT_FEASIBLE_RELAXED_INF + , (17, ) -- CPX_STAT_OPTIMAL_RELAXED_INF + , (18, ) -- CPX_STAT_FEASIBLE_RELAXED_QUAD + , (19, ) -- CPX_STAT_OPTIMAL_RELAXED_QUAD + , (20, ) -- CPX_STAT_OPTIMAL_FACE_UNBOUNDED + , (21, ) -- CPX_STAT_ABORT_PRIM_OBJ_LIM + , (22, ) -- CPX_STAT_ABORT_DUAL_OBJ_LIM + , (23, ) -- CPX_STAT_FEASIBLE +-} + , (24, MIP.StatusFeasible) -- CPX_STAT_FIRSTORDER +{- + , (25, ) -- CPX_STAT_ABORT_DETTIME_LIM + , (30, ) -- CPX_STAT_CONFLICT_FEASIBLE + , (31, ) -- CPX_STAT_CONFLICT_MINIMAL + , (32, ) -- CPX_STAT_CONFLICT_ABORT_CONTRADICTION + , (33, ) -- CPX_STAT_CONFLICT_ABORT_TIME_LIM + , (34, ) -- CPX_STAT_CONFLICT_ABORT_IT_LIM + , (35, ) -- CPX_STAT_CONFLICT_ABORT_NODE_LIM + , (36, ) -- CPX_STAT_CONFLICT_ABORT_OBJ_LIM + , (37, ) -- CPX_STAT_CONFLICT_ABORT_MEM_LIM + , (38, ) -- CPX_STAT_CONFLICT_ABORT_USER + , (39, ) -- CPX_STAT_CONFLICT_ABORT_DETTIME_LIM +-} + , (40, MIP.StatusInfeasibleOrUnbounded) -- CPX_STAT_BENDERS_MASTER_UNBOUNDED +-- , (41, ) -- CPX_STAT_BENDERS_NUM_BEST + , (101, MIP.StatusOptimal) -- CPXMIP_OPTIMAL + , (102, MIP.StatusOptimal) -- CPXMIP_OPTIMAL_TOL + , (103, MIP.StatusInfeasible) -- CPXMIP_INFEASIBLE +-- , (104, ) -- CPXMIP_SOL_LIM + , (105, MIP.StatusFeasible) -- CPXMIP_NODE_LIM_FEAS +-- , (106, ) -- CPXMIP_NODE_LIM_INFEAS + , (107, MIP.StatusFeasible) -- CPXMIP_TIME_LIM_FEAS +-- , (108, ) -- CPXMIP_TIME_LIM_INFEAS + , (109, MIP.StatusFeasible) -- CPXMIP_FAIL_FEAS +-- , (110, ) -- CPXMIP_FAIL_INFEAS + , (111, MIP.StatusFeasible) -- CPXMIP_MEM_LIM_FEAS +-- , (112, ) -- CPXMIP_MEM_LIM_INFEAS + , (113, MIP.StatusFeasible) -- CPXMIP_ABORT_FEAS +-- , (114, ) -- CPXMIP_ABORT_INFEAS + , (115, MIP.StatusOptimal) -- CPXMIP_OPTIMAL_INFEAS + , (116, MIP.StatusFeasible) -- CPXMIP_FAIL_FEAS_NO_TREE +-- , (117, ) -- CPXMIP_FAIL_INFEAS_NO_TREE + , (118, MIP.StatusUnbounded) -- CPXMIP_UNBOUNDED + , (119, MIP.StatusInfeasibleOrUnbounded) -- CPXMIP_INForUNBD +{- + , (120, ) -- CPXMIP_FEASIBLE_RELAXED_SUM + , (121, ) -- CPXMIP_OPTIMAL_RELAXED_SUM + , (122, ) -- CPXMIP_FEASIBLE_RELAXED_INF + , (123, ) -- CPXMIP_OPTIMAL_RELAXED_INF + , (124, ) -- CPXMIP_FEASIBLE_RELAXED_QUAD + , (125, ) -- CPXMIP_OPTIMAL_RELAXED_QUAD + , (126, ) -- CPXMIP_ABORT_RELAXED +-} + , (127, MIP.StatusFeasible) -- CPXMIP_FEASIBLE +-- , (128, ) -- CPXMIP_POPULATESOL_LIM + , (129, MIP.StatusOptimal) -- CPXMIP_OPTIMAL_POPULATED + , (130, MIP.StatusOptimal) -- CPXMIP_OPTIMAL_POPULATED_TOL + , (131, MIP.StatusFeasible) -- CPXMIP_DETTIME_LIM_FEAS +-- , (132, ) -- CPXMIP_DETTIME_LIM_INFEAS + , (133, MIP.StatusInfeasibleOrUnbounded) -- CPXMIP_ABORT_RELAXATION_UNBOUNDED + ] + +parse :: TL.Text -> MIP.Solution Scientific +parse t = parseDoc $ XML.parseText_ def t + +readFile :: FilePath -> IO (MIP.Solution Scientific) +readFile fname = parseDoc <$> XML.readFile def (fromString fname) + diff --git a/src/Numeric/Optimization/MIP/Solution/GLPK.hs b/src/Numeric/Optimization/MIP/Solution/GLPK.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solution/GLPK.hs @@ -0,0 +1,105 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solution.GLPK +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solution.GLPK + ( Solution (..) + , parse + , readFile + ) where + +import Prelude hiding (readFile, writeFile) +#if !MIN_VERSION_base(4,8,0) +import Control.Applicative +#endif +import Data.Interned (intern) +import Data.Map (Map) +import qualified Data.Map as Map +import Data.Scientific (Scientific) +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +import qualified Data.Text.Lazy.IO as TLIO +import Numeric.Optimization.MIP (Solution) +import qualified Numeric.Optimization.MIP as MIP + +parse :: TL.Text -> MIP.Solution Scientific +parse t = parse' $ TL.lines t + +parse' :: [TL.Text] -> MIP.Solution Scientific +parse' ls = + case parseHeaders ls of + (headers, ls2) -> + case parseColumns (skipRows ls2) of + (vs, _) -> + let status = Map.findWithDefault MIP.StatusUnknown (headers Map.! "Status") statusTable + objstr = headers Map.! "Objective" + objstr2 = + case T.findIndex ('='==) objstr of + Nothing -> objstr + Just idx -> T.drop (idx+1) objstr + obj = case reads (T.unpack objstr2) of + [] -> error "parse error" + (r,_):_ -> r + in MIP.Solution + { MIP.solStatus = status + , MIP.solObjectiveValue = Just obj + , MIP.solVariables = vs + } + +parseHeaders :: [TL.Text] -> (Map T.Text T.Text, [TL.Text]) +parseHeaders = f Map.empty + where + f _ [] = error "parse error" + f ret ("":ls) = (ret, ls) + f ret (l:ls) = + case TL.break (':'==) l of + (_, "") -> error "parse error" + (name, val) -> f (Map.insert (TL.toStrict name) (TL.toStrict (TL.strip (TL.tail val))) ret) ls + +skipRows :: [TL.Text] -> [TL.Text] +skipRows [] = error "parse error" +skipRows ("":ls) = ls +skipRows (_:ls) = skipRows ls + +parseColumns :: [TL.Text] -> (Map MIP.Var Scientific, [TL.Text]) +parseColumns (l1:l2:ls) + | l1 == " No. Column name Activity Lower bound Upper bound" + , l2 == "------ ------------ ------------- ------------- -------------" + = f [] ls + where + f :: [(MIP.Var, Scientific)] -> [TL.Text] -> (Map MIP.Var Scientific, [TL.Text]) + f _ [] = error "parse error" + f ret ("":ls2) = (Map.fromList ret, ls2) + f ret (l:ls2) = + case ws of + (_no : col : "*" : activity : _) -> f ((intern (TL.toStrict col), read (TL.unpack activity)) : ret) ls3 + (_no : col : activity : _) -> f ((intern (TL.toStrict col), read (TL.unpack activity)) : ret) ls3 + _ -> error "parse error" + where + (ws,ls3) = + case TL.words l of + ws1@(_:_:_:_) -> (ws1, ls2) + ws1@[_,_] -> (ws1 ++ TL.words (head ls2), tail ls2) + _ -> error "parse error" +parseColumns _ = error "parse error" + +statusTable :: Map T.Text MIP.Status +statusTable = Map.fromList + [ ("INTEGER OPTIMAL", MIP.StatusOptimal) + , ("INTEGER NON-OPTIMAL", MIP.StatusUnknown) + , ("INTEGER EMPTY", MIP.StatusInfeasible) + ] + +readFile :: FilePath -> IO (MIP.Solution Scientific) +readFile fname = parse <$> TLIO.readFile fname diff --git a/src/Numeric/Optimization/MIP/Solution/Gurobi.hs b/src/Numeric/Optimization/MIP/Solution/Gurobi.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solution/Gurobi.hs @@ -0,0 +1,79 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solution.Gurobi +-- Copyright : (c) Masahiro Sakai 2012,2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solution.Gurobi + ( Solution (..) + , render + , writeFile + , parse + , readFile + ) where + +import Prelude hiding (readFile, writeFile) +#if !MIN_VERSION_base(4,8,0) +import Control.Applicative +#endif +import Data.Default.Class +import Data.Interned (intern, unintern) +import Data.List (foldl') +import qualified Data.Map as Map +#if !MIN_VERSION_base(4,11,0) +import Data.Monoid +#endif +import Data.Scientific (Scientific) +import qualified Data.Text.Lazy as TL +import qualified Data.Text.Lazy.Builder as B +import qualified Data.Text.Lazy.Builder.Scientific as B +import qualified Data.Text.Lazy.IO as TLIO +import Numeric.Optimization.MIP (Solution) +import qualified Numeric.Optimization.MIP as MIP + +render :: MIP.Solution Scientific -> TL.Text +render sol = B.toLazyText $ ls1 <> mconcat ls2 + where + ls1 = case MIP.solObjectiveValue sol of + Nothing -> mempty + Just val -> "# Objective value = " <> B.scientificBuilder val <> B.singleton '\n' + ls2 = [ B.fromText (unintern name) <> B.singleton ' ' <> B.scientificBuilder val <> B.singleton '\n' + | (name,val) <- Map.toList (MIP.solVariables sol) + ] + +writeFile :: FilePath -> MIP.Solution Scientific -> IO () +writeFile fname sol = do + TLIO.writeFile fname (render sol) + +parse :: TL.Text -> MIP.Solution Scientific +parse t = + case foldl' f (Nothing,[]) $ TL.lines t of + (obj, vs) -> + def{ MIP.solStatus = MIP.StatusFeasible + , MIP.solObjectiveValue = obj + , MIP.solVariables = Map.fromList vs + } + where + f :: (Maybe Scientific, [(MIP.Var, Scientific)]) -> TL.Text -> (Maybe Scientific, [(MIP.Var, Scientific)]) + f (obj,vs) l + | Just l2 <- TL.stripPrefix "# " l + , Just l3 <- TL.stripPrefix "objective value = " (TL.toLower l2) + , (r:_) <- [r | (r,[]) <- reads (TL.unpack l3)] = + (Just r, vs) + | otherwise = + case TL.words (TL.takeWhile (/= '#') l) of + [w1, w2] -> (obj, (intern (TL.toStrict w1), read (TL.unpack w2)) : vs) + [] -> (obj, vs) + _ -> error ("Numeric.Optimization.MIP.Solution.Gurobi: invalid line " ++ show l) + +readFile :: FilePath -> IO (MIP.Solution Scientific) +readFile fname = parse <$> TLIO.readFile fname diff --git a/src/Numeric/Optimization/MIP/Solution/SCIP.hs b/src/Numeric/Optimization/MIP/Solution/SCIP.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solution/SCIP.hs @@ -0,0 +1,95 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solution.SCIP +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solution.SCIP + ( Solution (..) + , parse + , readFile + ) where + +import Prelude hiding (readFile, writeFile) +#if !MIN_VERSION_base(4,8,0) +import Control.Applicative +#endif +import Control.Monad.Except +import Data.Interned (intern) +import Data.Map (Map) +import qualified Data.Map as Map +import Data.Scientific (Scientific) +import qualified Data.Text as T +import qualified Data.Text.Lazy as TL +import qualified Data.Text.Lazy.IO as TLIO +import Numeric.Optimization.MIP (Solution) +import qualified Numeric.Optimization.MIP as MIP + +parse :: TL.Text -> MIP.Solution Scientific +parse t = + case parse' $ TL.lines t of + Left e -> error e + Right x -> x + +parse' :: [TL.Text] -> Either String (MIP.Solution Scientific) +parse' (t1:t2:ts) = do + status <- + case TL.stripPrefix "solution status:" t1 of + Nothing -> throwError "first line must start with \"solution status:\"" + Just s -> return $ Map.findWithDefault MIP.StatusUnknown (TL.toStrict $ TL.strip s) statusTable + if t2 == "no solution available" then do + return $ + MIP.Solution + { MIP.solStatus = status + , MIP.solObjectiveValue = Nothing + , MIP.solVariables = Map.empty + } + else do + obj <- + case TL.stripPrefix "objective value:" t2 of + Nothing -> throwError "second line must start with \"objective value:\"" + Just s -> return $ read $ TL.unpack $ TL.strip s + let f :: [(MIP.Var, Scientific)] -> TL.Text -> Either String [(MIP.Var, Scientific)] + f vs t = + case TL.words t of + (w1:w2:_) -> return $ (intern (TL.toStrict w1), read (TL.unpack w2)) : vs + [] -> return $ vs + _ -> throwError ("Numeric.Optimization.MIP.Solution.SCIP: invalid line " ++ show t) + vs <- foldM f [] ts + return $ + MIP.Solution + { MIP.solStatus = status + , MIP.solObjectiveValue = Just obj + , MIP.solVariables = Map.fromList vs + } +parse' _ = throwError "must have >=2 lines" + +statusTable :: Map T.Text MIP.Status +statusTable = Map.fromList + [ ("user interrupt", MIP.StatusUnknown) + , ("node limit reached", MIP.StatusUnknown) + , ("total node limit reached", MIP.StatusUnknown) + , ("stall node limit reached", MIP.StatusUnknown) + , ("time limit reached", MIP.StatusUnknown) + , ("memory limit reached", MIP.StatusUnknown) + , ("gap limit reached", MIP.StatusUnknown) + , ("solution limit reached", MIP.StatusUnknown) + , ("solution improvement limit reached", MIP.StatusUnknown) + , ("optimal solution found", MIP.StatusOptimal) + , ("infeasible", MIP.StatusInfeasible) + , ("unbounded", MIP.StatusUnbounded) + , ("infeasible or unbounded", MIP.StatusInfeasibleOrUnbounded) + -- , ("unknown", ) + ] + +readFile :: FilePath -> IO (MIP.Solution Scientific) +readFile fname = parse <$> TLIO.readFile fname diff --git a/src/Numeric/Optimization/MIP/Solver.hs b/src/Numeric/Optimization/MIP/Solver.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solver.hs @@ -0,0 +1,29 @@ +{-# OPTIONS_GHC -Wall #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solver +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solver + ( module Numeric.Optimization.MIP.Solver.Base + , module Numeric.Optimization.MIP.Solver.CBC + , module Numeric.Optimization.MIP.Solver.CPLEX + , module Numeric.Optimization.MIP.Solver.Glpsol + , module Numeric.Optimization.MIP.Solver.GurobiCl + , module Numeric.Optimization.MIP.Solver.LPSolve + , module Numeric.Optimization.MIP.Solver.SCIP + ) where + +import Numeric.Optimization.MIP.Solver.Base +import Numeric.Optimization.MIP.Solver.CBC +import Numeric.Optimization.MIP.Solver.CPLEX +import Numeric.Optimization.MIP.Solver.Glpsol +import Numeric.Optimization.MIP.Solver.GurobiCl +import Numeric.Optimization.MIP.Solver.LPSolve +import Numeric.Optimization.MIP.Solver.SCIP diff --git a/src/Numeric/Optimization/MIP/Solver/Base.hs b/src/Numeric/Optimization/MIP/Solver/Base.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solver/Base.hs @@ -0,0 +1,44 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE FunctionalDependencies #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solver.Base +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solver.Base + ( SolveOptions (..) + , Default (..) + , IsSolver (..) + ) where + +import Data.Default.Class +import Data.Scientific (Scientific) +import Numeric.Optimization.MIP.Base as MIP + +data SolveOptions + = SolveOptions + { solveTimeLimit :: Maybe Double + -- ^ time limit in seconds + , solveLogger :: String -> IO () + -- ^ invoked when a solver output a line + , solveErrorLogger :: String -> IO () + -- ^ invoked when a solver output a line to stderr + } + +instance Default SolveOptions where + def = + SolveOptions + { solveTimeLimit = Nothing + , solveLogger = const $ return () + , solveErrorLogger = const $ return () + } + +class Monad m => IsSolver s m | s -> m where + solve :: s -> SolveOptions -> MIP.Problem Scientific -> m (MIP.Solution Scientific) diff --git a/src/Numeric/Optimization/MIP/Solver/CBC.hs b/src/Numeric/Optimization/MIP/Solver/CBC.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solver/CBC.hs @@ -0,0 +1,67 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solver.CBC +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solver.CBC + ( CBC (..) + , cbc + ) where + +import Data.Default.Class +import qualified Data.Text.Lazy.IO as TLIO +import System.Exit +import System.IO +import System.IO.Temp +import qualified Numeric.Optimization.MIP.Base as MIP +import qualified Numeric.Optimization.MIP.LPFile as LPFile +import Numeric.Optimization.MIP.Solver.Base +import qualified Numeric.Optimization.MIP.Solution.CBC as CBCSol +import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback) + +data CBC + = CBC + { cbcPath :: String + } + +instance Default CBC where + def = cbc + +cbc :: CBC +cbc = CBC "cbc" + +instance IsSolver CBC IO where + solve solver opt prob = do + case LPFile.render def prob of + Left err -> ioError $ userError err + Right lp -> do + withSystemTempFile "cbc.lp" $ \fname1 h1 -> do + TLIO.hPutStr h1 lp + hClose h1 + withSystemTempFile "cbc.sol" $ \fname2 h2 -> do + hClose h2 + let args = [fname1] + ++ (case solveTimeLimit opt of + Nothing -> [] + Just sec -> ["sec", show sec]) + ++ ["solve", "solu", fname2] + onGetLine = solveLogger opt + onGetErrorLine = solveErrorLogger opt + exitcode <- runProcessWithOutputCallback (cbcPath solver) args "" onGetLine onGetErrorLine + case exitcode of + ExitFailure n -> ioError $ userError $ "exit with " ++ show n + ExitSuccess -> do + sol <- CBCSol.readFile fname2 + if MIP.objDir (MIP.objectiveFunction prob) == MIP.OptMax then + return $ sol{ MIP.solObjectiveValue = fmap negate (MIP.solObjectiveValue sol) } + else + return sol diff --git a/src/Numeric/Optimization/MIP/Solver/CPLEX.hs b/src/Numeric/Optimization/MIP/Solver/CPLEX.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solver/CPLEX.hs @@ -0,0 +1,83 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solver.CPLEX +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solver.CPLEX + ( CPLEX (..) + , cplex + ) where + +import Control.Monad +import Data.Default.Class +import Data.IORef +import qualified Data.Text.Lazy.IO as TLIO +import System.Exit +import System.IO +import System.IO.Temp +import qualified Numeric.Optimization.MIP.Base as MIP +import qualified Numeric.Optimization.MIP.LPFile as LPFile +import Numeric.Optimization.MIP.Solver.Base +import qualified Numeric.Optimization.MIP.Solution.CPLEX as CPLEXSol +import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback) + +data CPLEX + = CPLEX + { cplexPath :: String + } + +instance Default CPLEX where + def = cplex + +cplex :: CPLEX +cplex = CPLEX "cplex" + +instance IsSolver CPLEX IO where + solve solver opt prob = do + case LPFile.render def prob of + Left err -> ioError $ userError err + Right lp -> do + withSystemTempFile "cplex.lp" $ \fname1 h1 -> do + TLIO.hPutStr h1 lp + hClose h1 + withSystemTempFile "cplex.sol" $ \fname2 h2 -> do + hClose h2 + isInfeasibleRef <- newIORef False + let args = [] + input = unlines $ + (case solveTimeLimit opt of + Nothing -> [] + Just sec -> ["set timelimit ", show sec]) ++ + [ "read " ++ show fname1 + , "optimize" + , "write " ++ show fname2 + , "y" + , "quit" + ] + onGetLine s = do + when (s == "MIP - Integer infeasible.") $ do + writeIORef isInfeasibleRef True + solveLogger opt s + onGetErrorLine = solveErrorLogger opt + exitcode <- runProcessWithOutputCallback (cplexPath solver) args input onGetLine onGetErrorLine + case exitcode of + ExitFailure n -> ioError $ userError $ "exit with " ++ show n + ExitSuccess -> do + size <- withFile fname2 ReadMode $ hFileSize + if size == 0 then do + isInfeasible <- readIORef isInfeasibleRef + if isInfeasible then + return def{ MIP.solStatus = MIP.StatusInfeasible } + else + return def + else + CPLEXSol.readFile fname2 diff --git a/src/Numeric/Optimization/MIP/Solver/Glpsol.hs b/src/Numeric/Optimization/MIP/Solver/Glpsol.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solver/Glpsol.hs @@ -0,0 +1,83 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solver.Glpsol +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solver.Glpsol + ( Glpsol (..) + , glpsol + ) where + +import Algebra.PartialOrd +import Data.Default.Class +import Data.IORef +import qualified Data.Text.Lazy.IO as TLIO +import System.Exit +import System.IO +import System.IO.Temp +import qualified Numeric.Optimization.MIP.Base as MIP +import qualified Numeric.Optimization.MIP.LPFile as LPFile +import Numeric.Optimization.MIP.Solver.Base +import qualified Numeric.Optimization.MIP.Solution.GLPK as GLPKSol +import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback) + +data Glpsol + = Glpsol + { glpsolPath :: String + } + +instance Default Glpsol where + def = glpsol + +glpsol :: Glpsol +glpsol = Glpsol "glpsol" + +instance IsSolver Glpsol IO where + solve solver opt prob = do + case LPFile.render def prob of + Left err -> ioError $ userError err + Right lp -> do + withSystemTempFile "glpsol.lp" $ \fname1 h1 -> do + TLIO.hPutStr h1 lp + hClose h1 + withSystemTempFile "glpsol.sol" $ \fname2 h2 -> do + hClose h2 + isUnboundedRef <- newIORef False + isInfeasibleRef <- newIORef False + let args = ["--lp", fname1, "-o", fname2] ++ + (case solveTimeLimit opt of + Nothing -> [] + Just sec -> ["--tmlim", show (max 1 (floor sec) :: Int)]) + onGetLine s = do + case s of + "LP HAS UNBOUNDED PRIMAL SOLUTION" -> + writeIORef isUnboundedRef True + "PROBLEM HAS UNBOUNDED SOLUTION" -> + writeIORef isUnboundedRef True + "PROBLEM HAS NO PRIMAL FEASIBLE SOLUTION" -> + writeIORef isInfeasibleRef True + _ -> return () + solveLogger opt s + onGetErrorLine = solveErrorLogger opt + exitcode <- runProcessWithOutputCallback (glpsolPath solver) args "" onGetLine onGetErrorLine + case exitcode of + ExitFailure n -> ioError $ userError $ "exit with " ++ show n + ExitSuccess -> do + sol <- GLPKSol.readFile fname2 + isUnbounded <- readIORef isUnboundedRef + isInfeasible <- readIORef isInfeasibleRef + if isUnbounded && MIP.solStatus sol `leq` MIP.StatusInfeasibleOrUnbounded then + return $ sol{ MIP.solStatus = MIP.StatusInfeasibleOrUnbounded } + else if isInfeasible && MIP.solStatus sol `leq` MIP.StatusInfeasible then + return $ sol{ MIP.solStatus = MIP.StatusInfeasible } + else + return sol diff --git a/src/Numeric/Optimization/MIP/Solver/GurobiCl.hs b/src/Numeric/Optimization/MIP/Solver/GurobiCl.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solver/GurobiCl.hs @@ -0,0 +1,76 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solver.GurobiCl +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solver.GurobiCl + ( GurobiCl (..) + , gurobiCl + ) where + +import Data.Default.Class +import Data.IORef +import Data.List (isPrefixOf) +import qualified Data.Text.Lazy.IO as TLIO +import System.Exit +import System.IO +import System.IO.Temp +import qualified Numeric.Optimization.MIP.Base as MIP +import qualified Numeric.Optimization.MIP.LPFile as LPFile +import Numeric.Optimization.MIP.Solver.Base +import qualified Numeric.Optimization.MIP.Solution.Gurobi as GurobiSol +import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback) + +data GurobiCl + = GurobiCl + { gurobiClPath :: String + } + +instance Default GurobiCl where + def = gurobiCl + +gurobiCl :: GurobiCl +gurobiCl = GurobiCl "gurobi_cl" + +instance IsSolver GurobiCl IO where + solve solver opt prob = do + case LPFile.render def prob of + Left err -> ioError $ userError err + Right lp -> do + withSystemTempFile "gurobi.lp" $ \fname1 h1 -> do + TLIO.hPutStr h1 lp + hClose h1 + withSystemTempFile "gurobi.sol" $ \fname2 h2 -> do + hClose h2 + statusRef <- newIORef MIP.StatusUnknown + let args = ["ResultFile=" ++ fname2] + ++ (case solveTimeLimit opt of + Nothing -> [] + Just sec -> ["TimeLimit=" ++ show sec]) + ++ [fname1] + onGetLine s = do + case s of + -- "Time limit reached" -> writeIORef statusRef MIP.StatusUnknown + "Model is unbounded" -> writeIORef statusRef MIP.StatusUnbounded + "Model is infeasible" -> writeIORef statusRef MIP.StatusInfeasible + "Model is infeasible or unbounded" -> writeIORef statusRef MIP.StatusInfeasibleOrUnbounded + _ | isPrefixOf "Optimal solution found" s -> writeIORef statusRef MIP.StatusOptimal + _ -> return () + solveLogger opt s + onGetErrorLine = solveErrorLogger opt + exitcode <- runProcessWithOutputCallback (gurobiClPath solver) args "" onGetLine onGetErrorLine + case exitcode of + ExitFailure n -> ioError $ userError $ "exit with " ++ show n + ExitSuccess -> do + status <- readIORef statusRef + sol <- GurobiSol.readFile fname2 + return $ sol{ MIP.solStatus = status } diff --git a/src/Numeric/Optimization/MIP/Solver/LPSolve.hs b/src/Numeric/Optimization/MIP/Solver/LPSolve.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solver/LPSolve.hs @@ -0,0 +1,97 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solver.LPSolve +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solver.LPSolve + ( LPSolve (..) + , lpSolve + ) where + +import Control.Monad +import Data.Default.Class +import Data.IORef +import Data.List (stripPrefix) +import qualified Data.Map as Map +import Data.String +import qualified Data.Text.Lazy.IO as TLIO +import System.Exit +import System.IO +import System.IO.Temp +import qualified Numeric.Optimization.MIP.Base as MIP +import qualified Numeric.Optimization.MIP.MPSFile as MPSFile +import Numeric.Optimization.MIP.Solver.Base +import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback) + +data LPSolve + = LPSolve + { lpSolvePath :: String + } + +instance Default LPSolve where + def = lpSolve + +lpSolve :: LPSolve +lpSolve = LPSolve "lp_solve" + +instance IsSolver LPSolve IO where + solve solver opt prob = do + case MPSFile.render def prob of + Left err -> ioError $ userError err + Right lp -> do + withSystemTempFile "lp_solve.mps" $ \fname1 h1 -> do + TLIO.hPutStr h1 lp + hClose h1 + objRef <- newIORef Nothing + solRef <- newIORef [] + flagRef <- newIORef False + let args = (case solveTimeLimit opt of + Nothing -> [] + Just sec -> ["-timeout", show sec]) + ++ ["-fmps", fname1] + onGetLine s = do + case s of + "Actual values of the variables:" -> writeIORef flagRef True + _ | Just v <- stripPrefix "Value of objective function: " s -> do + writeIORef objRef (Just (read v)) + _ -> do + flag <- readIORef flagRef + when flag $ do + case words s of + [var,val] -> modifyIORef solRef ((fromString var, read val) :) + _ -> return () + return () + solveLogger opt s + onGetErrorLine = solveErrorLogger opt + exitcode <- runProcessWithOutputCallback (lpSolvePath solver) args "" onGetLine onGetErrorLine + status <- + case exitcode of + ExitSuccess -> return MIP.StatusOptimal + ExitFailure (-2) -> return MIP.StatusUnknown -- NOMEMORY + ExitFailure 1 -> return MIP.StatusFeasible -- SUBOPTIMAL + ExitFailure 2 -> return MIP.StatusInfeasible -- INFEASIBLE + ExitFailure 3 -> return MIP.StatusInfeasibleOrUnbounded -- UNBOUNDED + ExitFailure 4 -> return MIP.StatusUnknown -- DEGENERATE + ExitFailure 5 -> return MIP.StatusUnknown -- NUMFAILURE + ExitFailure 6 -> return MIP.StatusUnknown -- USERABORT + ExitFailure 7 -> return MIP.StatusUnknown -- TIMEOUT + ExitFailure 9 -> return MIP.StatusOptimal -- PRESOLVED + ExitFailure 25 -> return MIP.StatusUnknown -- NUMFAILURE + ExitFailure n -> ioError $ userError $ "unknown exit code: " ++ show n + obj <- readIORef objRef + sol <- readIORef solRef + return $ + MIP.Solution + { MIP.solStatus = status + , MIP.solObjectiveValue = obj + , MIP.solVariables = Map.fromList sol + } diff --git a/src/Numeric/Optimization/MIP/Solver/SCIP.hs b/src/Numeric/Optimization/MIP/Solver/SCIP.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Optimization/MIP/Solver/SCIP.hs @@ -0,0 +1,64 @@ +{-# OPTIONS_GHC -Wall #-} +{-# OPTIONS_HADDOCK show-extensions #-} +{-# LANGUAGE MultiParamTypeClasses #-} +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Optimization.MIP.Solver.SCIP +-- Copyright : (c) Masahiro Sakai 2017 +-- License : BSD-style +-- +-- Maintainer : masahiro.sakai@gmail.com +-- Stability : provisional +-- Portability : non-portable +-- +----------------------------------------------------------------------------- +module Numeric.Optimization.MIP.Solver.SCIP + ( SCIP (..) + , scip + ) where + +import Data.Default.Class +import qualified Data.Text.Lazy.IO as TLIO +import System.Exit +import System.IO +import System.IO.Temp +import qualified Numeric.Optimization.MIP.LPFile as LPFile +import Numeric.Optimization.MIP.Solver.Base +import qualified Numeric.Optimization.MIP.Solution.SCIP as ScipSol +import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback) + +data SCIP + = SCIP + { scipPath :: String + } + +instance Default SCIP where + def = scip + +scip :: SCIP +scip = SCIP "scip" + +instance IsSolver SCIP IO where + solve solver opt prob = do + case LPFile.render def prob of + Left err -> ioError $ userError err + Right lp -> do + withSystemTempFile "scip.lp" $ \fname1 h1 -> do + TLIO.hPutStr h1 lp + hClose h1 + withSystemTempFile "scip.sol" $ \fname2 h2 -> do + hClose h2 + let args = [ "-c", "read " ++ show fname1 ] + ++ (case solveTimeLimit opt of + Nothing -> [] + Just sec -> ["-c", "set limits time " ++ show sec]) + ++ [ "-c", "optimize" + , "-c", "write solution " ++ show fname2 + , "-c", "quit" + ] + onGetLine = solveLogger opt + onGetErrorLine = solveErrorLogger opt + exitcode <- runProcessWithOutputCallback (scipPath solver) args "" onGetLine onGetErrorLine + case exitcode of + ExitFailure n -> ioError $ userError $ "exit with " ++ show n + ExitSuccess -> ScipSol.readFile fname2 diff --git a/test/Test/LPFile.hs b/test/Test/LPFile.hs new file mode 100644 --- /dev/null +++ b/test/Test/LPFile.hs @@ -0,0 +1,64 @@ +{-# LANGUAGE TemplateHaskell #-} +module Test.LPFile (lpTestGroup) where + +import Control.Monad +import Data.Default.Class +import Data.List +import Data.Maybe +import Test.Tasty +import Test.Tasty.QuickCheck +import Test.Tasty.HUnit +import Test.Tasty.TH +import Numeric.Optimization.MIP.LPFile + +case_testdata = checkString "testdata" testdata +case_test_indicator = checkFile "samples/lp/test-indicator.lp" +case_test_qcp = checkFile "samples/lp/test-qcp.lp" +case_test_qcp2 = checkFile "samples/lp/test-qcp2.lp" +case_test_qp = checkFile "samples/lp/test-qp.lp" +case_empty_obj_1 = checkFile "samples/lp/empty_obj_1.lp" +case_empty_obj_2 = checkFile "samples/lp/empty_obj_2.lp" + +------------------------------------------------------------------------ +-- Sample data + +testdata :: String +testdata = unlines + [ "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" + ] + +------------------------------------------------------------------------ +-- Utilities + +checkFile :: FilePath -> Assertion +checkFile fname = do + lp <- parseFile def fname + case render def lp of + Left err -> assertFailure ("render failure: " ++ err) + Right _ -> return () + +checkString :: String -> String -> Assertion +checkString name str = do + case parseString def name str of + Left err -> assertFailure $ show err + Right lp -> + case render def lp of + Left err -> assertFailure ("render failure: " ++ err) + Right _ -> return () + +------------------------------------------------------------------------ +-- Test harness + +lpTestGroup :: TestTree +lpTestGroup = $(testGroupGenerator) diff --git a/test/Test/MIP.hs b/test/Test/MIP.hs new file mode 100644 --- /dev/null +++ b/test/Test/MIP.hs @@ -0,0 +1,154 @@ +{-# OPTIONS_GHC -Wall -fno-warn-orphans #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TemplateHaskell #-} +module Test.MIP (mipTestGroup) where + +import Algebra.PartialOrd +#if !MIN_VERSION_lattices(2,0,0) +import Algebra.Lattice +#endif +import Data.Maybe +import qualified Data.Map as Map +import Test.Tasty +import Test.Tasty.HUnit +import Test.Tasty.QuickCheck +import Test.Tasty.TH +import Numeric.Optimization.MIP (meetStatus) +import qualified Numeric.Optimization.MIP as MIP +import qualified Numeric.Optimization.MIP.Solution.CBC as CBCSol +import qualified Numeric.Optimization.MIP.Solution.CPLEX as CPLEXSol +import qualified Numeric.Optimization.MIP.Solution.GLPK as GLPKSol +import qualified Numeric.Optimization.MIP.Solution.Gurobi as GurobiSol +import qualified Numeric.Optimization.MIP.Solution.SCIP as SCIPSol + +prop_status_refl :: Property +prop_status_refl = forAll arbitrary $ \(x :: MIP.Status) -> do + x `leq` x + +prop_status_trans :: Property +prop_status_trans = + forAll arbitrary $ \(x :: MIP.Status) -> + forAll (upper x) $ \y -> + forAll (upper y) $ \z -> + x `leq` z + where + -- upper :: (PartialOrd a, Enum a, Bounded a) => a -> Gen a + upper a = elements [b | b <- [minBound .. maxBound], a `leq` b] + +prop_status_meet_idempotency :: Property +prop_status_meet_idempotency = + forAll arbitrary $ \(x :: MIP.Status) -> + x `meetStatus` x == x + +prop_status_meet_comm :: Property +prop_status_meet_comm = + forAll arbitrary $ \(x :: MIP.Status) y -> + x `meetStatus` y == y `meetStatus` x + +prop_status_meet_assoc :: Property +prop_status_meet_assoc = + forAll arbitrary $ \(x :: MIP.Status) y z -> + (x `meetStatus` y) `meetStatus` z == x `meetStatus` (y `meetStatus` z) + +prop_status_meet_leq :: Property +prop_status_meet_leq = + forAll arbitrary $ \(x :: MIP.Status) y -> +#if MIN_VERSION_lattices(2,0,0) + (x == (x `meetStatus` y)) == x `leq` y +#else + x `meetLeq` y == x `leq` y +#endif + +instance Arbitrary MIP.Status where + arbitrary = arbitraryBoundedEnum + +case_CBCSol :: Assertion +case_CBCSol = do + sol <- CBCSol.readFile "samples/lp/test-solution-cbc.txt" + 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_CBCSol_infeasible :: Assertion +case_CBCSol_infeasible = do + sol <- CBCSol.readFile "samples/lp/test-solution-cbc-infeasible.txt" + sol @?= + MIP.Solution + { MIP.solStatus = MIP.StatusInfeasible + , MIP.solObjectiveValue = Just 0.00000000 + , MIP.solVariables = Map.fromList [("x", 0.11111111), ("y", 0), ("z", 0.33333333)] + } + +case_CBCSol_unbounded :: Assertion +case_CBCSol_unbounded = do + sol <- CBCSol.readFile "samples/lp/test-solution-cbc-unbounded.txt" + sol @?= + MIP.Solution + { MIP.solStatus = MIP.StatusInfeasibleOrUnbounded + , MIP.solObjectiveValue = Just 0.00000000 + , MIP.solVariables = Map.fromList [("x", 0), ("y", 0)] + } + +case_CPLEXSol :: Assertion +case_CPLEXSol = do + sol <- CPLEXSol.readFile "samples/lp/test-solution-cplex.sol" + 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_CPLEXSol_unbounded :: Assertion +case_CPLEXSol_unbounded = do + sol <- CPLEXSol.readFile "samples/lp/test-solution-cplex-unbounded.sol" + sol @?= + MIP.Solution + { MIP.solStatus = MIP.StatusUnbounded + , MIP.solObjectiveValue = Just 3.0 + , MIP.solVariables = Map.fromList [("x", 1.0), ("y", 2.0)] + } + +case_GLPKSol :: Assertion +case_GLPKSol = do + sol <- GLPKSol.readFile "samples/lp/test-solution-glpk.sol" + 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_GLPKSol_long_var :: Assertion +case_GLPKSol_long_var = do + sol <- GLPKSol.readFile "samples/lp/test-solution-glpk-long.sol" + sol @?= + MIP.Solution + { MIP.solStatus = MIP.StatusOptimal + , MIP.solObjectiveValue = Just 122.5 + , MIP.solVariables = Map.fromList [("x1AAAAAAAAAAAAAAAAAAAAAAAAAAAA", 40), ("x2", 10.5), ("x3", 19.5), ("x4", 3)] + } + +case_GurobiSol :: Assertion +case_GurobiSol = do + sol <- GurobiSol.readFile "samples/lp/test-solution-gurobi.sol" + isJust (GurobiSol.solObjectiveValue sol) @?= True + GurobiSol.parse (GurobiSol.render sol) @?= sol + +case_SCIPSol :: Assertion +case_SCIPSol = do + sol <- SCIPSol.readFile "samples/lp/test-solution-scip.sol" + 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)] + } + +mipTestGroup :: TestTree +mipTestGroup = $(testGroupGenerator) diff --git a/test/Test/MIPSolver.hs b/test/Test/MIPSolver.hs new file mode 100644 --- /dev/null +++ b/test/Test/MIPSolver.hs @@ -0,0 +1,301 @@ +{-# OPTIONS_GHC -Wall -Wno-unused-top-binds #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE OverloadedStrings #-} +{-# LANGUAGE TemplateHaskell #-} +module Test.MIPSolver (mipSolverTestGroup) 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 + +-- ------------------------------------------------------------------------ + +case_cbc :: Assertion +case_cbc = do + prob <- MIP.readFile def "samples/lp/test.lp" + sol <- solve cbc 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_cbc_unbounded :: Assertion +case_cbc_unbounded = do + prob <- MIP.readFile def "samples/lp/unbounded-ip.lp" + sol <- solve cbc 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_cbc_infeasible :: Assertion +case_cbc_infeasible = do + prob <- MIP.readFile def "samples/lp/infeasible.lp" + sol <- solve cbc def prob + MIP.solStatus sol @?= MIP.StatusInfeasible + +case_cbc_infeasible2 :: Assertion +case_cbc_infeasible2 = do + prob <- MIP.readFile def "samples/lp/glpk-preprocess-bug.lp" + sol <- solve cbc def prob + let status = MIP.solStatus sol + unless (status == MIP.StatusInfeasible || status == MIP.StatusInfeasibleOrUnbounded) $ + assertFailure $ unlines $ + [ "expected: StatusInfeasible or StatusInfeasibleOrUnbounded" + , " but got: " ++ show status + ] + +-- ------------------------------------------------------------------------ + +case_cplex :: Assertion +case_cplex = do + prob <- MIP.readFile def "samples/lp/test.lp" + sol <- solve cplex 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_cplex_unbounded :: Assertion +case_cplex_unbounded = do + prob <- MIP.readFile def "samples/lp/unbounded-ip.lp" + sol <- solve cplex 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_cplex_infeasible :: Assertion +case_cplex_infeasible = do + prob <- MIP.readFile def "samples/lp/infeasible.lp" + sol <- solve cplex def prob + MIP.solStatus sol @?= MIP.StatusInfeasible + +case_cplex_infeasible2 :: Assertion +case_cplex_infeasible2 = do + prob <- MIP.readFile def "samples/lp/glpk-preprocess-bug.lp" + sol <- solve cplex def prob + let status = MIP.solStatus sol + unless (status == MIP.StatusInfeasible || status == MIP.StatusInfeasibleOrUnbounded) $ + assertFailure $ unlines $ + [ "expected: StatusInfeasible or StatusInfeasibleOrUnbounded" + , " but got: " ++ show status + ] + +-- ------------------------------------------------------------------------ + +case_glpsol :: Assertion +case_glpsol = do + prob <- MIP.readFile def "samples/lp/test.lp" + sol <- solve glpsol 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_glpsol_unbounded :: Assertion +case_glpsol_unbounded = do + prob <- MIP.readFile def "samples/lp/unbounded-ip.lp" + sol <- solve glpsol 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_glpsol_infeasible :: Assertion +case_glpsol_infeasible = do + prob <- MIP.readFile def "samples/lp/infeasible.lp" + sol <- solve glpsol def prob + MIP.solStatus sol @?= MIP.StatusInfeasible + +-- ------------------------------------------------------------------------ + +case_gurobiCl :: Assertion +case_gurobiCl = do + prob <- MIP.readFile def "samples/lp/test.lp" + sol <- solve gurobiCl def prob + sol @?= + MIP.Solution + { MIP.solStatus = MIP.StatusOptimal + , MIP.solObjectiveValue = Just 122.50000000000006 + , MIP.solVariables = Map.fromList [("x1", 40), ("x2", 10.5), ("x3", 19.500000000000018), ("x4", 3)] + } + +case_gurobiCl_unbounded :: Assertion +case_gurobiCl_unbounded = do + prob <- MIP.readFile def "samples/lp/unbounded-ip.lp" + sol <- solve gurobiCl 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_gurobiCl_infeasible :: Assertion +case_gurobiCl_infeasible = do + prob <- MIP.readFile def "samples/lp/infeasible.lp" + sol <- solve gurobiCl def prob + MIP.solStatus sol @?= MIP.StatusInfeasible + +case_gurobiCl_infeasible2 :: Assertion +case_gurobiCl_infeasible2 = do + prob <- MIP.readFile def "samples/lp/glpk-preprocess-bug.lp" + sol <- solve gurobiCl def prob + let status = MIP.solStatus sol + unless (status == MIP.StatusInfeasible || status == MIP.StatusInfeasibleOrUnbounded) $ + assertFailure $ unlines $ + [ "expected: StatusInfeasible or StatusInfeasibleOrUnbounded" + , " but got: " ++ show status + ] + +-- ------------------------------------------------------------------------ + +case_lpSolve :: Assertion +case_lpSolve = do + prob <- MIP.readFile def "samples/lp/test.lp" + sol <- solve lpSolve 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_lpSolve_unbounded :: Assertion +case_lpSolve_unbounded = do + prob <- MIP.readFile def "samples/lp/unbounded-ip.lp" + sol <- solve lpSolve 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_lpSolve_infeasible :: Assertion +case_lpSolve_infeasible = do + prob <- MIP.readFile def "samples/lp/infeasible.lp" + sol <- solve lpSolve def prob + MIP.solStatus sol @?= MIP.StatusInfeasible + +case_lpSolve_infeasible2 :: Assertion +case_lpSolve_infeasible2 = do + prob <- MIP.readFile def "samples/lp/glpk-preprocess-bug.lp" + sol <- solve lpSolve def prob + let status = MIP.solStatus sol + unless (status == MIP.StatusInfeasible || status == MIP.StatusInfeasibleOrUnbounded) $ + assertFailure $ unlines $ + [ "expected: StatusInfeasible or StatusInfeasibleOrUnbounded" + , " but got: " ++ show status + ] + +-- ------------------------------------------------------------------------ + +case_scip :: Assertion +case_scip = do + prob <- MIP.readFile def "samples/lp/test.lp" + sol <- solve scip 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_scip_unbounded :: Assertion +case_scip_unbounded = do + prob <- MIP.readFile def "samples/lp/unbounded-ip.lp" + sol <- solve scip 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_scip_infeasible :: Assertion +case_scip_infeasible = do + prob <- MIP.readFile def "samples/lp/infeasible.lp" + sol <- solve scip def prob + MIP.solStatus sol @?= MIP.StatusInfeasible + +case_scip_infeasible2 :: Assertion +case_scip_infeasible2 = do + prob <- MIP.readFile def "samples/lp/glpk-preprocess-bug.lp" + sol <- solve scip def prob + let status = MIP.solStatus sol + unless (status == MIP.StatusInfeasible || status == MIP.StatusInfeasibleOrUnbounded) $ + assertFailure $ unlines $ + [ "expected: StatusInfeasible or StatusInfeasibleOrUnbounded" + , " but got: " ++ show status + ] + +-- ------------------------------------------------------------------------ + +mipSolverTestGroup :: TestTree +mipSolverTestGroup = testGroup "Test.MIPSolver" $ [] +#ifdef TEST_CBC + ++ + [ testCase "cbc" case_cbc + , testCase "cbc unbounded" case_cbc_unbounded + , testCase "cbc infeasible" case_cbc_infeasible + , testCase "cbc infeasible2" case_cbc_infeasible2 + ] +#endif +#ifdef TEST_CPLEX + ++ + [ testCase "cplex" case_cplex + , testCase "cplex unbounded" case_cplex_unbounded + , testCase "cplex infeasible" case_cplex_infeasible + , testCase "cplex infeasible2" case_cplex_infeasible2 + ] +#endif +#ifdef TEST_GLPSOL + ++ + [ testCase "glpsol" case_glpsol + , testCase "glpsol unbounded" case_glpsol_unbounded + , testCase "glpsol infeasible" case_glpsol_infeasible + ] +#endif +#ifdef TEST_GUROBI_CL + ++ + [ testCase "gurobiCl" case_gurobiCl + , testCase "gurobiCl unbounded" case_gurobiCl_unbounded + , testCase "gurobiCl infeasible" case_gurobiCl_infeasible + , testCase "gurobiCl infeasible2" case_gurobiCl_infeasible2 + ] +#endif +#ifdef TEST_LP_SOLVE + ++ + [ testCase "lpSolve" case_lpSolve + , testCase "lpSolve unbounded" case_lpSolve_unbounded + , testCase "lpSolve infeasible" case_lpSolve_infeasible + , testCase "lpSolve infeasible2" case_lpSolve_infeasible2 + ] +#endif +#ifdef TEST_SCIP + ++ + [ testCase "scip" case_scip + , testCase "scip unbounded" case_scip_unbounded + , testCase "scip infeasible" case_scip_infeasible + , testCase "scip infeasible2" case_scip_infeasible2 + ] +#endif + diff --git a/test/Test/MPSFile.hs b/test/Test/MPSFile.hs new file mode 100644 --- /dev/null +++ b/test/Test/MPSFile.hs @@ -0,0 +1,67 @@ +{-# LANGUAGE TemplateHaskell #-} +module Test.MPSFile (mpsTestGroup) where + +import Control.Monad +import Data.Default.Class +import Data.List +import Data.Maybe +import Test.Tasty +import Test.Tasty.QuickCheck +import Test.Tasty.HUnit +import Test.Tasty.TH +import Numeric.Optimization.MIP.MPSFile + +case_testdata = checkString "testdata" testdata +case_example2 = checkFile "samples/mps/example2.mps" +case_ind1 = checkFile "samples/mps/ind1.mps" +case_intvar1 = checkFile "samples/mps/intvar1.mps" +case_intvar2 = checkFile "samples/mps/intvar2.mps" +case_quadobj1 = checkFile "samples/mps/quadobj1.mps" +case_quadobj2 = checkFile "samples/mps/quadobj2.mps" +case_ranges = checkFile "samples/mps/ranges.mps" +case_sos = checkFile "samples/mps/sos.mps" +case_sc = checkFile "samples/mps/sc.mps" + +------------------------------------------------------------------------ +-- Sample data + +testdata :: String +testdata = unlines + [ "NAME example2.mps" + , "ROWS" + , " N obj " + , " L c1 " + , " L c2 " + , "COLUMNS" + , " x1 obj -1 c1 -1" + , " x1 c2 1" + , " x2 obj -2 c1 1" + , " x2 c2 -3" + , " x3 obj -3 c1 1" + , " x3 c2 1" + , "RHS" + , " rhs c1 20 c2 30" + , "BOUNDS" + , " UP BOUND x1 40" + , "ENDATA" + ] + +------------------------------------------------------------------------ +-- Utilities + +checkFile :: FilePath -> Assertion +checkFile fname = do + _ <- parseFile def fname + return () + +checkString :: String -> String -> Assertion +checkString name str = do + case parseString def name str of + Left err -> assertFailure (show err) + Right lp -> return () + +------------------------------------------------------------------------ +-- Test harness + +mpsTestGroup :: TestTree +mpsTestGroup = $(testGroupGenerator) diff --git a/test/TestSuite.hs b/test/TestSuite.hs new file mode 100644 --- /dev/null +++ b/test/TestSuite.hs @@ -0,0 +1,16 @@ +module Main where + +import Test.Tasty (defaultMain, testGroup) + +import Test.LPFile +import Test.MIP +import Test.MIPSolver +import Test.MPSFile + +main :: IO () +main = defaultMain $ testGroup "ToySolver test suite" + [ lpTestGroup + , mipTestGroup + , mipSolverTestGroup + , mpsTestGroup + ]