diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -2,15 +2,38 @@
 
 ## Unreleased changes
 
-## 0.1.2.0
+## 0.2.0.0 (2025-02-03)
 
-* support mtl-2.3
+* `Problem` type
+  * Merge `varType` and `varBounds` fields into `varDomains` field 
+  * Rename `varType` into `varTypes` and deprecate `varType`
+  * Represent `Expr c` as `Seq (Term c)` instead of `[Term c]`
+  * Change `Var` into `newtype`
+* Add `Tol` data type and `Eval` type class
+* Add `continuousVariables :: Problem c -> Set Var`
+* Add `binaryVariables :: Problem c -> Set Var`
+* File I/O
+  * Allow CRLF in parsers
+  * Fix to print the right-hand side of indicators as integers when possible in LP files
+  * Accept numbers other than 0 or 1 as the right-hand side of indicators in LP files
+  * Add `optMPSWriteObjSense` to `FileOptions`
+  * Add `optMPSWriteObjName` to `FileOptions`
+* Solver supports
+  * Add omitted zeroes to solutions (thanks to @dpvanbalen)
+  * Allow passing arguments solver executables
+  * Support printemps solver
+  * Support HiGHS solver
+  * Add `solveTol` field to `SolveOptions`
+* Dependencies
+  * Require `base >=4.12` (i.e. GHC `>=8.6`)
+  * Require `megaparsec >=7`
+  * Support `mtl-2.3`
 
 ## 0.1.1.0
 
-* re-export `Default` class from `Numeric.Optimization.MIP`.
-* fix to work with recent version of CBC
-* generalize type of `parseLPString` and `parseMPSString`
+* re-export the `Default` class from `Numeric.Optimization.MIP`.
+* fix to work with recent versions of CBC
+* generalizethe  type of `parseLPString` and `parseMPSString`
 
 ## 0.1.0.0
 
diff --git a/MIP.cabal b/MIP.cabal
--- a/MIP.cabal
+++ b/MIP.cabal
@@ -1,11 +1,11 @@
 cabal-version: 1.18
 
--- This file has been generated from package.yaml by hpack version 0.36.0.
+-- This file has been generated from package.yaml by hpack version 0.37.0.
 --
 -- see: https://github.com/sol/hpack
 
 name:           MIP
-version:        0.1.2.0
+version:        0.2.0.0
 synopsis:       Library for using Mixed Integer Programming (MIP)
 description:    Please see the README on GitHub at <https://github.com/msakai/haskell-MIP/tree/master/MIP#readme>
 category:       Math, Algorithms, Optimisation, Optimization
@@ -21,6 +21,48 @@
     README.md
     ChangeLog.md
     doc-images/MIP-Status-diagram.tex
+    samples/mps/example2.mps
+    samples/mps/ind1.mps
+    samples/mps/intvar1.mps
+    samples/mps/intvar2.mps
+    samples/mps/quadobj1.mps
+    samples/mps/quadobj2.mps
+    samples/mps/ranges.mps
+    samples/mps/sc.mps
+    samples/mps/sos.mps
+    samples/mps/test-bounds-bv.mps
+    samples/mps/test-bounds-fixed.mps
+    samples/mps/test-bounds-free.mps
+    samples/mps/test-bounds-mi.mps
+    samples/mps/test-bounds-pl.mps
+    samples/mps/test-implicit-bv.mps
+    samples/mps/test-lazy-constraints.mps
+    samples/mps/test-negative-upper-bound.mps
+    samples/mps/test-qcp.mps
+    samples/mps/test-semicont.mps
+    samples/mps/test-semiint-sc.mps
+    samples/mps/test-semiint-si.mps
+    samples/mps/test-sos2.mps
+    samples/mps/test-user-cuts.mps
+    samples/lp/empty_obj_1.lp
+    samples/lp/empty_obj_2.lp
+    samples/lp/glpk-preprocess-bug.lp
+    samples/lp/infeasible.lp
+    samples/lp/test-bounds-fixed.lp
+    samples/lp/test-bounds-free.lp
+    samples/lp/test-indicator.lp
+    samples/lp/test-lazy-constraints.lp
+    samples/lp/test-problem-name.lp
+    samples/lp/test-qcp.lp
+    samples/lp/test-qcp2.lp
+    samples/lp/test-qp.lp
+    samples/lp/test-semicont.lp
+    samples/lp/test-semiint.lp
+    samples/lp/test-sos.lp
+    samples/lp/test-sos2.lp
+    samples/lp/test-user-cuts.lp
+    samples/lp/test.lp
+    samples/lp/unbounded-ip.lp
 extra-doc-files:
     doc-images/MIP-Status-diagram.png
 
@@ -48,11 +90,21 @@
   manual: True
   default: False
 
+flag TestHiGHS
+  description: run test cases that depend on highs command
+  manual: True
+  default: False
+
 flag TestLPSolve
   description: run test cases that depend on lp_solve command
   manual: True
   default: False
 
+flag TestPrintemps
+  description: run test cases that depend on mps_solver.exe command of printemps
+  manual: True
+  default: False
+
 flag TestSCIP
   description: run test cases that depend on scip command
   manual: True
@@ -74,6 +126,8 @@
       Numeric.Optimization.MIP.Solution.CPLEX
       Numeric.Optimization.MIP.Solution.GLPK
       Numeric.Optimization.MIP.Solution.Gurobi
+      Numeric.Optimization.MIP.Solution.HiGHS
+      Numeric.Optimization.MIP.Solution.Printemps
       Numeric.Optimization.MIP.Solution.SCIP
       Numeric.Optimization.MIP.Solver
       Numeric.Optimization.MIP.Solver.Base
@@ -81,7 +135,9 @@
       Numeric.Optimization.MIP.Solver.CPLEX
       Numeric.Optimization.MIP.Solver.Glpsol
       Numeric.Optimization.MIP.Solver.GurobiCl
+      Numeric.Optimization.MIP.Solver.HiGHS
       Numeric.Optimization.MIP.Solver.LPSolve
+      Numeric.Optimization.MIP.Solver.Printemps
       Numeric.Optimization.MIP.Solver.SCIP
   other-modules:
       Paths_MIP
@@ -93,14 +149,16 @@
       CPP
   build-depends:
       OptDir
-    , base >=4.7 && <5
+    , aeson >=1.4.2.0
+    , base >=4.12 && <5
     , containers >=0.5.0
     , data-default-class
     , extended-reals >=0.1 && <1.0
     , filepath
+    , hashable >=1.2.7.0 && <1.6.0.0
     , intern >=0.9.1.2 && <1.0.0.0
     , lattices
-    , megaparsec >=5 && <10
+    , megaparsec >=7 && <10
     , mtl >=2.1.2 && <2.4
     , process >=1.1.0.2
     , scientific
@@ -123,6 +181,7 @@
   type: exitcode-stdio-1.0
   main-is: TestSuite.hs
   other-modules:
+      IsClose
       Test.LPFile
       Test.MIP
       Test.MIPSolver
@@ -131,15 +190,19 @@
   hs-source-dirs:
       test
   build-depends:
-      MIP
-    , base >=4.7 && <5
+      HUnit
+    , MIP
+    , base >=4.12 && <5
     , containers >=0.5.0
     , data-default-class
+    , filepath
     , lattices
+    , quickcheck-instances >=0.3.19 && <0.4
     , tasty >=0.10.1
     , tasty-hunit >=0.9 && <0.11
     , tasty-quickcheck >=0.8 && <0.12
     , tasty-th
+    , temporary
   default-language: Haskell2010
   if flag(TestCBC)
     cpp-options: -DTEST_CBC
@@ -157,11 +220,23 @@
     cpp-options: -DTEST_GUROBI_CL
   else
     cpp-options:  
+  if flag(TestHiGHS)
+    cpp-options: -DTEST_HIGHS
+  else
+    cpp-options:  
   if flag(TestLPSolve)
     cpp-options: -DTEST_LP_SOLVE
   else
     cpp-options:  
+  if flag(TestPrintemps)
+    cpp-options: -DTEST_PRINTEMPS
+  else
+    cpp-options:  
   if flag(TestSCIP)
     cpp-options: -DTEST_SCIP
+  else
+    cpp-options:  
+  if flag(WithZlib)
+    cpp-options: -DWITH_ZLIB
   else
     cpp-options:  
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 
 [![Hackage](https://img.shields.io/hackage/v/MIP.svg)](https://hackage.haskell.org/package/MIP)
 [![Hackage Deps](https://img.shields.io/hackage-deps/v/MIP.svg)](https://packdeps.haskellers.com/feed?needle=MIP)
-[![License](https://img.shields.io/badge/License-BSD%203--Clause-blue.svg)](https://opensource.org/licenses/BSD-3-Clause)
+[![License: BSD 3-Clause](https://img.shields.io/badge/License-BSD%203--Clause-blue.svg)](https://opensource.org/licenses/BSD-3-Clause)
 
 Library for using Mixed Integer Programming (MIP) in Haskell.
 This library contains functions like:
@@ -10,3 +10,79 @@
 * 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.
+
+## Examples
+
+### Convert LP file into MPS file
+
+```haskell
+import qualified Numeric.Optimization.MIP as MIP
+
+main :: IO ()
+main = do
+  prob <- MIP.readFile MIP.def "samples/lp/test.lp"
+  MIP.writeFile MIP.def "test.mps" prob
+```
+
+### Solve LP file using the CbC solver
+
+
+```haskell
+import Control.Monad
+import qualified Data.Map.Lazy as Map
+import qualified Data.Text as T
+
+import qualified Numeric.Optimization.MIP as MIP
+import Numeric.Optimization.MIP.Solver
+
+main :: IO ()
+main = do
+  prob <- MIP.readFile MIP.def "samples/lp/test.lp"
+  sol <- solve cbc MIP.def{ solveTimeLimit = Just 10.0 } prob
+  print $ MIP.solStatus sol
+  putStrLn $ "Objective Value: " ++ show (MIP.solObjectiveValue sol)
+  forM_ (MIP.variables prob) $ \v -> do
+    putStrLn $ T.unpack (MIP.varName v) ++ " = " ++ show (MIP.solVariables sol Map.! v)
+```
+
+### Construcing a problem instance and solving it using the CbC solver
+
+```haskell
+{-# LANGUAGE OverloadedStrings #-}
+import Control.Monad
+import qualified Data.Map.Lazy as Map
+import qualified Data.Text as T
+
+import qualified Numeric.Optimization.MIP as MIP
+import Numeric.Optimization.MIP ((.<=.))
+import Numeric.Optimization.MIP.Solver
+
+-- Example from https://en.wikipedia.org/wiki/Integer_programming
+main :: IO ()
+main = do
+  let [x, y] = map MIP.varExpr ["x", "y"]
+      prob =
+        MIP.def
+        { MIP.objectiveFunction =
+            MIP.def
+            { MIP.objDir = MIP.OptMax
+            , MIP.objExpr = y
+            }
+        , MIP.constraints =
+            [ - x +   y .<=. 1
+            , 3*x + 2*y .<=. 12
+            , 2*x + 3*y .<=. 12
+            ]
+        , MIP.varDomains =
+            Map.fromList
+            [ ("x", (MIP.IntegerVariable, (0, MIP.PosInf)))
+            , ("y", (MIP.IntegerVariable, (0, MIP.PosInf)))
+            ]
+        }
+
+  sol <- solve cbc MIP.def{ solveTimeLimit = Just 10.0 } prob
+  print $ MIP.solStatus sol
+  putStrLn $ "Objective Value: " ++ show (MIP.solObjectiveValue sol)
+  forM_ (MIP.variables prob) $ \v -> do
+    putStrLn $ T.unpack (MIP.varName v) ++ " = " ++ show (MIP.solVariables sol Map.! v)
+```
diff --git a/samples/lp/empty_obj_1.lp b/samples/lp/empty_obj_1.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/empty_obj_1.lp
@@ -0,0 +1,15 @@
+\GLPK: parse error
+\Gurobi: OK
+\CPLEX: OK
+\SCIP: OK
+Maximize
+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
diff --git a/samples/lp/empty_obj_2.lp b/samples/lp/empty_obj_2.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/empty_obj_2.lp
@@ -0,0 +1,16 @@
+\GLPK: parse error
+\Gurobi: OK
+\CPLEX: OK
+\SCIP: OK
+Maximize
+obj: 
+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
diff --git a/samples/lp/glpk-preprocess-bug.lp b/samples/lp/glpk-preprocess-bug.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/glpk-preprocess-bug.lp
@@ -0,0 +1,20 @@
+\ http://lists.gnu.org/archive/html/help-glpk/2010-12/msg00055.html
+\ From:    Andrew Makhorin
+\ Subject: Re: [Help-glpk] Skip preprocessing
+\ Date:	   Thu, 16 Dec 2010 15:16:23 +0300
+
+\ Your instance has a property that lp relaxation to it has no finite
+\ maximum while the integer feasible region is empty (in fact, 1/3 <= x1
+\ - x2 <= 3/4, so no integers x1 and x2 can satisfy this constraint).
+
+Maximize
+  obj: x2
+Subject To
+  c1: 3 x1 - 3 x2 >= 1
+  c2: 4 x1 - 4 x2 <= 3
+Bounds 
+  1 <= x1
+  0 <= x2
+General
+  x1 x2
+End
diff --git a/samples/lp/infeasible.lp b/samples/lp/infeasible.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/infeasible.lp
@@ -0,0 +1,10 @@
+Minimize
+  0 x
+Subject To
+  3 x >= 1
+  3 x <= 2
+Bounds
+  -inf <= x <= +inf
+General
+  x
+End
diff --git a/samples/lp/test-bounds-fixed.lp b/samples/lp/test-bounds-fixed.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-bounds-fixed.lp
@@ -0,0 +1,7 @@
+Maximize
+ x1
+Subject To
+ x1 - x2 >= 0
+Bounds
+ x1 = 2
+End
diff --git a/samples/lp/test-bounds-free.lp b/samples/lp/test-bounds-free.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-bounds-free.lp
@@ -0,0 +1,7 @@
+Maximize
+ x1
+Subject To
+ x1 - x2 >= 0
+Bounds
+ x1 free
+End
diff --git a/samples/lp/test-indicator.lp b/samples/lp/test-indicator.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-indicator.lp
@@ -0,0 +1,14 @@
+\ GLPK: error
+\ Gurobi 5.0.1: error
+\ SCIP: OK
+Maximize
+ obj: x1 + 2 x2 + 3 x3 + x4
+Subject To
+ c1: - x1 + x2 + x3 + 10 x4 <= 20
+ c2: x4 = 1 -> x1 - 3 x2 + x3 <= 30
+ c3: x2 - 3.5 x4 = 0
+Bounds
+ 0 <= x1 <= 40
+Binaries
+ x4
+End
diff --git a/samples/lp/test-lazy-constraints.lp b/samples/lp/test-lazy-constraints.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-lazy-constraints.lp
@@ -0,0 +1,13 @@
+Maximize
+ obj: x1 + 2 x2 + 3 x3 + x4
+Subject To
+ c1: - x1 + x2 + x3 + 10 x4 <= 20
+ c2: x1 - 3 x2 + x3 <= 30
+Lazy Constraints
+ c3: x2 - 3.5 x4 = 0
+Bounds
+ 0 <= x1 <= 40
+ 2 <= x4 <= 3
+General
+ x4
+End
diff --git a/samples/lp/test-problem-name.lp b/samples/lp/test-problem-name.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-problem-name.lp
@@ -0,0 +1,13 @@
+\* Problem: foo bar *\
+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
diff --git a/samples/lp/test-qcp.lp b/samples/lp/test-qcp.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-qcp.lp
@@ -0,0 +1,16 @@
+\ http://www.gurobi.com/documentation/5.0/reference-manual/node746
+\ Gurobi 5.0.1: Ok
+\ SCIP 2.1.1: error
+
+Maximize
+  x + y + z
+Subject To
+  c0: x + y = 1
+  c1: x + 5 y + 2 z <= 10
+  qc0: x + y + [ x ^ 2 - 2 x * y + 3 y ^ 2 ] <= 5
+Bounds
+  0 <= x <= 5
+  z >= 2
+Generals
+  x y z
+End
diff --git a/samples/lp/test-qcp2.lp b/samples/lp/test-qcp2.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-qcp2.lp
@@ -0,0 +1,15 @@
+\ Gurobi 5.0.1: Ok
+\ SCIP 2.1.1: error
+
+Maximize
+  x + y + z
+Subject To
+  c0: x + y = 1
+  c1: x + 5 y + 2 z <= 10
+  qc0: x + y + [ x ^ 2 - 2 x * y + 3 y ^ 2 ] / 2 <= 5
+Bounds
+  0 <= x <= 5
+  z >= 2
+Generals
+  x y z
+End
diff --git a/samples/lp/test-qp.lp b/samples/lp/test-qp.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-qp.lp
@@ -0,0 +1,8 @@
+\ GLPK: error
+\ Gurobi 5.0.1: Ok
+\ SCIP: okay? but "Quadratic constraint handler does not have LAPACK for eigenvalue computation. Will assume that matrices (with size > 2x2) are indefinite."
+Minimize
+obj: a + b + [ a^2 + 4 a * b + 7 b^2 ] / 2
+Subject To
+c1: a + b >= 10
+End
diff --git a/samples/lp/test-semicont.lp b/samples/lp/test-semicont.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-semicont.lp
@@ -0,0 +1,13 @@
+Minimize
+ obj: -2 x3
+Subject To
+ c1: +x2 -x1 <= 10
+ c2: +x3 +x2 +x1 <= 20
+Bounds
+ x3 >= 2
+ x3 <= 3
+ x1 >= 2.1
+ x1 <= 30
+Semi-continuous
+ x3 x1
+End
diff --git a/samples/lp/test-semiint.lp b/samples/lp/test-semiint.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-semiint.lp
@@ -0,0 +1,16 @@
+\ Converted from Example 2 of http://lpsolve.sourceforge.net/5.5/Xpress-format.htm
+Minimize
+ obj: -2 x3
+Subject To
+ c1: +x2 -x1 <= 10
+ c2: +x3 +x2 +x1 <= 20
+Bounds
+ x3 >= 2
+ x3 <= 3
+ x1 >= 2.1
+ x1 <= 30
+General
+ x3 x1
+Semi-continuous
+ x3 x1
+End
diff --git a/samples/lp/test-sos.lp b/samples/lp/test-sos.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-sos.lp
@@ -0,0 +1,16 @@
+MINIMIZE
+obj: - x4 - 3.0 x3 - 2.0 x2 - x1
+SUBJECT TO
+c1: + 10.0 x4 + x3 + x2 - x1 <= 20.0
+c2: + x3 - 3.0 x2 + x1 <= 30.0
+c3: - 3.5 x4 + x2 = 0.0
+BOUNDS
+2.0 <= x4 <= 3.0
+0.0 <= x3 <= +inf
+0.0 <= x2 <= +inf
+0.0 <= x1 <= 40.0
+GENERALS
+x4
+SOS
+set1: S1 ::  x1 : 10000.0  x2 : 20000.0  x4 : 40000.0
+END
diff --git a/samples/lp/test-sos2.lp b/samples/lp/test-sos2.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-sos2.lp
@@ -0,0 +1,16 @@
+MINIMIZE
+  obj: x
+SUBJECT TO
+  c1: - x + 0 lam1 + 0.5 lam2 + 1 lam3 + 1.5 lam4 + 2.5 lam5 = 0
+  c2: - y + 0 lam1 + 0.125 lam2 + 1 lam3 + 3.375 lam4 + 15.625 lam5 = 0
+  c3: lam1 + lam2 + lam3 + lam4 + lam5 = 1
+BOUNDS
+  0 <= lam1 <= 1
+  0 <= lam2 <= 1
+  0 <= lam3 <= 1
+  0 <= lam4 <= 1
+  0 <= lam5 <= 1
+  5 <= y <= +inf
+SOS
+  set1: S2 :: lam1 : 1 lam2 : 2 lam3 : 3 lam4 : 4 lam5 : 5
+END
diff --git a/samples/lp/test-user-cuts.lp b/samples/lp/test-user-cuts.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test-user-cuts.lp
@@ -0,0 +1,14 @@
+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
+User Cuts
+ c4: x1 + x3 - 10.5 x4 <= 30
+Bounds
+ 0 <= x1 <= 40
+ 2 <= x4 <= 3
+General
+ x4
+End
diff --git a/samples/lp/test.lp b/samples/lp/test.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/test.lp
@@ -0,0 +1,13 @@
+\...
+Maximize
+ obj: x1 + 2 x2 + 3 x3 + x4
+Subject To
+ c1: - x1 + x2 + x3 + 10 x4 <= 20
+ c2: x1 - 3 x2 + x3 <= 30
+ c3: x2 - 3.5 x4 = 0
+Bounds
+ 0 <= x1 <= 40
+ 2 <= x4 <= 3
+General
+ x4
+End
diff --git a/samples/lp/unbounded-ip.lp b/samples/lp/unbounded-ip.lp
new file mode 100644
--- /dev/null
+++ b/samples/lp/unbounded-ip.lp
@@ -0,0 +1,10 @@
+Maximize
+  obj: x + y
+Subject To
+  c1: y - 2 x >= 0
+Bounds 
+  1 <= x
+  0 <= y
+General
+  x y
+End
diff --git a/samples/mps/example2.mps b/samples/mps/example2.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/example2.mps
@@ -0,0 +1,18 @@
+* http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_synopsis.html
+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
diff --git a/samples/mps/ind1.mps b/samples/mps/ind1.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/ind1.mps
@@ -0,0 +1,25 @@
+* http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_indicators.html
+NAME          ind1.mps
+ROWS
+ N  obj
+ L  row2
+ L  row4
+ E  row1
+ E  row3
+COLUMNS
+    x         obj                            -1
+    x         row2                            1
+    x         row4                            1
+    x         row1                            1
+    y         row4                            1
+    z         row4                            1
+    z         row3                            1
+RHS
+    rhs       row2                           10
+    rhs       row4                           15
+BOUNDS
+ UI bnd       y                               1
+INDICATORS
+ IF row1      y                               1
+ IF row3      y                               0
+ENDATA
diff --git a/samples/mps/intvar1.mps b/samples/mps/intvar1.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/intvar1.mps
@@ -0,0 +1,25 @@
+* http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_intvar.html
+NAME
+ROWS
+ N  obj     
+ L  c1      
+ L  c2      
+ E  c3      
+COLUMNS
+    x1        obj                 -1   c1                  -1
+    x1        c2                   1
+    x2        obj                 -2   c1                   1
+    x2        c2                  -3   c3                   1
+    x3        obj                 -3   c1                   1
+    x3        c2                   1
+    MARK0000  'MARKER'                 'INTORG'
+    x4        obj                 -1   c1                  10
+    x4        c3                -3.5
+    MARK0001  'MARKER'                 'INTEND'
+RHS
+    rhs       c1                  20   c2                  30
+BOUNDS
+ UP BOUND     x1                  40
+ LO BOUND     x4                   2
+ UP BOUND     x4                   3
+ENDATA
diff --git a/samples/mps/intvar2.mps b/samples/mps/intvar2.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/intvar2.mps
@@ -0,0 +1,23 @@
+* http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_intvar.html
+NAME
+ROWS
+ N  obj     
+ L  c1      
+ L  c2      
+ E  c3      
+COLUMNS
+    x1        obj                 -1   c1                  -1
+    x1        c2                   1
+    x2        obj                 -2   c1                   1
+    x2        c2                  -3   c3                   1
+    x3        obj                 -3   c1                   1
+    x3        c2                   1
+    x4        obj                 -1   c1                  10
+    x4        c3                 -3.5
+RHS
+    rhs       c1                  20   c2                  30
+BOUNDS
+ UP BOUND     x1                  40
+ LI BOUND     x4                   2
+ UI BOUND     x4                   3
+ENDATA
diff --git a/samples/mps/quadobj1.mps b/samples/mps/quadobj1.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/quadobj1.mps
@@ -0,0 +1,16 @@
+* http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_quadobj.html
+NAME          problem
+ROWS
+ N  obj
+ G  c1
+COLUMNS
+    a         obj                  1   c1                   1
+    b         obj                  1   c1                   1
+RHS
+    rhs       c1                  10
+QMATRIX
+    a         a                    1
+    a         b                    2
+    b         a                    2
+    b         b                    7
+ENDATA
diff --git a/samples/mps/quadobj2.mps b/samples/mps/quadobj2.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/quadobj2.mps
@@ -0,0 +1,15 @@
+* http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_quadobj.html
+NAME          problem
+ROWS
+ N  obj
+ G  c1
+COLUMNS
+    a         obj                  1   c1                   1
+    b         obj                  1   c1                   1
+RHS
+    rhs       c1                  10
+QUADOBJ
+    a         a                    1
+    a         b                    2
+    b         b                    7
+ENDATA
diff --git a/samples/mps/ranges.mps b/samples/mps/ranges.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/ranges.mps
@@ -0,0 +1,29 @@
+NAME          range.mps
+ROWS
+ N  obj     
+ E  c1      
+ E  c2      
+ G  c3      
+ L  c4      
+COLUMNS
+    x1        obj                 -1   c1                   1
+    x1        c2                  1    c3                   1
+    x1        c4                  1
+    x2        obj                 -2   c1                   -3
+    x2        c2                  -3   c3                   -3
+    x2        c4                  -3
+    x3        obj                 -3   c1                   -3
+    x3        c2                  -3   c3                  -3
+    x3        c4                  -3
+RHS
+    rhs       c1                  15
+    rhs       c2                  15
+    rhs       c3                  10
+    rhs       c4                  20
+RANGES
+    rhs       c1                  15
+    rhs       c2                  -15
+    rhs       c3                  10
+    rhs       c4                  10
+BOUNDS
+ENDATA
diff --git a/samples/mps/sc.mps b/samples/mps/sc.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/sc.mps
@@ -0,0 +1,11 @@
+NAME          sc.mps
+ROWS
+ N  obj     
+ E  c1      
+COLUMNS
+    x1        obj                 -1   c1                  1
+RHS
+    rhs       c1                  0.5
+BOUNDS
+ SC BOUND     x1                  40
+ENDATA
diff --git a/samples/mps/sos.mps b/samples/mps/sos.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/sos.mps
@@ -0,0 +1,28 @@
+* http://pic.dhe.ibm.com/infocenter/cosinfoc/v12r4/topic/ilog.odms.cplex.help/CPLEX/File_formats_reference/topics/MPS_ext_sos.html
+NAME
+ROWS
+ N  obj     
+ L  c1      
+ L  c2      
+ E  c3      
+COLUMNS
+    x1        obj                 -1   c1                  -1
+    x1        c2                   1
+    x2        obj                 -2   c1                   1
+    x2        c2                  -3   c3                   1
+    x3        obj                 -3   c1                   1
+    x3        c2                   1
+    x4        obj                 -1   c1                  10
+    x4        c3                 -3.5
+RHS
+    rhs       c1                  20   c2                  30
+BOUNDS
+ UP BOUND     x1                  40
+ LI BOUND     x4                   2
+ UI BOUND     x4                   3
+SOS
+ S1 set1
+    x1                  10000
+    x2                  20000
+    x4                  40000
+ENDATA
diff --git a/samples/mps/test-bounds-bv.mps b/samples/mps/test-bounds-bv.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-bounds-bv.mps
@@ -0,0 +1,12 @@
+NAME          prob
+ROWS
+ N  obj1
+ G  row1
+COLUMNS
+    x2        row1      -1.0
+    x1        obj1      1.0
+    x1        row1      1.0
+RHS
+BOUNDS
+ BV bound     x1
+ENDATA
diff --git a/samples/mps/test-bounds-fixed.mps b/samples/mps/test-bounds-fixed.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-bounds-fixed.mps
@@ -0,0 +1,12 @@
+NAME          prob
+ROWS
+ N  obj1
+ G  row1
+COLUMNS
+    x2        row1      -1.0
+    x1        obj1      1.0
+    x1        row1      1.0
+RHS
+BOUNDS
+ FX bound     x1        2.0
+ENDATA
diff --git a/samples/mps/test-bounds-free.mps b/samples/mps/test-bounds-free.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-bounds-free.mps
@@ -0,0 +1,12 @@
+NAME          prob
+ROWS
+ N  obj1
+ G  row1
+COLUMNS
+    x2        row1      -1.0
+    x1        obj1      1.0
+    x1        row1      1.0
+RHS
+BOUNDS
+ FR bound     x1
+ENDATA
diff --git a/samples/mps/test-bounds-mi.mps b/samples/mps/test-bounds-mi.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-bounds-mi.mps
@@ -0,0 +1,12 @@
+NAME          prob
+ROWS
+ N  obj1
+ G  row1
+COLUMNS
+    x2        row1      -1.0
+    x1        obj1      1.0
+    x1        row1      1.0
+RHS
+BOUNDS
+ MI bound     x1
+ENDATA
diff --git a/samples/mps/test-bounds-pl.mps b/samples/mps/test-bounds-pl.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-bounds-pl.mps
@@ -0,0 +1,12 @@
+NAME          prob
+ROWS
+ N  obj1
+ G  row1
+COLUMNS
+    x2        row1      -1.0
+    x1        obj1      1.0
+    x1        row1      1.0
+RHS
+BOUNDS
+ PL bound     x1
+ENDATA
diff --git a/samples/mps/test-implicit-bv.mps b/samples/mps/test-implicit-bv.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-implicit-bv.mps
@@ -0,0 +1,22 @@
+NAME
+ROWS
+ N  obj     
+ L  c1      
+ L  c2      
+ E  c3      
+COLUMNS
+    x1        obj                 -1   c1                  -1
+    x1        c2                   1
+    x2        obj                 -2   c1                   1
+    x2        c2                  -3   c3                   1
+    x3        obj                 -3   c1                   1
+    x3        c2                   1
+    MARK0000  'MARKER'                 'INTORG'
+    x4        obj                 -1   c1                  10
+    x4        c3                -3.5
+    MARK0001  'MARKER'                 'INTEND'
+RHS
+    rhs       c1                  20   c2                  30
+BOUNDS
+ UP BOUND     x1                  40
+ENDATA
diff --git a/samples/mps/test-lazy-constraints.mps b/samples/mps/test-lazy-constraints.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-lazy-constraints.mps
@@ -0,0 +1,31 @@
+NAME          prob
+ROWS
+ N  obj
+ L  c1
+ L  c2
+LAZYCONS
+ E  c3
+COLUMNS
+    x3        c1        1.0
+    x3        c2        1.0
+    x3        obj       3.0
+    x2        c1        1.0
+    x2        c2        -3.0
+    x2        c3        1.0
+    x2        obj       2.0
+    x1        c1        -1.0
+    x1        c2        1.0
+    x1        obj       1.0
+    MARK0000  'MARKER'                 'INTORG'
+    x4        c1        10.0
+    x4        c3        -3.5
+    x4        obj       1.0
+    MARK0001  'MARKER'                 'INTEND'
+RHS
+    rhs       c1        20.0
+    rhs       c2        30.0
+BOUNDS
+ LI bound     x4        2.0
+ UI bound     x4        3.0
+ UP bound     x1        40.0
+ENDATA
diff --git a/samples/mps/test-negative-upper-bound.mps b/samples/mps/test-negative-upper-bound.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-negative-upper-bound.mps
@@ -0,0 +1,12 @@
+NAME          prob
+ROWS
+ N  obj1
+ G  row1
+COLUMNS
+    x2        row1      -1.0
+    x1        obj1      1.0
+    x1        row1      1.0
+RHS
+BOUNDS
+ UP bound     x1        -2.0
+ENDATA
diff --git a/samples/mps/test-qcp.mps b/samples/mps/test-qcp.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-qcp.mps
@@ -0,0 +1,38 @@
+NAME          
+OBJSENSE
+ MAX
+ROWS
+ N  obj1
+ E  c0
+ L  c1
+ L  qc0
+COLUMNS
+    MARK0000  'MARKER'                 'INTORG'
+    x         c0        1.0
+    x         c1        1.0
+    x         obj1      1.0
+    x         qc0       1.0
+    y         c0        1.0
+    y         c1        5.0
+    y         obj1      1.0
+    y         qc0       1.0
+    z         c1        2.0
+    z         obj1      1.0
+    MARK0001  'MARKER'                 'INTEND'
+RHS
+    rhs       c0        1.0
+    rhs       c1        10.0
+    rhs       qc0       5.0
+BOUNDS
+ LI bound     x         0.0
+ UI bound     x         5.0
+ LI bound     y         0.0
+ PL bound     y
+ LI bound     z         2.0
+ PL bound     z
+QCMATRIX   qc0
+    x         x         1.0
+    x         y         -1.0
+    y         x         -1.0
+    y         y         3.0
+ENDATA
diff --git a/samples/mps/test-semicont.mps b/samples/mps/test-semicont.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-semicont.mps
@@ -0,0 +1,21 @@
+NAME          prob
+ROWS
+ N  obj
+ L  c1
+ L  c2
+COLUMNS
+    x3        c2        1.0
+    x3        obj       -2.0
+    x2        c1        1.0
+    x2        c2        1.0
+    x1        c1        -1.0
+    x1        c2        1.0
+RHS
+    rhs       c1        10.0
+    rhs       c2        20.0
+BOUNDS
+ LO bound     x3        2.0
+ SC bound     x3        3.0
+ LO bound     x1        2.1
+ SC bound     x1        30.0
+ENDATA
diff --git a/samples/mps/test-semiint-sc.mps b/samples/mps/test-semiint-sc.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-semiint-sc.mps
@@ -0,0 +1,23 @@
+NAME          prob
+ROWS
+ N  obj
+ L  c1
+ L  c2
+COLUMNS
+    x2        c1        1.0
+    x2        c2        1.0
+    MARK0000  'MARKER'                 'INTORG'
+    x3        c2        1.0
+    x3        obj       -2.0
+    x1        c1        -1.0
+    x1        c2        1.0
+    MARK0001  'MARKER'                 'INTEND'
+RHS
+    rhs       c1        10.0
+    rhs       c2        20.0
+BOUNDS
+ LO bound     x3        2.0
+ SC bound     x3        3.0
+ LO bound     x1        2.1
+ SC bound     x1        30.0
+ENDATA
diff --git a/samples/mps/test-semiint-si.mps b/samples/mps/test-semiint-si.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-semiint-si.mps
@@ -0,0 +1,21 @@
+NAME          prob
+ROWS
+ N  obj
+ L  c1
+ L  c2
+COLUMNS
+    x2        c1        1.0
+    x2        c2        1.0
+    x3        c2        1.0
+    x3        obj       -2.0
+    x1        c1        -1.0
+    x1        c2        1.0
+RHS
+    rhs       c1        10.0
+    rhs       c2        20.0
+BOUNDS
+ LO bound     x3        2.0
+ SI bound     x3        3.0
+ LO bound     x1        2.1
+ SI bound     x1        30.0
+ENDATA
diff --git a/samples/mps/test-sos2.mps b/samples/mps/test-sos2.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-sos2.mps
@@ -0,0 +1,40 @@
+NAME          prob
+ROWS
+ N  obj
+ E  c1
+ E  c2
+ E  c3
+COLUMNS
+    lam3      c1        1.0
+    lam3      c2        1.0
+    lam3      c3        1.0
+    lam2      c1        0.5
+    lam2      c2        0.125
+    lam2      c3        1.0
+    lam1      c3        1.0
+    lam5      c1        2.5
+    lam5      c2        15.625
+    lam5      c3        1.0
+    lam4      c1        1.5
+    lam4      c2        3.375
+    lam4      c3        1.0
+    x         c1        -1.0
+    x         obj       1.0
+    y         c2        -1.0
+RHS
+    rhs       c3        1.0
+BOUNDS
+ UP bound     lam3      1.0
+ UP bound     lam2      1.0
+ UP bound     lam1      1.0
+ UP bound     lam5      1.0
+ UP bound     lam4      1.0
+ LO bound     y         5.0
+SOS
+ S2 set1
+    lam1      1.0
+    lam2      2.0
+    lam3      3.0
+    lam4      4.0
+    lam5      5.0
+ENDATA
diff --git a/samples/mps/test-user-cuts.mps b/samples/mps/test-user-cuts.mps
new file mode 100644
--- /dev/null
+++ b/samples/mps/test-user-cuts.mps
@@ -0,0 +1,36 @@
+NAME          prob
+ROWS
+ N  obj
+ L  c1
+ L  c2
+ E  c3
+USERCUTS
+ L  c4
+COLUMNS
+    x3        c1        1.0
+    x3        c2        1.0
+    x3        c4        1.0
+    x3        obj       3.0
+    x2        c1        1.0
+    x2        c2        -3.0
+    x2        c3        1.0
+    x2        obj       2.0
+    x1        c1        -1.0
+    x1        c2        1.0
+    x1        c4        1.0
+    x1        obj       1.0
+    MARK0000  'MARKER'                 'INTORG'
+    x4        c1        10.0
+    x4        c3        -3.5
+    x4        c4        -10.5
+    x4        obj       1.0
+    MARK0001  'MARKER'                 'INTEND'
+RHS
+    rhs       c1        20.0
+    rhs       c2        30.0
+    rhs       c4        30.0
+BOUNDS
+ LI bound     x4        2.0
+ UI bound     x4        3.0
+ UP bound     x1        40.0
+ENDATA
diff --git a/src/Numeric/Optimization/MIP.hs b/src/Numeric/Optimization/MIP.hs
--- a/src/Numeric/Optimization/MIP.hs
+++ b/src/Numeric/Optimization/MIP.hs
@@ -4,6 +4,7 @@
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE OverloadedStrings #-}
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Optimization.MIP
@@ -19,11 +20,84 @@
 -----------------------------------------------------------------------------
 module Numeric.Optimization.MIP
   (
-  -- * The MIP Problem/Solution types
-    module Numeric.Optimization.MIP.Base
+  -- * Mixed-Integer Programming (MIP) problem specification
 
+  -- ** MIP problems
+    Problem (..)
+
+  -- *** Set of variables
+  , variables
+  , continuousVariables
+  , integerVariables
+  , binaryVariables
+  , semiContinuousVariables
+  , semiIntegerVariables
+
+  -- *** Variable's attributes
+  , varTypes
+  , varType
+  , getVarType
+  , varBounds
+  , getBounds
+
+  -- ** Variables
+  , Var (Var)
+  , varName
+  , toVar
+  , fromVar
+
+  -- *** Variable types
+  , VarType (..)
+
+  -- *** Variable bounds
+  , BoundExpr
+  , Extended (..)
+  , Bounds
+  , defaultBounds
+  , defaultLB
+  , defaultUB
+
+  -- ** Labels
+  , Label
+
+  -- ** Expressions
+  , Expr (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
+
+  -- * Evaluation
+  , Tol (..)
+  , zeroTol
+  , Eval (..)
+
   -- * File I/O
   -- $IO
+  , FileOptions (..)
+  , WriteSetting (..)
 
   -- ** Reading problem files
   , readFile
@@ -39,6 +113,11 @@
   , writeMPSFile
   , toLPString
   , toMPSString
+
+  -- * Utilities
+  , Default (..)
+  , Variables (..)
+  , intersectBounds
   ) where
 
 import Prelude hiding (readFile, writeFile)
@@ -118,19 +197,11 @@
 #endif
 
 -- | Parse a string containing LP file data.
-#if MIN_VERSION_megaparsec(6,0,0)
 parseLPString :: (Stream s, Token s ~ Char, IsString (Tokens s)) => FileOptions -> String -> s -> Either (ParseError s) (Problem Scientific)
-#else
-parseLPString :: (Stream s, Token s ~ Char) => FileOptions -> String -> s -> Either (ParseError s) (Problem Scientific)
-#endif
 parseLPString = LPFile.parseString
 
 -- | Parse a string containing MPS file data.
-#if MIN_VERSION_megaparsec(6,0,0)
 parseMPSString :: (Stream s, Token s ~ Char, IsString (Tokens s)) => FileOptions -> String -> s -> Either (ParseError s) (Problem Scientific)
-#else
-parseMPSString :: (Stream s, Token s ~ Char) => FileOptions -> String -> s -> Either (ParseError s) (Problem Scientific)
-#endif
 parseMPSString = MPSFile.parseString
 
 -- | Generate LP file or MPS file based on file extension.
diff --git a/src/Numeric/Optimization/MIP/Base.hs b/src/Numeric/Optimization/MIP/Base.hs
--- a/src/Numeric/Optimization/MIP/Base.hs
+++ b/src/Numeric/Optimization/MIP/Base.hs
@@ -1,6 +1,11 @@
 {-# OPTIONS_GHC -Wall #-}
 {-# OPTIONS_HADDOCK show-extensions #-}
 {-# LANGUAGE CPP #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE PatternSynonyms #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE ViewPatterns #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Optimization.MIP.Base
@@ -16,55 +21,67 @@
 -----------------------------------------------------------------------------
 module Numeric.Optimization.MIP.Base
   (
-  -- * The MIP Problem type
+  -- * Mixed-Integer Programming (MIP) problem specification
+
+  -- ** MIP problems
     Problem (..)
-  , Label
 
-  -- * Variables
-  , Var
+  -- *** Set of variables
+  , variables
+  , continuousVariables
+  , integerVariables
+  , binaryVariables
+  , semiContinuousVariables
+  , semiIntegerVariables
+
+  -- *** Variable's attributes
+  , varTypes
+  , varType
+  , getVarType
+  , varBounds
+  , getBounds
+
+  -- ** Variables
+  , Var (Var)
+  , varName
   , toVar
   , fromVar
 
-  -- ** Variable types
+  -- *** Variable types
   , VarType (..)
-  , getVarType
 
-  -- ** Variable bounds
+  -- *** Variable bounds
   , BoundExpr
   , Extended (..)
   , Bounds
   , defaultBounds
   , defaultLB
   , defaultUB
-  , getBounds
 
-  -- ** Variable getters
-  , variables
-  , integerVariables
-  , semiContinuousVariables
-  , semiIntegerVariables
+  -- ** Labels
+  , Label
 
-  -- * Expressions
-  , Expr (..)
+  -- ** Expressions
+  , Expr (Expr)
   , varExpr
   , constExpr
   , terms
   , Term (..)
 
-  -- * Objective function
+  -- ** Objective function
   , OptDir (..)
   , ObjectiveFunction (..)
 
-  -- * Constraints
+  -- ** Constraints
 
-  -- ** Linear (or Quadratic or Polynomial) constraints
+  -- *** Linear (or Quadratic or Polynomial) constraints
   , Constraint (..)
   , (.==.)
   , (.<=.)
   , (.>=.)
   , RelOp (..)
 
-  -- ** SOS constraints
+  -- *** SOS constraints
   , SOSType (..)
   , SOSConstraint (..)
 
@@ -73,8 +90,14 @@
   , Status (..)
   , meetStatus
 
-  -- * File I/O options
+  -- * Evaluation
+  , Tol (..)
+  , zeroTol
+  , Eval (..)
+
+  -- * File I/O
   , FileOptions (..)
+  , WriteSetting (..)
 
   -- * Utilities
   , Default (..)
@@ -87,12 +110,19 @@
 #endif
 import Algebra.PartialOrd
 import Control.Arrow ((***))
+import Control.Monad
 import Data.Default.Class
+import Data.Foldable (toList)
+import Data.Hashable
+import Data.List (sortBy)
 import Data.Map (Map)
 import qualified Data.Map as Map
+import Data.Ord (comparing)
+import Data.Sequence (Seq)
+import qualified Data.Sequence as Seq
 import Data.Set (Set)
 import qualified Data.Set as Set
-import Data.Interned (unintern)
+import Data.Interned (intern, unintern)
 import Data.Interned.Text
 import Data.ExtendedReal
 import Data.OptDir
@@ -104,16 +134,23 @@
 
 -- ---------------------------------------------------------------------------
 
--- | Problem
+-- | A problem instance
 data Problem c
   = Problem
   { name :: Maybe T.Text
+    -- ^ Problem name
   , objectiveFunction :: ObjectiveFunction c
+    -- ^ Objective functions of the problem
   , constraints :: [Constraint c]
+    -- ^ Constraints of the problem
+    --
+    -- Indicator constraints and lazy constraints are included in this list.
   , sosConstraints :: [SOSConstraint c]
+    -- ^ Special ordered sets
   , userCuts :: [Constraint c]
-  , varType :: Map Var VarType
-  , varBounds :: Map Var (Bounds c)
+    -- ^ User cuts
+  , varDomains :: Map Var (VarType, Bounds c)
+    -- ^ Variable types and their bounds
   }
   deriving (Show, Eq, Ord)
 
@@ -124,8 +161,7 @@
         , constraints = []
         , sosConstraints = []
         , userCuts = []
-        , varType = Map.empty
-        , varBounds = Map.empty
+        , varDomains = Map.empty
         }
 
 instance Functor Problem where
@@ -135,30 +171,91 @@
     , 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)
+    , varDomains        = fmap (id *** (fmap f *** fmap f)) (varDomains prob)
     }
 
--- | label
+-- | Types of variables
+--
+-- This is equivalent to:
+--
+-- @
+-- 'fmap' 'fst' . 'varDomains'
+-- @
+varTypes :: Problem c -> Map Var VarType
+varTypes = fmap fst . varDomains
+
+{-# DEPRECATED varType "Use varTypes instead" #-}
+-- | Types of variables
+--
+-- Deprecated alias of 'varTypes'.
+varType :: Problem c -> Map Var VarType
+varType = varTypes
+
+-- | Bounds of variables
+--
+-- This is equivalent to:
+--
+-- @
+-- 'fmap' 'snd' . 'varDomains'
+-- @
+varBounds :: Problem c -> Map Var (Bounds c)
+varBounds = fmap snd . varDomains
+
+-- | Label used for naming various elements of t'Problem'
 type Label = T.Text
 
 -- ---------------------------------------------------------------------------
 
--- | variable
-type Var = InternedText
+-- | Variables used in problems
+newtype Var = Var' InternedText
+  deriving Eq
 
+pattern Var :: T.Text -> Var
+pattern Var s <- Var' (unintern -> s) where
+  Var s = Var' (intern s)
+
+{-# COMPLETE Var #-}
+
+instance IsString Var where
+  fromString = Var' . fromString
+
+instance Ord Var where
+  compare (Var' a) (Var' b)
+    | a == b = EQ
+    | otherwise = compare (unintern a) (unintern b)
+
+instance Show Var where
+  showsPrec d (Var x) = showsPrec d x
+
+instance Hashable Var where
+#if MIN_VERSION_intern(0,9,3)
+  hashWithSalt salt (Var' x) = hashWithSalt salt x
+#else
+  hashWithSalt salt (Var' x) = hashWithSalt salt (internedTextId x)
+#endif
+
+-- | Variable's name
+varName :: Var -> T.Text
+varName (Var s) = s
+
+{-# DEPRECATED toVar "Use fromString function or Var pattern instead" #-}
 -- | convert a string into a variable
 toVar :: String -> Var
 toVar = fromString
 
+{-# DEPRECATED fromVar "Use varName function or Var pattern instead" #-}
 -- | convert a variable into a string
 fromVar :: Var -> String
-fromVar = T.unpack . unintern
+fromVar (Var s) = T.unpack s
 
+-- | Type of variables
+--
+-- Variables can take values depending on their types and their bounds ('Bounds').
 data VarType
-  = ContinuousVariable
-  | IntegerVariable
-  | SemiContinuousVariable
-  | SemiIntegerVariable
+  = ContinuousVariable     -- ^ can take values from \(\{x \in \mathbb{R} \mid L \le x \le U\}\)
+  | IntegerVariable        -- ^ can take values from \(\{x \in \mathbb{Z} \mid L \le x \le U\}\)
+  | SemiContinuousVariable -- ^ can take values from \(\{0\} \cup \{x \in \mathbb{R} \mid L \le x \le U\}\)
+  | SemiIntegerVariable    -- ^ can take values from \(\{0\} \cup \{x \in \mathbb{Z} \mid L \le x \le U\}\)
   deriving (Eq, Ord, Show)
 
 instance Default VarType where
@@ -166,7 +263,10 @@
 
 -- | looking up bounds for a variable
 getVarType :: Problem c -> Var -> VarType
-getVarType mip v = Map.findWithDefault def v (varType mip)
+getVarType mip v =
+  case Map.lookup v (varDomains mip) of
+    Just (vt, _) -> vt
+    Nothing -> def
 
 -- | type for representing lower/upper bound of variables
 type BoundExpr c = Extended c
@@ -188,44 +288,68 @@
 
 -- | looking up bounds for a variable
 getBounds :: Num c => Problem c -> Var -> Bounds c
-getBounds mip v = Map.findWithDefault defaultBounds v (varBounds mip)
+getBounds mip v =
+  case Map.lookup v (varDomains mip) of
+    Just (_, bs) -> bs
+    Nothing -> defaultBounds
 
+-- | Intersection of two 'Bounds'
 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)
+-- | Arithmetic expressions
+--
+-- Essentialy an expression is a sequence of t'Term's.
+newtype Expr c = Expr' (Seq (Term c))
+  deriving (Eq, Ord)
 
+pattern Expr :: [Term c] -> Expr c
+pattern Expr ts <- Expr' (toList -> ts) where
+  Expr ts = Expr' (Seq.fromList ts)
+
+{-# COMPLETE Expr #-}
+
+instance Show c => Show (Expr c) where
+  showsPrec d (Expr ts) = showParen (d > app_prec) $
+    showString "Expr " . showsPrec (app_prec+1) ts
+    where
+      app_prec = 10
+
+-- | Variable expression
 varExpr :: Num c => Var -> Expr c
-varExpr v = Expr [Term 1 [v]]
+varExpr v = Expr' $ Seq.singleton $ Term 1 [v]
 
+-- | Constant expression
 constExpr :: (Eq c, Num c) => c -> Expr c
-constExpr 0 = Expr []
-constExpr c = Expr [Term c []]
+constExpr 0 = Expr' Seq.empty
+constExpr c = Expr' $ Seq.singleton $ Term c []
 
+-- | Terms of an expression
 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' (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]
+  negate (Expr' e) = Expr' $ fmap (\(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
+  fmap f (Expr' ts) = Expr' $ fmap (fmap f) ts
 
+-- | Split an expression into an expression without constant term and a constant
 splitConst :: Num c => Expr c -> (Expr c, c)
-splitConst e = (e2, c2)
+splitConst (Expr' ts) = (e2, c2)
   where
-    e2 = Expr [t | t@(Term _ (_:_)) <- terms e]
-    c2 = sum [c | Term c [] <- terms e]
+    p (Term _ (_:_)) = True
+    p _ = False
+    e2 = Expr' $ Seq.filter p ts
+    c2 = sum [c | Term c [] <- toList ts]
 
 -- | terms
 data Term c = Term c [Var]
@@ -258,15 +382,23 @@
 
 -- ---------------------------------------------------------------------------
 
--- | constraint
+-- | Constraint
+--
+-- In the most general case, it is of the form @x = v → L ≤ e ≤ U@.
 data Constraint c
   = Constraint
   { constrLabel     :: Maybe Label
+    -- ^ name of the constraint
   , constrIndicator :: Maybe (Var, c)
+    -- ^ @x = v@ (v is either 0 or 1)
   , constrExpr      :: Expr c
+    -- ^ expression @e@
   , constrLB        :: BoundExpr c
+    -- ^ lower bound @L@
   , constrUB        :: BoundExpr c
+    -- ^ upper bound @U@
   , constrIsLazy    :: Bool
+    -- ^ if it is set to @True@, solver can delay adding the constraint until the constraint is violated.
   }
   deriving (Eq, Ord, Show)
 
@@ -308,7 +440,10 @@
     }
 
 -- | relational operators
-data RelOp = Le | Ge | Eql
+data RelOp
+  = Le  -- ^ (≤)
+  | Ge  -- ^ (≥)
+  | Eql -- ^ (=)
   deriving (Eq, Ord, Enum, Show)
 
 -- ---------------------------------------------------------------------------
@@ -397,8 +532,11 @@
 data Solution r
   = Solution
   { solStatus :: Status
+    -- ^ status
   , solObjectiveValue :: Maybe r
+    -- ^ value of the objective function
   , solVariables :: Map Var r
+    -- ^ variable assignments
   }
   deriving (Eq, Ord, Show)
 
@@ -414,6 +552,120 @@
 
 -- ---------------------------------------------------------------------------
 
+-- | Tolerance for evaluating solutions against t'Problem'.
+data Tol r
+  = Tol
+  { integralityTol :: r
+    -- ^ If a value of integer variable is within this amount from its nearest
+    -- integer, it is considered feasible.
+  , feasibilityTol :: r
+    -- ^ If the amount of violation of constraints is within this amount, it is
+    -- considered feasible.
+  , optimalityTol :: r
+    -- ^ Feasiblity tolerance of dual constraints.
+  }
+
+-- | Defautl is @1e-6@ for the feasibility and optimality tolerances, and @1e-5@ for the integrality tolerance.
+instance Fractional r => Default (Tol r) where
+  def =
+    Tol
+    { integralityTol = 1e-5
+    , feasibilityTol = 1e-6
+    , optimalityTol = 1e-6
+    }
+
+-- | t'Tol' value with all tolerances are zero
+zeroTol :: Fractional r => Tol r
+zeroTol =
+  Tol
+  { integralityTol = 1e-5
+  , feasibilityTol = 1e-6
+  , optimalityTol = 1e-6
+  }
+
+-- | Type class for evaluation various elements of t'Problem' under
+-- the given variable assignments.
+class Eval r a where
+  -- | Result type of 'eval'
+  type Evaluated r a
+
+  -- | Evaluate a value of type @a@ under given assignments and the tolerance
+  eval :: Tol r -> Map Var r -> a -> Evaluated r a
+
+instance Num r => Eval r Var where
+  type Evaluated r Var = r
+  eval _tol sol v =
+    case Map.lookup v sol of
+      Just val -> val
+      Nothing -> 0
+
+instance Num r => Eval r (Term r) where
+  type Evaluated r (Term r) = r
+  eval tol sol (Term c vs) = product (c : [eval tol sol v | v <- vs])
+
+instance Num r => Eval r (Expr r) where
+  type Evaluated r (Expr r) = r
+  eval tol sol expr = sum [eval tol sol t | t <- terms expr]
+
+instance Num r => Eval r (ObjectiveFunction r) where
+  type Evaluated r (ObjectiveFunction r) = r
+  eval tol sol obj = eval tol sol (objExpr obj)
+
+instance (Num r, Ord r) => Eval r (Constraint r) where
+  type Evaluated r (Constraint r) = Bool
+  eval tol sol constr =
+    not (evalIndicator (constrIndicator constr)) ||
+    isInBounds tol (constrLB constr, constrUB constr) (eval tol sol (constrExpr constr))
+    where
+      evalIndicator Nothing = True
+      evalIndicator (Just (v, val')) = isInBounds tol (Finite val', Finite val') (eval tol sol v)
+
+instance (Num r, Ord r) => Eval r (SOSConstraint r) where
+  type Evaluated r (SOSConstraint r) = Bool
+  eval tol sol sos =
+    case sosType sos of
+      S1 -> length [() | val <- body, val] <= 1
+      S2 -> f body
+    where
+      body = map (not . isInBounds tol (0, 0) . eval tol sol . fst) $ sortBy (comparing snd) $ (sosBody sos)
+      f [] = True
+      f [_] = True
+      f (x1 : x2 : xs)
+        | x1 = all not xs
+        | otherwise = f (x2 : xs)
+
+instance (RealFrac r) => Eval r (Problem r) where
+  type Evaluated r (Problem r) = Maybe r
+  eval tol sol prob = do
+    forM_ (Map.toList (varDomains prob)) $ \(v, (vt, bounds)) -> do
+      let val = eval tol sol v
+      case vt of
+        ContinuousVariable -> return ()
+        SemiContinuousVariable -> return ()
+        IntegerVariable -> guard $ isIntegral tol val
+        SemiIntegerVariable -> guard $ isIntegral tol val
+      case vt of
+        ContinuousVariable -> guard $ isInBounds tol bounds val
+        IntegerVariable -> guard $ isInBounds tol bounds val
+        SemiIntegerVariable -> guard $ isInBounds tol (0,0) val || isInBounds tol bounds val
+        SemiContinuousVariable -> guard $ isInBounds tol (0,0) val || isInBounds tol bounds val
+    forM_ (constraints prob) $ \constr -> do
+      guard $ eval tol sol constr
+    forM_ (sosConstraints prob) $ \constr -> do
+      guard $ eval tol sol constr
+    return $ eval tol sol (objectiveFunction prob)
+
+isIntegral :: RealFrac r => Tol r -> r -> Bool
+isIntegral tol x = abs (x - fromIntegral (floor (x + 0.5) :: Integer)) <= integralityTol tol
+
+isInBounds :: (Num r, Ord r) => Tol r -> Bounds r -> r -> Bool
+isInBounds tol (lb, ub) x =
+  lb - Finite (feasibilityTol tol) <= Finite x &&
+  Finite x <= ub + Finite (feasibilityTol tol)
+
+-- ---------------------------------------------------------------------------
+
+-- | Type class for types that contain variables.
 class Variables a where
   vars :: a -> Set Var
 
@@ -433,6 +685,9 @@
 instance Variables (Term c) where
   vars (Term _ xs) = Set.fromList xs
 
+instance Variables Var where
+  vars v = Set.singleton v
+
 instance Variables (ObjectiveFunction c) where
   vars ObjectiveFunction{ objExpr = e } = vars e
 
@@ -446,27 +701,74 @@
 
 -- ---------------------------------------------------------------------------
 
+-- | Set of variables of a t'Problem'
 variables :: Problem c -> Set Var
-variables mip = Map.keysSet $ varType mip
+variables mip = Map.keysSet $ varDomains mip
 
+-- | Set of continuous variables of a t'Problem'
+continuousVariables :: Problem c -> Set Var
+continuousVariables mip = Map.keysSet $ Map.filter ((ContinuousVariable ==) . fst) (varDomains mip)
+
+-- | Set of integer variables of a t'Problem'
 integerVariables :: Problem c -> Set Var
-integerVariables mip = Map.keysSet $ Map.filter (IntegerVariable ==) (varType mip)
+integerVariables mip = Map.keysSet $ Map.filter ((IntegerVariable ==) . fst) (varDomains mip)
 
+-- | Set of binary variables (integers variables with lower bound 0 and upper bound 1) of a t'Problem'
+binaryVariables :: (Num c, Eq c) => Problem c -> Set Var
+binaryVariables mip = Map.keysSet $ Map.filter p (varDomains mip)
+  where
+    p (IntegerVariable, (Finite 0, Finite 1)) = True
+    p (_, _) = False
+
+-- | Set of semi-continuous variables of a t'Problem'
 semiContinuousVariables :: Problem c -> Set Var
-semiContinuousVariables mip = Map.keysSet $ Map.filter (SemiContinuousVariable ==) (varType mip)
+semiContinuousVariables mip = Map.keysSet $ Map.filter ((SemiContinuousVariable ==) . fst) (varDomains mip)
 
+-- | Set of semi-integer variables of a t'Problem'
 semiIntegerVariables :: Problem c -> Set Var
-semiIntegerVariables mip = Map.keysSet $ Map.filter (SemiIntegerVariable ==) (varType mip)
+semiIntegerVariables mip = Map.keysSet $ Map.filter ((SemiIntegerVariable ==) . fst) (varDomains mip)
 
 -- ---------------------------------------------------------------------------
 
+-- | Options for reading/writing problem files
 data FileOptions
   = FileOptions
   { optFileEncoding :: Maybe TextEncoding
+    -- ^ Text encoding used for file input/output
+  , optMPSWriteObjSense :: WriteSetting
+    -- ^ @OBJSENSE@ section in MPS file is an extention of MPS file
+    -- format for specifying the direction of the objective function
+    -- in MPS file. But not all solvers support it (e.g. GLPK-4.48
+    -- does not support it).
+    --
+    -- This option controls whether the @OBJSENSE@ sections is written.
+    -- If 'WriteIfNotDefault' is used, @OBJSENSE@ is written when the
+    -- objective is maximization and @OBJSENSE@ is not written written
+    -- when the objective is minimizing.
+    --
+    -- (Default: 'WriteIfNotDefault')
+  , optMPSWriteObjName :: Bool
+    -- ^ @OBJNAME@ section is an extention of MPS file format for
+    -- selecting an objective function from among the free rows within
+    -- a MPS file. Not all solver support it (e.g. GLPK-4.48
+    -- does not support @OBJNAME@ it).
+    --
+    -- This option controls whether the @OBJNAME@ section is written.
+    --
+    -- (Default: 'True')
   } deriving (Show)
 
 instance Default FileOptions where
   def =
     FileOptions
     { optFileEncoding = Nothing
+    , optMPSWriteObjSense = WriteIfNotDefault
+    , optMPSWriteObjName = True
     }
+
+-- | Options for writing something of not
+data WriteSetting
+  = WriteAlways
+  | WriteIfNotDefault
+  | WriteNever
+  deriving (Eq, Ord, Enum, Bounded, Show, Read)
diff --git a/src/Numeric/Optimization/MIP/FileUtils.hs b/src/Numeric/Optimization/MIP/FileUtils.hs
--- a/src/Numeric/Optimization/MIP/FileUtils.hs
+++ b/src/Numeric/Optimization/MIP/FileUtils.hs
@@ -1,6 +1,5 @@
 {-# OPTIONS_GHC -Wall #-}
 {-# OPTIONS_HADDOCK show-extensions #-}
-{-# LANGUAGE CPP #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Optimization.MIP.FileUtils
@@ -16,18 +15,10 @@
   ( ParseError
   ) where
 
-#if MIN_VERSION_megaparsec(6,0,0)
 import Data.Void
-#endif
 import qualified Text.Megaparsec as MP
 
 -- | Error type for parsing.
 --
 -- The definition is slightly different based on the @megaparsec@ version.
-#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
--- a/src/Numeric/Optimization/MIP/Internal/ProcessUtil.hs
+++ b/src/Numeric/Optimization/MIP/Internal/ProcessUtil.hs
@@ -34,17 +34,22 @@
 runProcessWithOutputCallback
   :: FilePath -- ^ Filename of the executable (see 'proc' for details)
   -> [String] -- ^ any arguments
+  -> Maybe FilePath -- ^ optional path to the working directory for the new process
   -> 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
+runProcessWithOutputCallback cmd args cwd input putMsg putErr = do
   (Just inh, Just outh, Just errh, processh) <- createProcess
     (proc cmd args)
     { std_in  = CreatePipe
     , std_out = CreatePipe
     , std_err = CreatePipe
+    , cwd = cwd
     }
+  hSetBuffering outh LineBuffering
+  hSetBuffering errh LineBuffering
+
   req <- newEmptyTMVarIO
   let f act = atomically (putTMVar req act)
       m1 = forever (hGetLine outh >>= \s -> f (putMsg s))
@@ -59,8 +64,6 @@
       -- 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
diff --git a/src/Numeric/Optimization/MIP/Internal/Util.hs b/src/Numeric/Optimization/MIP/Internal/Util.hs
--- a/src/Numeric/Optimization/MIP/Internal/Util.hs
+++ b/src/Numeric/Optimization/MIP/Internal/Util.hs
@@ -15,81 +15,8 @@
 
 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
--- a/src/Numeric/Optimization/MIP/LPFile.hs
+++ b/src/Numeric/Optimization/MIP/LPFile.hs
@@ -1,12 +1,12 @@
 {-# 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 #-}
+{-# LANGUAGE TypeOperators #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Optimization.MIP.LPFile
@@ -43,13 +43,9 @@
 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.Scientific (Scientific, floatingOrInteger)
 import Data.Map (Map)
 import qualified Data.Map as Map
 import Data.Set (Set)
@@ -60,19 +56,14 @@
 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.Int 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)
@@ -80,19 +71,11 @@
 
 -- ---------------------------------------------------------------------------
 
-#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.
@@ -109,10 +92,8 @@
 
 -- ---------------------------------------------------------------------------
 
-#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
@@ -126,7 +107,7 @@
 comment :: C e s m => m ()
 comment = do
   char '\\'
-  skipManyTill anyChar (try newline)
+  skipManyTill anyChar (try eol)
 
 tok :: C e s m => m a -> m a
 tok p = do
@@ -146,7 +127,7 @@
     syms2 = '.' : syms1
 
 variable :: C e s m => m MIP.Var
-variable = liftM MIP.toVar ident
+variable = liftM fromString ident
 
 label :: C e s m => m MIP.Label
 label = do
@@ -168,16 +149,12 @@
 -- ---------------------------------------------------------------------------
 
 -- | 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"))
+    liftM fromString $ manyTill anyChar (try (string " *\\" >> eol))
   sep
   obj <- problem
 
@@ -214,17 +191,18 @@
     , 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]
+    , MIP.varDomains        = Map.fromAscList
+       [ (v, (t, bs))
+       | v <- Set.toAscList vs
+       , let t =
+               if isInt v then
+                 if isSemi v then MIP.SemiIntegerVariable
+                 else MIP.IntegerVariable
+               else
+                 if isSemi v then MIP.SemiContinuousVariable
+                 else MIP.ContinuousVariable
+       , let bs = Map.findWithDefault MIP.defaultBounds v bnds2
+       ]
     }
 
 problem :: C e s m => m (MIP.ObjectiveFunction Scientific)
@@ -296,7 +274,7 @@
 indicator = do
   var <- variable
   tok (char '=')
-  val <- tok ((char '0' >> return 0) <|> (char '1' >> return 1))
+  val <- number  -- numbers other than 0 or 1 should be error?
   tok $ string "->"
   return (var, val)
 
@@ -458,11 +436,7 @@
        ]
 
 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
@@ -489,7 +463,7 @@
 render _ mip = Right $ execM $ render' $ normalize mip
 
 writeVar :: MIP.Var -> M ()
-writeVar v = writeString $ unintern v
+writeVar (MIP.Var v) = writeString v
 
 render' :: MIP.Problem Scientific -> M ()
 render' mip = do
@@ -589,13 +563,13 @@
 
     f :: MIP.Term Scientific -> T.Text
     f (MIP.Term c [])  = showConstTerm c
-    f (MIP.Term c [v]) = showCoeff c <> fromString (MIP.fromVar v)
+    f (MIP.Term c [v]) = showCoeff c <> MIP.varName 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))
+      mconcat (intersperse " * " (map MIP.varName vs))
 
 showValue :: Scientific -> T.Text
 showValue = fromString . show
@@ -633,7 +607,10 @@
     Just (v,vval) -> do
       writeVar v
       writeString " = "
-      tell $ B.scientificBuilder vval
+      tell $
+        case floatingOrInteger vval of
+          Right (i :: Integer) -> B.decimal i
+          Left (_ :: Double) -> B.scientificBuilder vval  -- should be error?
       writeString " -> "
 
   renderExpr False e
@@ -654,7 +631,7 @@
 renderBoundExpr MIP.PosInf = writeString "+inf"
 
 renderVariableList :: [MIP.Var] -> M ()
-renderVariableList vs = fill 80 (map unintern vs) >> writeChar '\n'
+renderVariableList vs = fill 80 (map MIP.varName vs) >> writeChar '\n'
 
 fill :: Int -> [T.Text] -> M ()
 fill width str = go str 0
@@ -692,7 +669,7 @@
   let gensym = do
         vs <- readSTRef vsRef
         let loop !c = do
-              let v = MIP.toVar ("~r_" ++ show c)
+              let v = fromString ("~r_" ++ show c)
               if v `Set.member` vs then
                 loop (c+1)
               else do
@@ -720,8 +697,7 @@
   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)
+    , MIP.varDomains = MIP.varDomains prob `Map.union` Map.fromList [(v, (MIP.ContinuousVariable, bs)) | (v,bs) <- newvs]
     }
 
 removeEmptyExpr :: Num r => MIP.Problem r -> MIP.Problem r
@@ -734,7 +710,7 @@
   where
     obj = MIP.objectiveFunction prob
 
-    convertExpr (MIP.Expr []) = MIP.Expr [MIP.Term 0 [MIP.toVar "x0"]]
+    convertExpr (MIP.Expr []) = MIP.Expr [MIP.Term 0 [fromString "x0"]]
     convertExpr e = e
 
     convertConstr constr =
diff --git a/src/Numeric/Optimization/MIP/MPSFile.hs b/src/Numeric/Optimization/MIP/MPSFile.hs
--- a/src/Numeric/Optimization/MIP/MPSFile.hs
+++ b/src/Numeric/Optimization/MIP/MPSFile.hs
@@ -1,11 +1,11 @@
 {-# OPTIONS_GHC -Wall -fno-warn-unused-do-bind #-}
 {-# OPTIONS_HADDOCK show-extensions #-}
-{-# LANGUAGE CPP #-}
 {-# LANGUAGE ConstraintKinds #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE OverloadedStrings #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.Optimization.MIP.MPSFile
@@ -37,17 +37,11 @@
   , 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)
@@ -62,17 +56,10 @@
 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 Text.Megaparsec.Char hiding (string', eol)
 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
@@ -97,19 +84,11 @@
 
 -- ---------------------------------------------------------------------------
 
-#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.
@@ -127,10 +106,8 @@
 -- ---------------------------------------------------------------------------
 
 
-#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']
@@ -144,49 +121,41 @@
 commentline :: C e s m => m ()
 commentline = do
   _ <- char '*'
-  _ <- manyTill anyChar P.newline
+  _ <- manyTill anyChar P.eol
   return ()
 
-newline' :: C e s m => m ()
-newline' = do
+eol' :: C e s m => m ()
+eol' = do
   spaces'
-  _ <- P.newline
+  _ <- P.eol
   skipMany commentline
   return ()
 
 tok :: C e s m => m a -> m a
 tok p = do
   x <- p
-  msum [eof, lookAhead (char '\n' >> return ()), spaces1']
+  msum [eof, lookAhead (P.eol >> 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
+column = liftM MIP.Var $ ident
 
 ident :: C e s m => m T.Text
-ident = liftM fromString $ tok $ some $ noneOf [' ', '\t', '\n']
+ident = liftM fromString $ tok $ some $ noneOf [' ', '\t', '\r', '\n']
 
 stringLn :: C e s m => String -> m ()
-stringLn s = string (fromString s) >> newline'
+stringLn s = string (fromString s) >> eol'
 
 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
 
@@ -353,21 +322,22 @@
         , 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]
+        , MIP.varDomains            = Map.fromAscList
+            [ (v, (t, bs))
+            | v <- Set.toAscList vs
+            , let t =
+                    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
+            , let bs = Map.findWithDefault MIP.defaultBounds v bounds
+            ]
         }
 
   return mip
@@ -378,7 +348,7 @@
   n <- optional $ try $ do
     spaces1'
     ident
-  newline'
+  eol'
   return n
 
 objSenseSection :: C e s m => m OptDir
@@ -394,7 +364,7 @@
   try $ stringLn "OBJNAME"
   spaces1'
   name <- ident
-  newline'
+  eol'
   return name
 
 rowsSection :: C e s m => m [(Maybe MIP.RelOp, Row)]
@@ -423,7 +393,7 @@
         ]
   spaces1'
   name <- row
-  newline'
+  eol'
   return (op, name)
 
 colsSection :: forall e s m. C e s m => m (Map Column (Map Row Scientific), Set Column)
@@ -452,15 +422,15 @@
       spaces1'
       b <-  (try (string "'INTORG'") >> return True)
         <|> (string "'INTEND'" >> return False)
-      newline'
+      eol'
       return b
 
     entry :: T.Text -> m (Column, Map Row Scientific)
     entry x = do
-      let col = intern x
+      let col = MIP.Var x
       rv1 <- rowAndVal
       opt <- optional rowAndVal
-      newline'
+      eol'
       case opt of
         Nothing -> return (col, rv1)
         Just rv2 ->  return (col, Map.union rv1 rv2)
@@ -481,7 +451,7 @@
       _name <- ident
       rv1 <- rowAndVal
       opt <- optional rowAndVal
-      newline'
+      eol'
       case opt of
         Nothing  -> return rv1
         Just rv2 -> return $ Map.union rv1 rv2
@@ -496,7 +466,7 @@
       _name <- ident
       rv1 <- rowAndVal
       opt <- optional rowAndVal
-      newline'
+      eol'
       case opt of
         Nothing  -> return rv1
         Just rv2 -> return $ Map.union rv1 rv2
@@ -514,7 +484,7 @@
       val   <- if typ `elem` [FR, BV, MI, PL]
                then return 0
                else number
-      newline'
+      eol'
       return (typ, col, val)
 
 boundType :: C e s m => m BoundType
@@ -532,7 +502,7 @@
           <|> (string "S2" >> return MIP.S2)
       spaces1'
       name <- ident
-      newline'
+      eol'
       xs <- many (try identAndVal)
       return $ MIP.SOSConstraint{ MIP.sosLabel = Just name, MIP.sosType = typ, MIP.sosBody = xs }
 
@@ -541,7 +511,7 @@
       spaces1'
       col <- column
       val <- number
-      newline'
+      eol'
       return (col, val)
 
 quadObjSection :: C e s m => m [MIP.Term Scientific]
@@ -554,7 +524,7 @@
       col1 <- column
       col2 <- column
       val  <- number
-      newline'
+      eol'
       return $ MIP.Term (if col1 /= col2 then val else val / 2) [col1, col2]
 
 qMatrixSection :: C e s m => m [MIP.Term Scientific]
@@ -567,7 +537,7 @@
       col1 <- column
       col2 <- column
       val  <- number
-      newline'
+      eol'
       return $ MIP.Term (val / 2) [col1, col2]
 
 qcMatrixSection :: C e s m => m (Row, [MIP.Term Scientific])
@@ -575,7 +545,7 @@
   try $ string "QCMATRIX"
   spaces1'
   r <- row
-  newline'
+  eol'
   xs <- many entry
   return (r, xs)
   where
@@ -584,7 +554,7 @@
       col1 <- column
       col2 <- column
       val  <- number
-      newline'
+      eol'
       return $ MIP.Term val [col1, col2]
 
 indicatorsSection :: C e s m => m (Map Row (Column, Scientific))
@@ -599,7 +569,7 @@
       r <- row
       var <- column
       val <- number
-      newline'
+      eol'
       return (r, (var, val))
 
 -- ---------------------------------------------------------------------------
@@ -620,10 +590,10 @@
 -- | Render a problem into a 'TL.Text' containing MPS file data.
 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 opt mip = Right $ execM $ render' opt $ nameRows mip
 
-render' :: MIP.Problem Scientific -> M ()
-render' mip = do
+render' :: MIP.FileOptions -> MIP.Problem Scientific -> M ()
+render' opt mip = do
   let probName = fromMaybe "" (MIP.name mip)
 
   -- NAME section
@@ -637,18 +607,18 @@
        } = 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"]
+  when (MIP.optMPSWriteObjSense opt == MIP.WriteAlways ||
+        MIP.optMPSWriteObjSense opt == MIP.WriteIfNotDefault && dir /= OptMin) $ do
+    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]
--}
+  when (MIP.optMPSWriteObjName opt) $ do
+    writeSectionHeader "OBJNAME"
+    writeFields [objName]
 
   let splitRange c =
         case (MIP.constrLB c, MIP.constrUB c) of
@@ -696,7 +666,7 @@
              ]
       f col xs =
         forM_ (Map.toList xs) $ \(r, d) -> do
-          writeFields ["", unintern col, r, showValue d]
+          writeFields ["", MIP.varName 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
@@ -719,33 +689,33 @@
 
   -- BOUNDS section
   writeSectionHeader "BOUNDS"
-  forM_ (Map.toList (MIP.varType mip)) $ \(col, vt) -> do
+  forM_ (Map.toList (MIP.varDomains 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]
+        writeFields ["FR", "bound", MIP.varName col]
 
       (MIP.Finite 0, MIP.Finite 1) | vt == MIP.IntegerVariable -> do
         -- variable is binary (equal 0 or 1)
-        writeFields ["BV", "bound", unintern col]
+        writeFields ["BV", "bound", MIP.varName col]
 
       (MIP.Finite a, MIP.Finite b) | a == b -> do
         -- variable is fixed at the specified value
-        writeFields ["FX", "bound", unintern col, showValue a]
+        writeFields ["FX", "bound", MIP.varName col, showValue a]
 
       _ -> do
         case lb of
           MIP.PosInf -> error "should not happen"
           MIP.NegInf -> do
             -- Minus infinity
-            writeFields ["MI", "bound", unintern col]
+            writeFields ["MI", "bound", MIP.varName 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]
+            writeFields [t, "bound", MIP.varName col, showValue a]
 
         case ub of
           MIP.NegInf -> error "should not happen"
@@ -753,7 +723,7 @@
           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
+            writeFields ["PL", "bound", MIP.varName col] -- Plus infinity
           MIP.Finite a -> do
             let t = case vt of
                       MIP.SemiContinuousVariable -> "SC" -- Upper bound for semi-continuous variable
@@ -762,7 +732,7 @@
                         "SC"
                       MIP.IntegerVariable -> "UI" -- Upper bound for integer variable
                       _ -> "UP" -- Upper bound
-            writeFields [t, "bound", unintern col, showValue a]
+            writeFields [t, "bound", MIP.varName col, showValue a]
 
   -- QMATRIX section
   -- Gurobiは対称行列になっていないと "qmatrix isn't symmetric" というエラーを発生させる
@@ -770,7 +740,7 @@
      unless (Map.null qm) $ do
        writeSectionHeader "QMATRIX"
        forM_ (Map.toList qm) $ \(((v1,v2), val)) -> do
-         writeFields ["", unintern v1, unintern v2, showValue val]
+         writeFields ["", MIP.varName v1, MIP.varName v2, showValue val]
 
   -- SOS section
   unless (null (MIP.sosConstraints mip)) $ do
@@ -781,7 +751,7 @@
                 MIP.S2 -> "S2"
       writeFields $ t : maybeToList (MIP.sosLabel sos)
       forM_ (MIP.sosBody sos) $ \(var,val) -> do
-        writeFields ["", unintern var, showValue val]
+        writeFields ["", MIP.varName var, showValue val]
 
   -- QCMATRIX section
   let xs = [ (fromJust $ MIP.constrLabel c, qm)
@@ -794,7 +764,7 @@
       -- 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]
+        writeFields ["", MIP.varName v1, MIP.varName v2, showValue val]
 
   -- INDICATORS section
   -- Note: Gurobi-5.6.3 does not support this section.
@@ -803,7 +773,7 @@
     writeSectionHeader "INDICATORS"
     forM_ ics $ \c -> do
       let Just (var,val) = MIP.constrIndicator c
-      writeFields ["IF", fromJust (MIP.constrLabel c), unintern var, showValue val]
+      writeFields ["IF", fromJust (MIP.constrLabel c), MIP.varName var, showValue val]
 
   -- ENDATA section
   writeSectionHeader "ENDATA"
diff --git a/src/Numeric/Optimization/MIP/Solution/CBC.hs b/src/Numeric/Optimization/MIP/Solution/CBC.hs
--- a/src/Numeric/Optimization/MIP/Solution/CBC.hs
+++ b/src/Numeric/Optimization/MIP/Solution/CBC.hs
@@ -1,6 +1,5 @@
 {-# OPTIONS_GHC -Wall #-}
 {-# OPTIONS_HADDOCK show-extensions #-}
-{-# LANGUAGE CPP #-}
 {-# LANGUAGE OverloadedStrings #-}
 -----------------------------------------------------------------------------
 -- |
@@ -20,13 +19,9 @@
   ) where
 
 import Prelude hiding (readFile, writeFile)
-#if !MIN_VERSION_base(4,8,0)
-import Control.Applicative
-#endif
 
 import Control.Monad (foldM)
 import Control.Monad.Except
-import Data.Interned (intern)
 import Data.Map (Map)
 import qualified Data.Map as Map
 import Data.Scientific (Scientific)
@@ -63,8 +58,8 @@
   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
+          ("**":_no:var:val:_) -> return $ (MIP.Var (TL.toStrict var), read (TL.unpack val)) : vs
+          (_no:var:val:_) -> return $ (MIP.Var (TL.toStrict var), read (TL.unpack val)) : vs
           [] -> return $ vs
           _ -> throwError ("Numeric.Optimization.MIP.Solution.CBC: invalid line " ++ show t)
   vs <- foldM f [] ls
diff --git a/src/Numeric/Optimization/MIP/Solution/CPLEX.hs b/src/Numeric/Optimization/MIP/Solution/CPLEX.hs
--- a/src/Numeric/Optimization/MIP/Solution/CPLEX.hs
+++ b/src/Numeric/Optimization/MIP/Solution/CPLEX.hs
@@ -1,6 +1,5 @@
 {-# OPTIONS_GHC -Wall #-}
 {-# OPTIONS_HADDOCK show-extensions #-}
-{-# LANGUAGE CPP #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE OverloadedStrings #-}
 -----------------------------------------------------------------------------
@@ -21,11 +20,7 @@
   ) 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
@@ -69,7 +64,7 @@
           let m = XML.elementAttributes e
           name <- Map.lookup "name" m
           value <- read . T.unpack <$> Map.lookup "value" m
-          return (intern name, value)
+          return (MIP.Var name, value)
       | otherwise = []
     vs = fromDocument doc
       $| element "CPLEXSolution"
diff --git a/src/Numeric/Optimization/MIP/Solution/GLPK.hs b/src/Numeric/Optimization/MIP/Solution/GLPK.hs
--- a/src/Numeric/Optimization/MIP/Solution/GLPK.hs
+++ b/src/Numeric/Optimization/MIP/Solution/GLPK.hs
@@ -1,6 +1,5 @@
 {-# OPTIONS_GHC -Wall #-}
 {-# OPTIONS_HADDOCK show-extensions #-}
-{-# LANGUAGE CPP #-}
 {-# LANGUAGE OverloadedStrings #-}
 -----------------------------------------------------------------------------
 -- |
@@ -20,10 +19,6 @@
   ) 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)
@@ -83,8 +78,8 @@
     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
+        (_no : col : "*" : activity : _) -> f ((MIP.Var (TL.toStrict col), read (TL.unpack activity)) : ret) ls3
+        (_no : col : activity : _) -> f ((MIP.Var (TL.toStrict col), read (TL.unpack activity)) : ret) ls3
         _ -> error "parse error"
       where
         (ws,ls3) =
diff --git a/src/Numeric/Optimization/MIP/Solution/Gurobi.hs b/src/Numeric/Optimization/MIP/Solution/Gurobi.hs
--- a/src/Numeric/Optimization/MIP/Solution/Gurobi.hs
+++ b/src/Numeric/Optimization/MIP/Solution/Gurobi.hs
@@ -22,16 +22,11 @@
   ) 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)
+#if !MIN_VERSION_base(4,20,0)
 import Data.List (foldl')
-import qualified Data.Map as Map
-#if !MIN_VERSION_base(4,11,0)
-import Data.Monoid
 #endif
+import qualified Data.Map as Map
 import Data.Scientific (Scientific)
 import qualified Data.Text.Lazy as TL
 import qualified Data.Text.Lazy.Builder as B
@@ -46,7 +41,7 @@
     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'
+    ls2 = [ B.fromText (MIP.varName name) <> B.singleton ' ' <> B.scientificBuilder val <> B.singleton '\n'
           | (name,val) <- Map.toList (MIP.solVariables sol)
           ]
 
@@ -71,7 +66,7 @@
           (Just r, vs)
       | otherwise =
           case TL.words (TL.takeWhile (/= '#') l) of
-            [w1, w2] -> (obj, (intern (TL.toStrict w1), read (TL.unpack w2)) : vs)
+            [w1, w2] -> (obj, (MIP.Var (TL.toStrict w1), read (TL.unpack w2)) : vs)
             [] -> (obj, vs)
             _ -> error ("Numeric.Optimization.MIP.Solution.Gurobi: invalid line " ++ show l)
 
diff --git a/src/Numeric/Optimization/MIP/Solution/HiGHS.hs b/src/Numeric/Optimization/MIP/Solution/HiGHS.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Optimization/MIP/Solution/HiGHS.hs
@@ -0,0 +1,133 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# OPTIONS_HADDOCK show-extensions #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE TypeOperators #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Optimization.MIP.Solution.HiGHS
+-- Copyright   :  (c) Masahiro Sakai 2025
+-- License     :  BSD-style
+--
+-- Maintainer  :  masahiro.sakai@gmail.com
+-- Stability   :  provisional
+-- Portability :  non-portable
+--
+-----------------------------------------------------------------------------
+module Numeric.Optimization.MIP.Solution.HiGHS
+  ( Solution (..)
+  , parse
+  , readFile
+  ) where
+
+import Prelude hiding (readFile, writeFile)
+
+import Control.Exception
+import Control.Monad
+import Data.Char (isSpace)
+import qualified Data.Map as Map
+import Data.Scientific (Scientific)
+import Data.String
+import qualified Data.Text.Lazy as TL
+import qualified Data.Text.Lazy.IO as TLIO
+import System.IO hiding (readFile, writeFile)
+import qualified Text.Megaparsec as Megaparsec
+import Text.Megaparsec hiding (label, skipManyTill, parse, ParseError)
+import Text.Megaparsec.Char hiding (string', char')
+import qualified Text.Megaparsec.Char.Lexer as P
+
+import qualified Numeric.Optimization.MIP.Base as MIP
+import Numeric.Optimization.MIP.Base (Solution (..))
+import Numeric.Optimization.MIP.FileUtils (ParseError)
+
+parser :: forall e s m. (MonadParsec e s m, Token s ~ Char, IsString (Tokens s)) => m (MIP.Solution Scientific)
+parser = do
+  _ <- strippedLine (string "Model status")
+  status <- stripped $ msum $ map Megaparsec.try
+    [ string "Not Set" *> pure MIP.StatusUnknown
+    , string "Load error" *> pure MIP.StatusUnknown
+    , string "Model error" *> pure MIP.StatusUnknown
+    , string "Presolve error" *> pure MIP.StatusUnknown
+    , string "Solve error" *> pure MIP.StatusUnknown
+    , string "Postsolve error" *> pure MIP.StatusUnknown
+    , string "Empty" *> pure MIP.StatusUnknown
+    , string "Memory limit reached" *> pure MIP.StatusUnknown
+    , string "Optimal" *> pure MIP.StatusOptimal
+    , string "Infeasible" *> pure MIP.StatusInfeasible
+    , string "Primal infeasible or unbounded" *> pure MIP.StatusInfeasibleOrUnbounded
+    , string "Unbounded" *> pure MIP.StatusUnbounded
+    , string "Bound on objective reached" *> pure MIP.StatusFeasible -- ???
+    , string "Target for objective reached" *> pure MIP.StatusFeasible -- ???
+    , string "Time limit reached" *> pure MIP.StatusUnknown
+    , string "Iteration limit reached" *> pure MIP.StatusUnknown
+    , string "Solution limit reached" *> pure MIP.StatusFeasible -- ???
+    , string "Interrupted by user" *> pure MIP.StatusUnknown
+    , string "Unknown" *> pure MIP.StatusUnknown
+    ]
+  _ <- eol
+
+  space
+
+  _ <- strippedLine (string "# Primal solution values")
+  hspace
+  msum
+    [ do
+        _ <- string "None" <* hspace <* eol
+        return $
+          MIP.Solution
+          { MIP.solStatus = status
+          , MIP.solObjectiveValue = Nothing
+          , MIP.solVariables = Map.empty
+          }
+    , do
+        _ <- (string "Feasible" <|> string "Infeasible") <* hspace <* eol
+        obj <- strippedLine (string "Objective" *> hspace1 *> P.scientific)
+        n <- strippedLine (string "# Columns " *> P.decimal)
+        values <- replicateM n $ strippedLine ((,) <$> (ident <* hspace1) <*> P.scientific)
+        return $
+          MIP.Solution
+          { MIP.solStatus = status
+          , MIP.solObjectiveValue = Just obj
+          , MIP.solVariables = Map.fromList values
+          }
+    ]
+
+  where
+    stripped :: m a -> m a
+    stripped p = hspace *> p <* hspace
+
+    strippedLine :: m a -> m a
+    strippedLine p = stripped p <* eol
+
+    ident :: m MIP.Var
+    ident = fromString <$> some (satisfy (not . isSpace)) <?> "identifier"
+
+parse :: TL.Text -> Either (ParseError TL.Text) (MIP.Solution Scientific)
+parse = Megaparsec.parse parser "<string>"
+
+readFile :: FilePath -> IO (MIP.Solution Scientific)
+readFile fname = do
+  h <- openFile fname ReadMode
+  hSetEncoding h utf8
+  ret <- Megaparsec.parse parser fname <$> TLIO.hGetContents h
+  case ret of
+    Left e -> throwIO (e :: ParseError TL.Text)
+    Right a -> return a
+
+#if !MIN_VERSION_megaparsec(9,0,0)
+
+hspace :: (MonadParsec e s m, Token s ~ Char) => m ()
+hspace = void $ takeWhileP (Just "white space") isHSpace
+{-# INLINE hspace #-}
+
+hspace1 :: (MonadParsec e s m, Token s ~ Char) => m ()
+hspace1 = void $ takeWhile1P (Just "white space") isHSpace
+{-# INLINE hspace1 #-}
+
+isHSpace :: Char -> Bool
+isHSpace x = isSpace x && x /= '\n' && x /= '\r'
+
+#endif
diff --git a/src/Numeric/Optimization/MIP/Solution/Printemps.hs b/src/Numeric/Optimization/MIP/Solution/Printemps.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Optimization/MIP/Solution/Printemps.hs
@@ -0,0 +1,71 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# OPTIONS_HADDOCK show-extensions #-}
+{-# LANGUAGE DeriveGeneric #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Optimization.MIP.Solution.Printemps
+-- Copyright   :  (c) Masahiro Sakai 2025
+-- License     :  BSD-style
+--
+-- Maintainer  :  masahiro.sakai@gmail.com
+-- Stability   :  provisional
+-- Portability :  non-portable
+--
+-----------------------------------------------------------------------------
+module Numeric.Optimization.MIP.Solution.Printemps
+  ( Solution (..)
+  , readFile
+  ) where
+
+import Prelude hiding (readFile, writeFile)
+
+import qualified Data.Aeson as J
+import Data.Map (Map)
+import qualified Data.Map as Map
+import Data.Scientific (Scientific)
+import qualified Data.Text as T
+import GHC.Generics
+
+import qualified Numeric.Optimization.MIP.Base as MIP
+import Numeric.Optimization.MIP.Base (Solution (..))
+
+data Incumbent
+  = Incumbent
+  { incumbentVersion :: T.Text
+  , incumbentName :: T.Text
+  , incumbentNumberOfVariables :: Int
+  , incumbentIsFoundFeasibleSolution :: Bool
+  , incumbentObjective :: Scientific
+  , incumbentTotalViolation :: Scientific
+  , incumbentVariables :: Map T.Text Scientific
+  , incumbentExpressions :: Map T.Text Scientific
+  , incumbentConstraints :: Map T.Text Scientific
+  , incumbentViolations :: Map T.Text Scientific
+  }
+  deriving (Generic, Show)
+
+customOptions :: J.Options
+customOptions =
+  J.defaultOptions
+  { J.fieldLabelModifier = (J.camelTo2 '_' . drop (length "Incumbent"))
+  }
+
+instance J.FromJSON Incumbent where
+  parseJSON = J.genericParseJSON customOptions
+
+readFile :: FilePath -> IO (MIP.Solution Scientific)
+readFile fname = do
+  ret <- J.eitherDecodeFileStrict' fname
+  case ret of
+    Left err -> ioError $ userError err
+    Right incumbent -> return $
+      MIP.Solution
+      { MIP.solStatus =
+          if incumbentIsFoundFeasibleSolution incumbent
+          then MIP.StatusFeasible
+          else MIP.StatusUnknown
+      , MIP.solObjectiveValue =
+          Just (incumbentObjective incumbent)
+      , MIP.solVariables =
+          Map.fromList [(MIP.Var v, val) | (v, val) <- Map.toList (incumbentVariables incumbent)]
+      }
diff --git a/src/Numeric/Optimization/MIP/Solution/SCIP.hs b/src/Numeric/Optimization/MIP/Solution/SCIP.hs
--- a/src/Numeric/Optimization/MIP/Solution/SCIP.hs
+++ b/src/Numeric/Optimization/MIP/Solution/SCIP.hs
@@ -1,6 +1,5 @@
 {-# OPTIONS_GHC -Wall #-}
 {-# OPTIONS_HADDOCK show-extensions #-}
-{-# LANGUAGE CPP #-}
 {-# LANGUAGE OverloadedStrings #-}
 -----------------------------------------------------------------------------
 -- |
@@ -20,12 +19,8 @@
   ) where
 
 import Prelude hiding (readFile, writeFile)
-#if !MIN_VERSION_base(4,8,0)
-import Control.Applicative
-#endif
 import Control.Monad (foldM)
 import Control.Monad.Except
-import Data.Interned (intern)
 import Data.Map (Map)
 import qualified Data.Map as Map
 import Data.Scientific (Scientific)
@@ -62,7 +57,7 @@
     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
+            (w1:w2:_) -> return $ (MIP.Var (TL.toStrict w1), read (TL.unpack w2)) : vs
             [] -> return $ vs
             _ -> throwError ("Numeric.Optimization.MIP.Solution.SCIP: invalid line " ++ show t)
     vs <- foldM f [] ts
diff --git a/src/Numeric/Optimization/MIP/Solver.hs b/src/Numeric/Optimization/MIP/Solver.hs
--- a/src/Numeric/Optimization/MIP/Solver.hs
+++ b/src/Numeric/Optimization/MIP/Solver.hs
@@ -11,19 +11,29 @@
 --
 -----------------------------------------------------------------------------
 module Numeric.Optimization.MIP.Solver
-  ( module Numeric.Optimization.MIP.Solver.Base
+  (
+  -- * Solver type
+    IsSolver (..)
+  , SolveOptions (..)
+  -- * Individual solvers
   , 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.HiGHS
   , module Numeric.Optimization.MIP.Solver.LPSolve
+  , module Numeric.Optimization.MIP.Solver.Printemps
   , module Numeric.Optimization.MIP.Solver.SCIP
+  -- * Utilities
+  , Default (..)
   ) where
 
-import Numeric.Optimization.MIP.Solver.Base
+import Numeric.Optimization.MIP.Solver.Base hiding (solve')
 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.HiGHS
 import Numeric.Optimization.MIP.Solver.LPSolve
+import Numeric.Optimization.MIP.Solver.Printemps
 import Numeric.Optimization.MIP.Solver.SCIP
diff --git a/src/Numeric/Optimization/MIP/Solver/Base.hs b/src/Numeric/Optimization/MIP/Solver/Base.hs
--- a/src/Numeric/Optimization/MIP/Solver/Base.hs
+++ b/src/Numeric/Optimization/MIP/Solver/Base.hs
@@ -13,32 +13,68 @@
 --
 -----------------------------------------------------------------------------
 module Numeric.Optimization.MIP.Solver.Base
-  ( SolveOptions (..)
+  (
+  -- * Solver type
+    IsSolver (..)
+  , SolveOptions (..)
+  -- * Utilities
   , Default (..)
-  , IsSolver (..)
   ) where
 
 import Data.Default.Class
 import Data.Scientific (Scientific)
 import Numeric.Optimization.MIP.Base as MIP
+import qualified Data.Map as Map
 
+-- | Options for 'solve' function
 data SolveOptions
   = SolveOptions
   { solveTimeLimit :: Maybe Double
     -- ^ time limit in seconds
+  , solveTol :: Maybe (MIP.Tol Scientific)
+    -- ^ tolerance
   , solveLogger :: String -> IO ()
     -- ^ invoked when a solver output a line
   , solveErrorLogger :: String -> IO ()
     -- ^ invoked when a solver output a line to stderr
+  , solveCondensedSolution :: Bool
+    -- ^ potentially omit variables set to zero from the solution
   }
 
 instance Default SolveOptions where
   def =
     SolveOptions
     { solveTimeLimit = Nothing
+    , solveTol = Nothing
     , solveLogger = const $ return ()
     , solveErrorLogger = const $ return ()
+    , solveCondensedSolution = False
     }
 
+
+-- | Type class for solvers
+--
+-- 
 class Monad m => IsSolver s m | s -> m where
-  solve :: s -> SolveOptions -> MIP.Problem Scientific -> m (MIP.Solution Scientific)
+  -- | Low level version of 'solve'' that allows omission of variables with a value 0.
+  --
+  -- Implementor of the type class must implement this method.
+  solve' :: s -> SolveOptions -> MIP.Problem Scientific -> m (MIP.Solution Scientific)
+
+  -- | A method for solving 'MIP.Problem'
+  --
+  -- This method is a bit higher level than 'solve'' in that it does not omit variables
+  -- with a value @0@ unless 'solveCondensedSolution' is set to @True@.
+  -- Implementor of the type class can override this method as @solve = solve'@ if the
+  -- solver always returns all variables.
+  solve  :: s -> SolveOptions -> MIP.Problem Scientific -> m (MIP.Solution Scientific)
+  solve s opts problem = (if solveCondensedSolution opts then id else addZeroes problem) <$> solve' s opts problem
+  {-# MINIMAL solve' #-}
+
+-- Several solvers (at least CBC) do not include any variables set to 0 in their solution.
+-- TODO: for solvers that do return all variables, add @solve = solve'@
+-- for a minor performance improvement.
+addZeroes :: MIP.Problem Scientific -> MIP.Solution Scientific -> MIP.Solution Scientific
+addZeroes problem (Solution stat obj solmap) = 
+  -- Map.union is left-biased: only values not present in the solution are added.
+  Solution stat obj $ Map.union solmap (Map.fromSet (const 0) (vars problem))
diff --git a/src/Numeric/Optimization/MIP/Solver/CBC.hs b/src/Numeric/Optimization/MIP/Solver/CBC.hs
--- a/src/Numeric/Optimization/MIP/Solver/CBC.hs
+++ b/src/Numeric/Optimization/MIP/Solver/CBC.hs
@@ -28,19 +28,24 @@
 import qualified Numeric.Optimization.MIP.Solution.CBC as CBCSol
 import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback)
 
+-- | A solver instance for calling @cbc@ command from [CBC (COIN-OR Branch-and-Cut solver)](https://github.com/coin-or/Cbc).
+--
+-- Use 'cbc' and record update syntax to modify its field.
 data CBC
   = CBC
   { cbcPath :: String
+  , cbcArgs :: [String]
   }
 
 instance Default CBC where
   def = cbc
 
+-- | Default value of t'CBC'
 cbc :: CBC
-cbc = CBC "cbc"
+cbc = CBC "cbc" []
 
 instance IsSolver CBC IO where
-  solve solver opt prob = do
+  solve' solver opt prob = do
     case LPFile.render def prob{ MIP.objectiveFunction = obj' } of
       Left err -> ioError $ userError err
       Right lp -> do
@@ -49,14 +54,22 @@
           hClose h1
           withSystemTempFile "cbc.sol" $ \fname2 h2 -> do
             hClose h2
-            let args = [fname1]
+            let args = cbcArgs solver
+                    ++ [fname1]
                     ++ (case solveTimeLimit opt of
                           Nothing -> []
                           Just sec -> ["sec", show sec])
+                    ++ (case solveTol opt of
+                          Nothing -> []
+                          Just tol ->
+                            [ "integerTolerance", show (MIP.integralityTol tol)
+                            , "primalTolerance", show (MIP.feasibilityTol tol)
+                            , "dualTolerance", show (MIP.optimalityTol tol)
+                            ])
                     ++ ["solve", "solu", fname2]
                 onGetLine = solveLogger opt
                 onGetErrorLine = solveErrorLogger opt
-            exitcode <- runProcessWithOutputCallback (cbcPath solver) args "" onGetLine onGetErrorLine
+            exitcode <- runProcessWithOutputCallback (cbcPath solver) args Nothing "" onGetLine onGetErrorLine
             case exitcode of
               ExitFailure n -> ioError $ userError $ "exit with " ++ show n
               ExitSuccess -> do
diff --git a/src/Numeric/Optimization/MIP/Solver/CPLEX.hs b/src/Numeric/Optimization/MIP/Solver/CPLEX.hs
--- a/src/Numeric/Optimization/MIP/Solver/CPLEX.hs
+++ b/src/Numeric/Optimization/MIP/Solver/CPLEX.hs
@@ -30,19 +30,25 @@
 import qualified Numeric.Optimization.MIP.Solution.CPLEX as CPLEXSol
 import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback)
 
+-- | A solver instance for calling @cplex@ command from [IBM ILOG CPLEX Optimization Studio](https://www.ibm.com/products/ilog-cplex-optimization-studio).
+--
+-- Use 'cplex' and record update syntax to modify its fields.
 data CPLEX
   = CPLEX
   { cplexPath :: String
+  , cplexArgs :: [String]                 
+  , cplexCommands :: [String]
   }
 
 instance Default CPLEX where
   def = cplex
 
+-- | Default value of t'CPLEX'.
 cplex :: CPLEX
-cplex = CPLEX "cplex"
+cplex = CPLEX "cplex" [] []
 
 instance IsSolver CPLEX IO where
-  solve solver opt prob = do
+  solve' solver opt prob = do
     case LPFile.render def prob of
       Left err -> ioError $ userError err
       Right lp -> do
@@ -52,13 +58,21 @@
           withSystemTempFile "cplex.sol" $ \fname2 h2 -> do
             hClose h2
             isInfeasibleRef <- newIORef False
-            let args = []
-                input = unlines $
+            let input = unlines $
                   (case solveTimeLimit opt of
-                          Nothing -> []
-                          Just sec -> ["set timelimit ", show sec]) ++
-                  [ "read " ++ show fname1
-                  , "optimize"
+                     Nothing -> []
+                     Just sec -> ["set timelimit " ++ show sec]) ++
+                  (case solveTol opt of
+                     Nothing -> []
+                     Just tol ->
+                       [ "set mip tolerances integrality " ++ show (MIP.integralityTol tol)
+                       , "set simplex tolerances feasibility " ++ show (MIP.feasibilityTol tol)
+                       , "set simplex tolerances optimality " ++ show (MIP.optimalityTol tol)
+                       ]
+                  ) ++
+                  [ "read " ++ show fname1 ] ++
+                  cplexCommands solver ++
+                  [ "optimize"
                   , "write " ++ show fname2
                   , "y"
                   , "quit"
@@ -68,7 +82,7 @@
                     writeIORef isInfeasibleRef True
                   solveLogger opt s
                 onGetErrorLine = solveErrorLogger opt
-            exitcode <- runProcessWithOutputCallback (cplexPath solver) args input onGetLine onGetErrorLine
+            exitcode <- runProcessWithOutputCallback (cplexPath solver) (cplexArgs solver) Nothing input onGetLine onGetErrorLine
             case exitcode of
               ExitFailure n -> ioError $ userError $ "exit with " ++ show n
               ExitSuccess -> do
diff --git a/src/Numeric/Optimization/MIP/Solver/Glpsol.hs b/src/Numeric/Optimization/MIP/Solver/Glpsol.hs
--- a/src/Numeric/Optimization/MIP/Solver/Glpsol.hs
+++ b/src/Numeric/Optimization/MIP/Solver/Glpsol.hs
@@ -30,19 +30,24 @@
 import qualified Numeric.Optimization.MIP.Solution.GLPK as GLPKSol
 import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback)
 
+-- | A solver instance for calling @glpsol@ command from [GLPK (GNU Linear Programming Kit)](https://www.gnu.org/software/glpk/).
+--
+-- Use 'glpsol' and record update syntax to modify its field.
 data Glpsol
   = Glpsol
   { glpsolPath :: String
+  , glpsolArgs :: [String]
   }
 
 instance Default Glpsol where
   def = glpsol
 
+-- | Default value of t'Glpsol'
 glpsol :: Glpsol
-glpsol = Glpsol "glpsol"
+glpsol = Glpsol "glpsol" []
 
 instance IsSolver Glpsol IO where
-  solve solver opt prob = do
+  solve' solver opt prob = do
     case LPFile.render def prob of
       Left err -> ioError $ userError err
       Right lp -> do
@@ -53,7 +58,8 @@
             hClose h2
             isUnboundedRef <- newIORef False
             isInfeasibleRef <- newIORef False
-            let args = ["--lp", fname1, "-o", fname2] ++
+            let args = glpsolArgs solver ++
+                       ["--lp", fname1, "-o", fname2] ++
                        (case solveTimeLimit opt of
                           Nothing -> []
                           Just sec -> ["--tmlim", show (max 1 (floor sec) :: Int)])
@@ -68,7 +74,7 @@
                     _ -> return ()
                   solveLogger opt s
                 onGetErrorLine = solveErrorLogger opt
-            exitcode <- runProcessWithOutputCallback (glpsolPath solver) args "" onGetLine onGetErrorLine
+            exitcode <- runProcessWithOutputCallback (glpsolPath solver) args Nothing "" onGetLine onGetErrorLine
             case exitcode of
               ExitFailure n -> ioError $ userError $ "exit with " ++ show n
               ExitSuccess -> do
diff --git a/src/Numeric/Optimization/MIP/Solver/GurobiCl.hs b/src/Numeric/Optimization/MIP/Solver/GurobiCl.hs
--- a/src/Numeric/Optimization/MIP/Solver/GurobiCl.hs
+++ b/src/Numeric/Optimization/MIP/Solver/GurobiCl.hs
@@ -30,19 +30,24 @@
 import qualified Numeric.Optimization.MIP.Solution.Gurobi as GurobiSol
 import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback)
 
+-- | A solver instance for calling @gurobi_cl@ command from [Gurobi Optimizer](https://www.gurobi.com/solutions/gurobi-optimizer/).
+--
+-- Use 'gurobiCl' and record update syntax to modify its field.
 data GurobiCl
   = GurobiCl
   { gurobiClPath :: String
+  , gurobiClArgs :: [String]
   }
 
 instance Default GurobiCl where
   def = gurobiCl
 
+-- | Default value of t'GurobiCl'
 gurobiCl :: GurobiCl
-gurobiCl = GurobiCl "gurobi_cl"
+gurobiCl = GurobiCl "gurobi_cl" []
 
 instance IsSolver GurobiCl IO where
-  solve solver opt prob = do
+  solve' solver opt prob = do
     case LPFile.render def prob of
       Left err -> ioError $ userError err
       Right lp -> do
@@ -52,10 +57,18 @@
           withSystemTempFile "gurobi.sol" $ \fname2 h2 -> do
             hClose h2
             statusRef <- newIORef MIP.StatusUnknown
-            let args = ["ResultFile=" ++ fname2]
+            let args = gurobiClArgs solver
+                    ++ ["ResultFile=" ++ fname2]
                     ++ (case solveTimeLimit opt of
                           Nothing -> []
                           Just sec -> ["TimeLimit=" ++ show sec])
+                    ++ (case solveTol opt of
+                          Nothing -> []
+                          Just tol ->
+                            [ "IntFeasTol=" ++ show (MIP.integralityTol tol)
+                            , "FeasibilityTol=" ++ show (MIP.feasibilityTol tol)
+                            , "OptimalityTol=" ++ show (MIP.optimalityTol tol)
+                            ])
                     ++ [fname1]
                 onGetLine s = do
                   case s of
@@ -67,7 +80,7 @@
                     _ -> return ()
                   solveLogger opt s
                 onGetErrorLine = solveErrorLogger opt
-            exitcode <- runProcessWithOutputCallback (gurobiClPath solver) args "" onGetLine onGetErrorLine
+            exitcode <- runProcessWithOutputCallback (gurobiClPath solver) args Nothing "" onGetLine onGetErrorLine
             case exitcode of
               ExitFailure n -> ioError $ userError $ "exit with " ++ show n
               ExitSuccess -> do
diff --git a/src/Numeric/Optimization/MIP/Solver/HiGHS.hs b/src/Numeric/Optimization/MIP/Solver/HiGHS.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Optimization/MIP/Solver/HiGHS.hs
@@ -0,0 +1,79 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# OPTIONS_HADDOCK show-extensions #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Optimization.MIP.Solver.HiGHS
+-- Copyright   :  (c) Masahiro Sakai 2025
+-- License     :  BSD-style
+--
+-- Maintainer  :  masahiro.sakai@gmail.com
+-- Stability   :  provisional
+-- Portability :  non-portable
+--
+-----------------------------------------------------------------------------
+module Numeric.Optimization.MIP.Solver.HiGHS
+  ( HiGHS (..)
+  , highs
+  ) where
+
+import Data.Default.Class
+import qualified Data.Text.Lazy.IO as TLIO
+import System.IO
+import System.IO.Temp
+import qualified Numeric.Optimization.MIP.LPFile as LPFile
+import qualified Numeric.Optimization.MIP.Base as MIP
+import Numeric.Optimization.MIP.Solver.Base
+import qualified Numeric.Optimization.MIP.Solution.HiGHS as HiGHSSol
+import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback)
+
+-- | A solver instance for calling @highs@ command from [HiGHS](https://github.com/ERGO-Code/HiGHS).
+--
+-- Use 'highs' and record update syntax to modify its field.
+data HiGHS
+  = HiGHS
+  { highsPath :: String
+  , highsArgs :: [String]
+  }
+
+instance Default HiGHS where
+  def = highs
+
+-- | Default value of t'HiGHS'
+highs :: HiGHS
+highs = HiGHS "highs" []
+
+instance IsSolver HiGHS IO where
+  solve' solver opt prob = do
+    case LPFile.render def prob of
+      Left err -> ioError $ userError err
+      Right lp -> do
+        withSystemTempFile "highs.lp" $ \fname1 h1 -> do
+          TLIO.hPutStr h1 lp
+          hClose h1
+
+          withSystemTempFile "highs.sol" $ \fname2 h2 -> do
+            hClose h2
+
+            withSystemTempFile "highs_options.txt" $ \fname3 h3 -> do
+              -- XXX: HiGHS does not support multiple options files
+              options_args <-
+                case solveTol opt of
+                  Nothing -> return []
+                  Just tol -> do
+                    hPutStrLn h3 $ "mip_feasibility_tolerance = " ++ show (MIP.integralityTol tol)
+                    hPutStrLn h3 $ "primal_feasibility_tolerance = " ++ show (MIP.feasibilityTol tol)
+                    hPutStrLn h3 $ "dual_feasibility_tolerance = " ++ show (MIP.optimalityTol tol)
+                    return ["--options_file", fname3]
+              hClose h3
+
+              let args = highsArgs solver ++
+                         ["--model_file", fname1, "--solution_file", fname2] ++
+                         options_args ++
+                         (case solveTimeLimit opt of
+                            Nothing -> []
+                            Just sec -> ["--time_limit", show sec])
+                  onGetLine s = solveLogger opt s
+                  onGetErrorLine = solveErrorLogger opt
+              _exitcode <- runProcessWithOutputCallback (highsPath solver) args Nothing "" onGetLine onGetErrorLine
+              HiGHSSol.readFile fname2
diff --git a/src/Numeric/Optimization/MIP/Solver/LPSolve.hs b/src/Numeric/Optimization/MIP/Solver/LPSolve.hs
--- a/src/Numeric/Optimization/MIP/Solver/LPSolve.hs
+++ b/src/Numeric/Optimization/MIP/Solver/LPSolve.hs
@@ -32,19 +32,24 @@
 import Numeric.Optimization.MIP.Solver.Base
 import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback)
 
+-- | A solver instance for calling @lp_solve@ command from [lp_solve](https://lpsolve.sourceforge.net/5.5/) package.
+--
+-- Use 'lpSolve' and record update syntax to modify its field.
 data LPSolve
   = LPSolve
   { lpSolvePath :: String
+  , lpSolveArgs :: [String]
   }
 
 instance Default LPSolve where
   def = lpSolve
 
+-- | Default value of t'LPSolve'
 lpSolve :: LPSolve
-lpSolve = LPSolve "lp_solve"
+lpSolve = LPSolve "lp_solve" []
 
 instance IsSolver LPSolve IO where
-  solve solver opt prob = do
+  solve' solver opt prob = do
     case MPSFile.render def prob of
       Left err -> ioError $ userError err
       Right lp -> do
@@ -54,9 +59,17 @@
           objRef <- newIORef Nothing
           solRef <- newIORef []
           flagRef <- newIORef False
-          let args = (case solveTimeLimit opt of
+          let args = lpSolveArgs solver
+                  ++ (case solveTimeLimit opt of
                         Nothing -> []
                         Just sec -> ["-timeout", show sec])
+                  ++ (case solveTol opt of
+                        Nothing -> []
+                        Just tol ->
+                          [ "-e", show (MIP.integralityTol tol)
+                          , "-epsb", show (MIP.feasibilityTol tol)
+                          , "-epsd", show (MIP.optimalityTol tol)
+                          ])
                   ++ ["-fmps", fname1]
               onGetLine s = do
                 case s of
@@ -72,7 +85,7 @@
                     return ()
                 solveLogger opt s
               onGetErrorLine = solveErrorLogger opt
-          exitcode <- runProcessWithOutputCallback (lpSolvePath solver) args "" onGetLine onGetErrorLine
+          exitcode <- runProcessWithOutputCallback (lpSolvePath solver) args Nothing "" onGetLine onGetErrorLine
           status <-
             case exitcode of
               ExitSuccess      -> return MIP.StatusOptimal
diff --git a/src/Numeric/Optimization/MIP/Solver/Printemps.hs b/src/Numeric/Optimization/MIP/Solver/Printemps.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Optimization/MIP/Solver/Printemps.hs
@@ -0,0 +1,107 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# OPTIONS_HADDOCK show-extensions #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE OverloadedStrings #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Optimization.MIP.Solver.Printemps
+-- Copyright   :  (c) Masahiro Sakai 2025
+-- License     :  BSD-style
+--
+-- Maintainer  :  masahiro.sakai@gmail.com
+-- Stability   :  provisional
+-- Portability :  non-portable
+--
+-----------------------------------------------------------------------------
+module Numeric.Optimization.MIP.Solver.Printemps
+  ( Printemps (..)
+  , printemps
+  ) where
+
+import qualified Data.Aeson as J
+import Data.Default.Class
+import Data.Map.Lazy (Map)
+import qualified Data.Map.Lazy as Map
+import qualified Data.Text as T
+import qualified Data.Text.Lazy.IO as TLIO
+import System.IO.Temp
+import qualified Numeric.Optimization.MIP.MPSFile as MPSFile
+import Numeric.Optimization.MIP.Base
+import Numeric.Optimization.MIP.Solver.Base
+import qualified Numeric.Optimization.MIP.Solution.Printemps as PrintempsSol
+import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback)
+import System.Exit
+import System.FilePath ((</>))
+
+-- | A solver instance for calling @mps_solver.exe@ command from [PRINTEMPS](https://snowberryfield.github.io/printemps/).
+--
+-- It requires PRINTEMPS version 2.6.0 or later.
+--
+-- Use 'printemps' and record update syntax to modify its field.
+data Printemps
+  = Printemps
+  { printempsPath :: String
+  , printempsArgs :: [String]
+  }
+
+instance Default Printemps where
+  def = printemps
+
+-- | Default value of t'Printemps'
+printemps :: Printemps
+printemps = Printemps "mps_solver" []
+
+instance IsSolver Printemps IO where
+  solve' solver opt prob = do
+    let prob' = prob
+                { name =
+                    case name prob of
+                      Just s | not (T.null s) -> Just s
+                      _ -> Just "problem"
+                }
+
+    let (orig_option_file, args') = removeOptionArgs (printempsArgs solver)
+    orig_option <-
+      case orig_option_file of
+        Nothing -> return Map.empty
+        Just fname -> do
+          ret <- J.eitherDecodeFileStrict' fname
+          case ret of
+            Left err -> ioError $ userError err
+            Right option -> return option
+    orig_general <-
+      case J.fromJSON (Map.findWithDefault (J.object []) "general" orig_option) of
+        J.Error err -> ioError $ userError err
+        J.Success val -> return val
+    let general :: Map T.Text J.Value
+        general =
+          case solveTimeLimit opt of
+            Nothing -> orig_general
+            Just t -> Map.insert "time_max" (J.toJSON t) orig_general
+        option :: Map T.Text J.Value
+        option = Map.insert "general" (J.toJSON general) orig_option
+
+    withSystemTempDirectory "printemps" $ \path ->
+      case MPSFile.render def{ optMPSWriteObjName = False } prob' of
+        Left err -> ioError $ userError err
+        Right s -> do
+          let problem_file = path </> "input.mps"
+          TLIO.writeFile problem_file s
+          let option_file = path </> "option.json"
+          J.encodeFile option_file option
+
+          let args = ["-p", option_file] ++ args' ++ [problem_file]
+              onGetLine s = solveLogger opt s
+              onGetErrorLine = solveErrorLogger opt
+          exitcode <- runProcessWithOutputCallback (printempsPath solver) args (Just path) "" onGetLine onGetErrorLine
+          if exitcode /= ExitSuccess then do
+            return $ def{ solStatus = StatusUnknown }
+          else do
+            PrintempsSol.readFile (path </> "incumbent.json")
+
+removeOptionArgs :: [String] -> (Maybe FilePath, [String])
+removeOptionArgs = f Nothing []
+  where
+    f optionFile args [] = (optionFile, reverse args)
+    f _ args ("-p" : fname : xs) = f (Just fname) args xs
+    f optionFile args (x : xs) = f optionFile (x : args) xs
diff --git a/src/Numeric/Optimization/MIP/Solver/SCIP.hs b/src/Numeric/Optimization/MIP/Solver/SCIP.hs
--- a/src/Numeric/Optimization/MIP/Solver/SCIP.hs
+++ b/src/Numeric/Optimization/MIP/Solver/SCIP.hs
@@ -23,23 +23,30 @@
 import System.IO
 import System.IO.Temp
 import qualified Numeric.Optimization.MIP.LPFile as LPFile
+import qualified Numeric.Optimization.MIP.Base as MIP
 import Numeric.Optimization.MIP.Solver.Base
 import qualified Numeric.Optimization.MIP.Solution.SCIP as ScipSol
 import Numeric.Optimization.MIP.Internal.ProcessUtil (runProcessWithOutputCallback)
 
+-- | A solver instance for calling @scip@ command from [SCIP Optimization Suite](https://www.scipopt.org/).
+--
+-- Use 'scip' and record update syntax to modify its field.
 data SCIP
   = SCIP
   { scipPath :: String
+  , scipArgs :: [String]
+  , scipCommands :: [String]
   }
 
 instance Default SCIP where
   def = scip
 
+-- | Default value of t'SCIP'
 scip :: SCIP
-scip = SCIP "scip"
+scip = SCIP "scip" [] []
 
 instance IsSolver SCIP IO where
-  solve solver opt prob = do
+  solve' solver opt prob = do
     case LPFile.render def prob of
       Left err -> ioError $ userError err
       Right lp -> do
@@ -48,17 +55,29 @@
           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"
+            let commands = 
+                  [ "read " ++ show fname1 ] ++
+                  (case solveTimeLimit opt of
+                     Nothing -> []
+                     Just sec -> ["set limits time " ++ show sec]) ++
+                  (case solveTol opt of
+                     Nothing -> []
+                     Just tol ->
+                       [ "set numeric feastol " ++ show (MIP.feasibilityTol tol)
+                       , "set numeric dualfeastol " ++ show (MIP.optimalityTol tol)
+                       , -- https://listserv.zib.de/pipermail/scip/2010-September/000552.html
+                         "set numeric epsilon " ++ show (MIP.integralityTol tol)
                        ]
+                  ) ++
+                  scipCommands solver ++
+                  [ "optimize"
+                  , "write solution " ++ show fname2
+                  , "quit"
+                  ]
+                args = scipArgs solver ++ concat [["-c", cmd] | cmd <- commands]
                 onGetLine = solveLogger opt
                 onGetErrorLine = solveErrorLogger opt
-            exitcode <- runProcessWithOutputCallback (scipPath solver) args "" onGetLine onGetErrorLine
+            exitcode <- runProcessWithOutputCallback (scipPath solver) args Nothing "" onGetLine onGetErrorLine
             case exitcode of
               ExitFailure n -> ioError $ userError $ "exit with " ++ show n
               ExitSuccess -> ScipSol.readFile fname2
diff --git a/test/IsClose.hs b/test/IsClose.hs
new file mode 100644
--- /dev/null
+++ b/test/IsClose.hs
@@ -0,0 +1,142 @@
+{-# OPTIONS_GHC -Wall #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+module IsClose
+  (
+  -- Tolerance type
+    Tol (..)
+
+  -- AllClose class
+  , AllClose (..)
+  , allCloseRawUnit
+  , allCloseRawRealFrac
+  , allCloseRawRealFloat
+
+  -- * Re-exports
+  , Default (..)
+
+  -- * HUnit
+  , assertAllClose
+  ) where
+
+import Data.Default.Class
+import Data.List.NonEmpty (NonEmpty (..))
+import Data.Map (Map)
+import qualified Data.Map as Map
+import Data.Monoid
+import Data.Semigroup
+import qualified Numeric.Optimization.MIP as MIP
+import GHC.Stack (HasCallStack)
+import Test.HUnit
+import Text.Printf
+
+-- ------------------------------------------------------------------------
+
+-- | Tolerance
+--
+-- Values @a@ and @b@ are considered /close/ if @abs (a - b) <= atol + rtol * abs b@.
+data Tol a
+  = Tol
+  { rtol :: a -- ^ The relative tolerance parameter (default: @1e-05@)
+  , atol :: a -- ^ The absolute tolerance parameter (default: @1e-08@)
+  , equalNan :: Bool -- ^ Whether to compare NaN’s as equal (default: @False@)
+  } deriving (Show)
+
+instance RealFrac a => Default (Tol a) where
+  def = Tol
+    { rtol = 1e-05
+    , atol = 1e-08
+    , equalNan = False
+    }
+
+-- ------------------------------------------------------------------------
+
+class Real r => AllClose r a where
+  -- | Returns number of mismatches, number of elements, maximal absolute difference, and maximal relative difference.
+  -- Returns @'Ap' 'Nothing'@ if given values are incomparable.
+  allCloseRaw :: Tol r -> a -> a -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)
+
+  -- | Returns 'True' if the two arrays are equal within the given tolerance; 'False' otherwise.
+  allClose :: Tol r -> a -> a -> Bool
+  allClose tol x y =
+    case getAp (allCloseRaw tol x y) of
+      Nothing -> False
+      Just (Sum numMismatched, _, _, _) -> numMismatched == 0
+
+allCloseRawRealFrac :: RealFrac r => Tol r -> r -> r -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)
+allCloseRawRealFrac t a b = Ap $ Just $
+  ( Sum $ if abs (a - b) <= atol t + rtol t * abs b then 0 else 1
+  , Sum 1
+  , Max (abs (a - b))
+  , Max (abs (a - b) / abs b)
+  )
+
+allCloseRawRealFloat :: RealFloat r => Tol r -> r -> r -> Ap Maybe (Sum Int, Sum Int, Max r, Max r)
+allCloseRawRealFloat t a b
+  | isNaN a /= isNaN b = Ap Nothing
+  | otherwise = Ap $ Just $
+      ( Sum $ if (equalNan t && isNaN a && isNaN b) || a == b || abs (a - b) <= atol t + rtol t * abs b then 0 else 1
+      , Sum 1
+      , Max (abs (a - b))
+      , Max (abs (a - b) / abs b)
+      )
+
+allCloseRawUnit :: Num r => Ap Maybe (Sum Int, Sum Int, Max r, Max r)
+allCloseRawUnit = Ap (Just (Sum 0, Sum 0, Max 0, Max 0))
+
+instance AllClose Rational Rational where
+  allCloseRaw = allCloseRawRealFrac
+
+instance AllClose Double Double where
+  allCloseRaw = allCloseRawRealFloat
+
+instance (AllClose r a) => AllClose r (Maybe a) where
+  allCloseRaw tol (Just a) (Just b) = allCloseRaw tol a b
+  allCloseRaw _ Nothing Nothing = allCloseRawUnit
+  allCloseRaw _ _ _ = Ap Nothing
+
+instance (AllClose r v) => AllClose r [v] where
+  allCloseRaw tol xs ys
+    | length xs == length ys = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip xs ys])
+    | otherwise = Ap Nothing
+
+instance (Ord k, AllClose r v) => AllClose r (Map k v) where
+  allCloseRaw tol m1 m2
+    | Map.keys m1 == Map.keys m2 = sconcat (allCloseRawUnit :| [allCloseRaw tol a b | (a,b) <- zip (Map.elems m1) (Map.elems m2)])
+    | otherwise = Ap Nothing
+
+instance (Real r, AllClose r r) => AllClose r (MIP.Solution r) where
+  allCloseRaw tol a b
+    | MIP.solStatus a == MIP.solStatus b =
+        allCloseRaw tol (MIP.solObjectiveValue a) (MIP.solObjectiveValue b) <>
+        allCloseRaw tol (MIP.solVariables a) (MIP.solVariables b)
+    | otherwise = Ap Nothing
+
+-- ------------------------------------------------------------------------
+
+-- | Assert that two objects are equal up to desired tolerance.
+assertAllClose
+  :: (HasCallStack, AllClose r a, Show r, Show a)
+  => Tol r
+  -> a -- ^ actual
+  -> a -- ^ desired
+  -> Assertion
+assertAllClose tol a b =
+  case getAp (allCloseRaw tol a b) of
+    Nothing ->
+      assertString $ unlines $ header ++ ["x and y nan location mismatch:"] ++ footer
+    Just (Sum numMismatch, Sum numTotal, Max absDiff, Max relDiff)
+      | numMismatch == 0 -> return ()
+      | otherwise ->
+          assertString $ unlines $
+            header ++
+            [ printf "Mismatched elements: %d / %d (%f%%)" numMismatch numTotal (fromIntegral numMismatch * 100 / fromIntegral numTotal :: Double)
+            , " Max absolute difference: " ++ show absDiff
+            , " Max relative difference: " ++ show relDiff
+            ] ++ footer
+   where
+     header, footer :: [String]
+     header = [printf "Not equal to tolerance rtol=%s, atol=%s" (show (rtol tol)) (show (atol tol)), ""]
+     footer = [" x: " ++ show a, " y: " ++ show b]
+
+-- ------------------------------------------------------------------------
diff --git a/test/Test/LPFile.hs b/test/Test/LPFile.hs
--- a/test/Test/LPFile.hs
+++ b/test/Test/LPFile.hs
@@ -12,10 +12,19 @@
 import Numeric.Optimization.MIP.LPFile
 
 case_testdata       = checkString "testdata" testdata
+case_test_bounds_free = checkFile "samples/lp/test-bounds-free.lp"
+case_test_bounds_fixed = checkFile "samples/lp/test-bounds-fixed.lp"
 case_test_indicator = checkFile "samples/lp/test-indicator.lp"
+case_test_problem_name = checkFile "samples/lp/test-problem-name.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_test_semicont  = checkFile "samples/lp/test-semicont.lp"
+case_test_semiint   = checkFile "samples/lp/test-semiint.lp"
+case_test_sos       = checkFile "samples/lp/test-sos.lp"
+case_test_sos2      = checkFile "samples/lp/test-sos2.lp"
+case_test_lazy_constraints = checkFile "samples/lp/test-lazy-constraints.lp"
+case_test_user_cuts = checkFile "samples/lp/test-user-cuts.lp"
 case_empty_obj_1    = checkFile "samples/lp/empty_obj_1.lp"
 case_empty_obj_2    = checkFile "samples/lp/empty_obj_2.lp"
 
@@ -46,7 +55,10 @@
   lp <- parseFile def fname
   case render def lp of
     Left err -> assertFailure ("render failure: " ++ err)
-    Right _ -> return ()
+    Right str -> -- parseString def fname str @?= Right lp
+      case parseString def fname str of
+        Left err -> print str >> undefined
+        Right _ -> return ()
 
 checkString :: String -> String -> Assertion
 checkString name str = do
@@ -55,7 +67,7 @@
     Right lp ->
       case render def lp of
         Left err -> assertFailure ("render failure: " ++ err)
-        Right _ -> return ()
+        Right str -> parseString def name str @?= Right lp
 
 ------------------------------------------------------------------------
 -- Test harness
diff --git a/test/Test/MIP.hs b/test/Test/MIP.hs
--- a/test/Test/MIP.hs
+++ b/test/Test/MIP.hs
@@ -11,6 +11,10 @@
 #endif
 import Data.Maybe
 import qualified Data.Map as Map
+import System.FilePath ((</>))
+import System.IO (utf8)
+import System.IO.Temp
+import Test.QuickCheck.Instances.Text ()
 import Test.Tasty
 import Test.Tasty.HUnit
 import Test.Tasty.QuickCheck
@@ -21,8 +25,23 @@
 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.HiGHS as HiGHSSol
+import qualified Numeric.Optimization.MIP.Solution.Printemps as PrintempsSol
 import qualified Numeric.Optimization.MIP.Solution.SCIP as SCIPSol
 
+case_var_show :: Assertion
+case_var_show = show (MIP.Var "x") @?= show ("x" :: String)
+
+prop_var_name :: Property
+prop_var_name =
+  forAll arbitrary $ \x -> do
+    MIP.varName (MIP.Var x) === x
+
+prop_var_compare :: Property
+prop_var_compare =
+  forAll arbitrary $ \(x1, x2) -> do
+    compare (MIP.Var x1) (MIP.Var x2) === compare x1 x2
+
 prop_status_refl :: Property
 prop_status_refl = forAll arbitrary $ \(x :: MIP.Status) -> do
   x `leq` x
@@ -64,6 +83,143 @@
 instance Arbitrary MIP.Status where
   arbitrary = arbitraryBoundedEnum
 
+case_eval_expr :: Assertion
+case_eval_expr = do
+  MIP.eval MIP.def sol (MIP.varExpr "x" + 2 * MIP.varExpr "y" :: MIP.Expr Rational) @?= 8
+  where
+    sol :: Map.Map MIP.Var Rational
+    sol = Map.fromList [("x", 2), ("y", 3)]
+
+case_eval_constraint :: Assertion
+case_eval_constraint = do
+  let constr1 = MIP.varExpr "x" MIP..<=. MIP.constExpr (0 :: Double)
+  MIP.eval MIP.def (Map.singleton "x" (-1 :: Double)) constr1 @?= True
+  MIP.eval MIP.def (Map.singleton "x" (0 :: Double)) constr1 @?= True
+  MIP.eval MIP.def (Map.singleton "x" (1e-10 :: Double)) constr1 @?= True
+  MIP.eval MIP.def (Map.singleton "x" (1 :: Double)) constr1 @?= False
+
+  let constr2 :: MIP.Constraint Double
+      constr2 = MIP.def{ MIP.constrExpr = MIP.varExpr "x", MIP.constrLB = MIP.NegInf, MIP.constrUB = MIP.PosInf }
+  MIP.eval MIP.def (Map.singleton "x" (0 :: Double)) constr2 @?= True
+  MIP.eval MIP.def (Map.singleton "x" (1 :: Double)) constr2 @?= True
+
+  let constr3 :: MIP.Constraint Double
+      constr3 = MIP.def{ MIP.constrIndicator = Just ("y", 0), MIP.constrExpr = MIP.varExpr "x", MIP.constrLB = MIP.Finite 0, MIP.constrUB = MIP.Finite 1 }
+  MIP.eval MIP.def (Map.fromList [("x", 0 :: Double), ("y", 0)]) constr3 @?= True
+  MIP.eval MIP.def (Map.fromList [("x", 2 :: Double), ("y", 0)]) constr3 @?= False
+  MIP.eval MIP.def (Map.fromList [("x", 2 :: Double), ("y", 1)]) constr3 @?= True
+
+case_eval_sos_constraint :: Assertion
+case_eval_sos_constraint = do
+  let constr1 :: MIP.SOSConstraint Double
+      constr1 =
+        MIP.SOSConstraint
+        { MIP.sosLabel = Nothing
+        , MIP.sosType = MIP.S1
+        , MIP.sosBody = [("x1", 1), ("x2", 2), ("x3", 3)]
+        }
+  MIP.eval MIP.def (Map.fromList [("x1", 0 :: Double), ("x2", 0), ("x3", 0)]) constr1 @?= True
+  MIP.eval MIP.def (Map.fromList [("x1", 1 :: Double), ("x2", 0), ("x3", 0)]) constr1 @?= True
+  MIP.eval MIP.def (Map.fromList [("x1", 1 :: Double), ("x2", 1), ("x3", 0)]) constr1 @?= False
+
+  let constr2 = constr1{ MIP.sosType = MIP.S2 }
+  MIP.eval MIP.def (Map.fromList [("x1", 0 :: Double), ("x2", 0), ("x3", 0)]) constr2 @?= True
+  MIP.eval MIP.def (Map.fromList [("x1", 1 :: Double), ("x2", 0), ("x3", 0)]) constr2 @?= True
+  MIP.eval MIP.def (Map.fromList [("x1", 1 :: Double), ("x2", 1), ("x3", 0)]) constr2 @?= True
+  MIP.eval MIP.def (Map.fromList [("x1", 0 :: Double), ("x2", 1), ("x3", 1)]) constr2 @?= True
+  MIP.eval MIP.def (Map.fromList [("x1", 1 :: Double), ("x2", 0), ("x3", 1)]) constr2 @?= False
+  MIP.eval MIP.def (Map.fromList [("x1", 1 :: Double), ("x2", 1), ("x3", 1)]) constr2 @?= False
+
+case_eval_continuous_variable :: Assertion
+case_eval_continuous_variable = do
+  MIP.eval MIP.def (Map.singleton "x" (0 :: Double)) prob @?= Nothing
+  MIP.eval MIP.def (Map.singleton "x" (1 - 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (1.5 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (2 + 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (3 :: Double)) prob @?= Nothing
+  where
+    prob :: MIP.Problem Double
+    prob =
+      MIP.def
+      { MIP.varDomains = Map.fromList [("x", (MIP.ContinuousVariable, (MIP.Finite 1, MIP.Finite 2)))]
+      }
+
+case_eval_semi_continuous_variable :: Assertion
+case_eval_semi_continuous_variable = do
+  MIP.eval MIP.def (Map.singleton "x" (0 - 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (0 + 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (0.5 :: Double)) prob @?= Nothing
+  MIP.eval MIP.def (Map.singleton "x" (1 - 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (1.5 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (2 + 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (3 :: Double)) prob @?= Nothing
+  where
+    prob :: MIP.Problem Double
+    prob =
+      MIP.def
+      { MIP.varDomains = Map.fromList [("x", (MIP.SemiContinuousVariable, (MIP.Finite 1, MIP.Finite 2)))]
+      }
+
+case_eval_integer_variable :: Assertion
+case_eval_integer_variable = do
+  MIP.eval MIP.def (Map.singleton "x" (0 - 1e-10 :: Double)) prob @?= Nothing
+  MIP.eval MIP.def (Map.singleton "x" (0 + 1e-10 :: Double)) prob @?= Nothing
+  MIP.eval MIP.def (Map.singleton "x" (0.5 :: Double)) prob @?= Nothing
+  MIP.eval MIP.def (Map.singleton "x" (1 - 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (1.5 :: Double)) prob @?= Nothing
+  MIP.eval MIP.def (Map.singleton "x" (2 + 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (3 :: Double)) prob @?= Nothing
+  where
+    prob :: MIP.Problem Double
+    prob =
+      MIP.def
+      { MIP.varDomains = Map.fromList [("x", (MIP.IntegerVariable, (MIP.Finite 1, MIP.Finite 2)))]
+      }
+
+case_eval_semi_integer_variable :: Assertion
+case_eval_semi_integer_variable = do
+  MIP.eval MIP.def (Map.singleton "x" (0 - 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (0 + 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (0.5 :: Double)) prob @?= Nothing
+  MIP.eval MIP.def (Map.singleton "x" (1 - 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (1.5 :: Double)) prob @?= Nothing
+  MIP.eval MIP.def (Map.singleton "x" (2 + 1e-10 :: Double)) prob @?= Just 0
+  MIP.eval MIP.def (Map.singleton "x" (3 :: Double)) prob @?= Nothing
+  where
+    prob :: MIP.Problem Double
+    prob =
+      MIP.def
+      { MIP.varDomains = Map.fromList [("x", (MIP.SemiIntegerVariable, (MIP.Finite 1, MIP.Finite 2)))]
+      }
+
+case_file_io_lp :: Assertion
+case_file_io_lp = do
+  let opt = MIP.def{ MIP.optFileEncoding = Just utf8 }
+  prob <- MIP.readFile opt "samples/lp/test.lp"
+  withSystemTempDirectory "MIP" $ \dir -> do
+    MIP.writeFile opt (dir </> "test.lp") prob
+    prob2 <- MIP.readFile opt (dir </> "test.lp")
+    prob2 @?= prob
+#ifdef WITH_ZLIB
+    MIP.writeFile opt (dir </> "test.lp.gz") prob
+    prob3 <- MIP.readFile opt (dir </> "test.lp.gz")
+    prob3 @?= prob
+#endif
+
+case_file_io_mps :: Assertion
+case_file_io_mps = do
+  let opt = MIP.def{ MIP.optFileEncoding = Just utf8 }
+  prob <- MIP.readFile opt "samples/lp/test.lp"
+  withSystemTempDirectory "MIP" $ \dir -> do
+    MIP.writeFile opt (dir </> "test.mps") prob
+    prob2 <- MIP.readFile opt (dir </> "test.mps")
+    prob2 @?= prob
+#ifdef WITH_ZLIB
+    MIP.writeFile opt (dir </> "test.mps.gz") prob
+    prob3 <- MIP.readFile opt (dir </> "test.mps.gz")
+    prob3 @?= prob
+#endif
+
 case_CBCSol :: Assertion
 case_CBCSol = do
   sol <- CBCSol.readFile "samples/lp/test-solution-cbc.txt"
@@ -139,6 +295,56 @@
   sol <- GurobiSol.readFile "samples/lp/test-solution-gurobi.sol"
   isJust (GurobiSol.solObjectiveValue sol) @?= True
   GurobiSol.parse (GurobiSol.render sol) @?= sol
+
+case_HiGHSSol :: Assertion
+case_HiGHSSol = do
+  sol <- HiGHSSol.readFile "samples/lp/test-solution-highs.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_HiGHSSol_infeasible :: Assertion
+case_HiGHSSol_infeasible = do
+  sol <- HiGHSSol.readFile "samples/lp/test-solution-highs-infeasible.sol"
+  sol @?=
+    MIP.Solution
+    { MIP.solStatus = MIP.StatusInfeasible
+    , MIP.solObjectiveValue = Nothing
+    , MIP.solVariables = Map.empty
+    }
+
+case_HiGHSSol_unbounded :: Assertion
+case_HiGHSSol_unbounded = do
+  sol <- HiGHSSol.readFile "samples/lp/test-solution-highs-unbounded.sol"
+  sol @?=
+    MIP.Solution
+    { MIP.solStatus = MIP.StatusInfeasibleOrUnbounded
+    , MIP.solObjectiveValue = Nothing
+    , MIP.solVariables = Map.empty
+    }
+
+case_PrintempsSol :: Assertion
+case_PrintempsSol = do
+  sol <- PrintempsSol.readFile "samples/lp/test-solution-printemps.json"
+  sol @?=
+    MIP.Solution
+    { MIP.solStatus = MIP.StatusFeasible
+    , MIP.solObjectiveValue = Just (-1.070000e+02)
+    , MIP.solVariables = Map.fromList [("x1", 29), ("x2", 7), ("x3", 22), ("x4", 2)]
+    }
+
+case_PrintempsSol_infeasible :: Assertion
+case_PrintempsSol_infeasible = do
+  sol <- PrintempsSol.readFile "samples/lp/test-solution-printemps-infeasible.json"
+  sol @?=
+    MIP.Solution
+    { MIP.solStatus = MIP.StatusUnknown
+    , MIP.solObjectiveValue = Just 0
+    , MIP.solVariables = Map.fromList [("x", 0), ("y", 0)]
+    }
 
 case_SCIPSol :: Assertion
 case_SCIPSol = do
diff --git a/test/Test/MIPSolver.hs b/test/Test/MIPSolver.hs
--- a/test/Test/MIPSolver.hs
+++ b/test/Test/MIPSolver.hs
@@ -4,21 +4,24 @@
 {-# LANGUAGE TemplateHaskell #-}
 module Test.MIPSolver (mipSolverTestGroup) where
 
+import Control.Arrow ((***))
 import Control.Monad
 import Data.Default.Class
 import qualified Data.Map as Map
+import qualified Data.Set as Set
 import Test.Tasty
 import Test.Tasty.HUnit
 import qualified Numeric.Optimization.MIP as MIP
 import Numeric.Optimization.MIP.Solver
+import IsClose
 
 -- ------------------------------------------------------------------------
 
 case_cbc :: Assertion
 case_cbc = do
   prob <- MIP.readFile def "samples/lp/test.lp"
-  sol <- solve cbc def prob
-  sol @?=
+  sol <- solve cbc def{ solveTol = Just def } prob
+  assertAllClose (def :: Tol Rational) (fmap toRational sol)
     MIP.Solution
     { MIP.solStatus = MIP.StatusOptimal
     , MIP.solObjectiveValue = Just 122.5
@@ -58,8 +61,8 @@
 case_cplex :: Assertion
 case_cplex = do
   prob <- MIP.readFile def "samples/lp/test.lp"
-  sol <- solve cplex def prob
-  sol @?=
+  sol <- solve cplex def{ solveTol = Just def } prob
+  assertAllClose (def :: Tol Rational) (fmap toRational sol)
     MIP.Solution
     { MIP.solStatus = MIP.StatusOptimal
     , MIP.solObjectiveValue = Just 122.5
@@ -99,8 +102,8 @@
 case_glpsol :: Assertion
 case_glpsol = do
   prob <- MIP.readFile def "samples/lp/test.lp"
-  sol <- solve glpsol def prob
-  sol @?=
+  sol <- solve glpsol def{ solveTol = Just def } prob
+  assertAllClose (def :: Tol Rational) (fmap toRational sol)
     MIP.Solution
     { MIP.solStatus = MIP.StatusOptimal
     , MIP.solObjectiveValue = Just 122.5
@@ -129,12 +132,12 @@
 case_gurobiCl :: Assertion
 case_gurobiCl = do
   prob <- MIP.readFile def "samples/lp/test.lp"
-  sol <- solve gurobiCl def prob
-  sol @?=
+  sol <- solve gurobiCl def{ solveTol = Just def } prob
+  assertAllClose (def :: Tol Rational) (fmap toRational 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)]
+    , MIP.solObjectiveValue = Just 122.5
+    , MIP.solVariables = Map.fromList [("x1", 40), ("x2", 10.5), ("x3", 19.5), ("x4", 3)]
     }
 
 case_gurobiCl_unbounded :: Assertion
@@ -167,11 +170,52 @@
 
 -- ------------------------------------------------------------------------
 
+case_highs :: Assertion
+case_highs = do
+  prob <- MIP.readFile def "samples/lp/test.lp"
+  sol <- solve highs def{ solveTol = Just def } prob
+  assertAllClose (def :: Tol Rational) (fmap toRational 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_highs_unbounded :: Assertion
+case_highs_unbounded = do
+  prob <- MIP.readFile def "samples/lp/unbounded-ip.lp"
+  sol <- solve highs 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_highs_infeasible :: Assertion
+case_highs_infeasible = do
+  prob <- MIP.readFile def "samples/lp/infeasible.lp"
+  sol <- solve highs def prob
+  MIP.solStatus sol @?= MIP.StatusInfeasible
+
+case_highs_infeasible2 :: Assertion
+case_highs_infeasible2 = do
+  prob <- MIP.readFile def "samples/lp/glpk-preprocess-bug.lp"
+  sol <- solve highs 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 @?=
+  sol <- solve lpSolve def{ solveTol = Just def } prob
+  assertAllClose (def :: Tol Rational) (fmap toRational sol)
     MIP.Solution
     { MIP.solStatus = MIP.StatusOptimal
     , MIP.solObjectiveValue = Just 122.5
@@ -208,11 +252,28 @@
 
 -- ------------------------------------------------------------------------
 
+case_printemps :: Assertion
+case_printemps = do
+  prob <- MIP.readFile def "samples/lp/test.lp"
+  sol <- solve printemps def{ solveTol = Just def } prob{ MIP.varDomains = fmap ((const MIP.IntegerVariable) *** id) (MIP.varDomains prob) }
+  MIP.solStatus sol @?= MIP.StatusFeasible
+  let vs = MIP.solVariables sol
+  Map.keysSet vs @?= Set.fromList ["x1", "x2", "x3", "x4"]
+  Just (sum [c * (vs Map.! v) | (c, v) <- [(1, "x1"), (2, "x2"), (3, "x3"), (1, "x4")]]) @?= MIP.solObjectiveValue sol 
+
+case_printemps_with_time_limit :: Assertion
+case_printemps_with_time_limit = do
+  prob <- MIP.readFile def "samples/lp/test.lp"
+  _ <- solve printemps def{ solveTimeLimit = Just 1 } prob{ MIP.varDomains = fmap ((const MIP.IntegerVariable) *** id) (MIP.varDomains prob) }
+  return ()
+
+-- ------------------------------------------------------------------------
+
 case_scip :: Assertion
 case_scip = do
   prob <- MIP.readFile def "samples/lp/test.lp"
-  sol <- solve scip def prob
-  sol @?=
+  sol <- solve scip def{ solveTol = Just def } prob
+  assertAllClose (def :: Tol Rational) (fmap toRational sol)
     MIP.Solution
     { MIP.solStatus = MIP.StatusOptimal
     , MIP.solObjectiveValue = Just 122.5
@@ -282,12 +343,26 @@
   , testCase "gurobiCl infeasible2" case_gurobiCl_infeasible2
   ]
 #endif
+#ifdef TEST_HIGHS
+  ++
+  [ testCase "highs" case_highs
+  , testCase "highs unbounded" case_highs_unbounded
+  , testCase "highs infeasible" case_highs_infeasible
+  , testCase "highs infeasible2" case_highs_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_PRINTEMPS
+  ++
+  [ testCase "printemps" case_printemps
+  , testCase "printemps with time limit" case_printemps_with_time_limit
   ]
 #endif
 #ifdef TEST_SCIP
diff --git a/test/Test/MPSFile.hs b/test/Test/MPSFile.hs
--- a/test/Test/MPSFile.hs
+++ b/test/Test/MPSFile.hs
@@ -3,6 +3,7 @@
 
 import Control.Monad
 import Data.Default.Class
+import Data.Either
 import Data.List
 import Data.Maybe
 import Test.Tasty
@@ -20,7 +21,21 @@
 case_quadobj2 = checkFile "samples/mps/quadobj2.mps"
 case_ranges   = checkFile "samples/mps/ranges.mps"
 case_sos      = checkFile "samples/mps/sos.mps"
+case_sos2     = checkFile "samples/mps/test-sos2.mps"
 case_sc       = checkFile "samples/mps/sc.mps"
+case_semicont = checkFile "samples/mps/test-semicont.mps"
+case_semiint_sc = checkFile "samples/mps/test-semiint-sc.mps"
+case_semiint_si = checkFile "samples/mps/test-semiint-si.mps"
+case_bounds_fixed = checkFile "samples/mps/test-bounds-fixed.mps"
+case_bounds_free = checkFile "samples/mps/test-bounds-free.mps"
+case_bounds_MI = checkFile "samples/mps/test-bounds-mi.mps"
+case_bounds_PL = checkFile "samples/mps/test-bounds-pl.mps"
+case_bounds_BV = checkFile "samples/mps/test-bounds-bv.mps"
+case_implicit_bv = checkFile "samples/mps/test-implicit-bv.mps"
+case_negative_upper_bound = checkFile "samples/mps/test-negative-upper-bound.mps"
+case_negative_upper_qcp = checkFile "samples/mps/test-qcp.mps"
+case_lazy_constraints = checkFile "samples/mps/test-lazy-constraints.mps"
+case_user_cuts = checkFile "samples/mps/test-user-cuts.mps"
 
 ------------------------------------------------------------------------
 -- Sample data
@@ -51,14 +66,19 @@
 
 checkFile :: FilePath -> Assertion
 checkFile fname = do
-  _ <- parseFile def fname
-  return ()
+  lp <- parseFile def fname
+  case render def lp of
+    Left err -> assertFailure ("render failure: " ++ err)
+    Right str -> assertBool ("failed to parse " ++ show str)  $ isRight $ parseString def fname str
 
 checkString :: String -> String -> Assertion
 checkString name str = do
   case parseString def name str of
     Left err -> assertFailure (show err)
-    Right lp -> return ()
+    Right lp ->
+      case render def lp of
+        Left err -> assertFailure ("render failure: " ++ err)
+        Right str -> assertBool ("failed to parse " ++ show str)  $ isRight $ parseString def name str
 
 ------------------------------------------------------------------------
 -- Test harness