From a9736643a780ec1b2b2e5072115137034486d3c2 Mon Sep 17 00:00:00 2001 From: mckib2 Date: Mon, 17 Aug 2020 16:00:50 -0700 Subject: [PATCH 1/2] Use linprog-style bounds for glpk --- README.rst | 9 +++---- glpk/_glpk.py | 73 +++++++++++++++++++++++++++++++++++++-------------- setup.py | 20 +++++++++----- 3 files changed, 71 insertions(+), 31 deletions(-) diff --git a/README.rst b/README.rst index 8d99e13..d3a2116 100644 --- a/README.rst +++ b/README.rst @@ -12,7 +12,8 @@ Should be an easy pip installation: pip install scikit-glpk -A Python-compatible C compiler is required to build GLPK from source. +A Python-compatible C compiler is required to build GLPK from source. By default a precompiled +wheel is used during compilation, so you don't have to compile if you don't want to. Usage ----- @@ -40,15 +41,13 @@ There's lots of information in the docstrings for these functions, please check Notice that `glpk` is the wrapper and `GLPK` acts as a namespace that holds constants. -`bounds` is also behaves a little differently for both of these: +`bounds` behaves the same as the `scipy.optimize.linprog` `bounds` argument. They are converted to GLPK-style bounds first thing. -- as an input to `glpk`, `bounds` is a list of triplets in the style of GLPK (probably should be converted to `linprog` conventions) -- as an output of `mpsread`, `bounds` is a list of tuples in the style of `linprog` GLPK stuffs ----------- -GLPK is installed with the module and a `linprog`-like wrapper is provided with a ctypes backend. A pared-down version of glpk-4.65 is vendored from `here `_ and compile instructions are scraped from the makefiles. I'll try my best to support cross-platform pip installations. +GLPK is installed with the module and a `linprog`-like wrapper is provided with a ctypes backend. A pared-down version of glpk-4.65 is vendored from `here `_ and compile instructions are scraped from the makefiles. I'll try my best to support cross-platform pip installations. Wheels are current being built for Linux/Mac/Windows. Background diff --git a/glpk/_glpk.py b/glpk/_glpk.py index f17375b..c775016 100644 --- a/glpk/_glpk.py +++ b/glpk/_glpk.py @@ -9,6 +9,7 @@ from ._glpk_defines import GLPK, glp_smcp, glp_iptcp, glp_bfcp, glp_iocp + def glpk( c, A_ub=None, @@ -39,15 +40,13 @@ def glpk( A_eq : 2-D array (k, n) scipy.sparse.coo_matrix b_eq : 1-D array (k,) - bounds : list (n,) of tuple (3,) + bounds : None or list (n,) of tuple (2,) or tuple (2,) The jth entry in the list corresponds to the jth objective coefficient. - Each entry is made up of a tuple describing the bounds: - - - ``type`` : {GLP_FR, GLP_LO, GLP_UP, GLP_DB, GLP_FX} - - ``lb`` : double - - ``up`` : double - - If the entry is ``None``, then ``type=GLP_FX`` and ``ub==lb==0``. + Each entry is made up of a tuple describing the bounds. + Use None to indicate that there is no bound. By default, bounds are + (0, None) (all decision variables are non-negative). If a single tuple + (min, max) is provided, then min and max will serve as bounds for all + decision variables. solver : { 'simplex', 'interior', 'mip' } Use simplex (LP/MIP) or interior point method (LP only). Default is ``simplex``. @@ -137,7 +136,8 @@ def glpk( - ``qmd`` : quotient minimum degree ordering - ``amd`` : approximate minimum degree ordering - ``symamd`` : approximate minimum degree ordering - algorithm for Cholesky factorization of symmetric matrices. + algorithm for Cholesky factorization of symmetric + matrices. mip_options : dict Options specific to MIP solver. @@ -163,8 +163,8 @@ def glpk( - ``last`` : branch on last integer variable - ``mostf`` : branch on most fractional variable - ``drtom`` : branch using heuristic by Driebeck and Tomlin - - ``pcost`` : branch using hybrid pseudocost heuristic (may be - useful for hard instances) + - ``pcost`` : branch using hybrid pseudocost heuristic + (may be useful for hard instances) - backtrack : { 'dfs', 'bfs', 'bestp', 'bestb' } Backtracking rule. Default is ``bestb``. @@ -218,9 +218,12 @@ def glpk( ''' # Housekeeping + c = np.array(c).flatten() if bounds is None: # Default bound is [0, inf) - bounds = [(GLPK.GLP_LO, 0, 0)]*len(c) + bounds = [(0, None)]*len(c) + elif np.array(bounds).size == 2: + bounds = [(bounds[0], bounds[1])]*len(c) if A_ub is None: # we need numpy arrays A_ub = np.empty((0, len(c))) @@ -238,6 +241,38 @@ def glpk( if mip_options is None: mip_options = {} + # Make sure input is good + assert len(A_ub) == len(b_ub), 'A_ub, b_ub must have same number of rows!' + assert len(A_eq) == len(b_eq), 'A_eq, b_eq must have same number of rows!' + if A_ub: + assert len(A_ub[0]) == len(c), 'A_ub must have as many columns as c!' + if A_eq: + assert len(A_eq[0]) == len(c), 'A_eq must have as many columns as c!' + + # coo for (i, j, val) format + A = coo_matrix(np.concatenate((A_ub, A_eq), axis=0)) + + assert len(bounds) == len(c), 'Must have as many bounds as coefficients!' + assert all(len(bnd) == 2 for bnd in bounds) + + # Convert linprog-style bounds to GLPK-style bounds + for ii, (lb, ub) in enumerate(bounds): + if lb in {-np.inf, None} and ub in {np.inf, None}: + # -inf < x < inf + bounds[ii] = (GLPK.GLP_FR, 0, 0) + elif lb in {-np.inf, None}: + # -inf < x <= ub + bounds[ii] = (GLPK.GLP_UP, 0, ub) + elif ub in {np.inf, None}: + # lb <= x < inf + bounds[ii] = (GLPK.GLP_LO, lb, 0) + elif ub < lb: + # lb <= x <= ub + bounds[ii] = (GLPK.GLP_DB, lb, ub) + else: + # lb == x == up + bounds[ii] = (GLPK.GLP_FX, lb, ub) + # Get the library _lib = GLPK()._lib @@ -264,8 +299,6 @@ def glpk( # else: default is GLP_FX with lb=0, ub=0 # Need to load both matrices at the same time - A = coo_matrix(np.concatenate((A_ub, A_eq), axis=0)) # coo for (i, j, val) format - b = np.concatenate((b_ub, b_eq)) first_row = _lib.glp_add_rows(prob, A.shape[0]) # prepend an element and make 1-based index @@ -376,8 +409,8 @@ def glpk( res.x = np.array([_lib.glp_get_col_prim(prob, ii) for ii in range(1, len(c)+1)]) res.dual = np.array([_lib.glp_get_col_dual(prob, ii) for ii in range(1, len(b_ub)+1)]) - # We don't get slack without doing sensitivity analysis since GLPK uses - # auxiliary variables instead of slack! + # We don't get slack without doing sensitivity analysis since GLPK + # uses auxiliary variables instead of slack! res.slack = b_ub - A_ub @ res.x res.con = b_eq - A_eq @ res.x @@ -528,7 +561,6 @@ def glpk( res.fun = _lib.glp_mip_obj_val(prob) res.x = np.array([_lib.glp_mip_col_val(prob, ii) for ii in range(1, len(c)+1)]) - else: raise ValueError('"%s" is not a recognized solver.' % solver) @@ -536,11 +568,12 @@ def glpk( _lib.glp_delete_prob(prob) # Map status codes to scipy: - #res.status = { - # GLPK.GLP_OPT: 0, - #}[res.status] + # res.status = { + # GLPK.GLP_OPT: 0, + # }[res.status] return res + if __name__ == '__main__': pass diff --git a/setup.py b/setup.py index 071be41..840d64f 100644 --- a/setup.py +++ b/setup.py @@ -1,13 +1,15 @@ '''Install scikit-glpk bindings.''' -from setuptools import find_packages +import pathlib from distutils.core import Extension, setup from distutils.command.build_ext import build_ext as _build_ext -import pathlib +from setuptools import find_packages GLPK_SRC_DIR = pathlib.Path('glpk-4.65/src') + def scrape_makefile_list(filename, START_TOKEN, END_TOKEN): + '''Grab tags from GLPK makefile.''' with open(str(filename), 'r', encoding='utf-8') as f: _contents = f.read() sidx = _contents.find(START_TOKEN) @@ -15,25 +17,31 @@ def scrape_makefile_list(filename, START_TOKEN, END_TOKEN): lines = _contents[sidx+len(START_TOKEN):eidx].splitlines() return [str(_l.replace('\\', '').strip()) for _l in lines] + class build_ext(_build_ext): + '''Overide get_export_symbols to provide them for Windows DLL.''' def get_export_symbols(self, ext): '''Only for generating Windows DLL.''' def_file = GLPK_SRC_DIR / '../w64/glpk_4_65.def' return scrape_makefile_list(def_file, 'EXPORTS\n', ';; end of file ;;') + # Get sources for GLPK makefile = GLPK_SRC_DIR / 'Makefile.am' -sources = scrape_makefile_list(makefile, 'libglpk_la_SOURCES = \\\n', '\n## eof ##') +sources = scrape_makefile_list( + makefile, 'libglpk_la_SOURCES = \\\n', '\n## eof ##') sources = [str(GLPK_SRC_DIR / _s) for _s in sources] # Get include dirs for GLPK -include_dirs = scrape_makefile_list(makefile, 'libglpk_la_CPPFLAGS = \\\n', '\nlibglpk_la_LDFLAGS') -include_dirs = [str(GLPK_SRC_DIR / _d[len('-I($srcdir)/'):]) for _d in include_dirs] +include_dirs = scrape_makefile_list( + makefile, 'libglpk_la_CPPFLAGS = \\\n', '\nlibglpk_la_LDFLAGS') +include_dirs = [ + str(GLPK_SRC_DIR / _d[len('-I($srcdir)/'):]) for _d in include_dirs] setup( name='scikit-glpk', - version='0.1.4', + version='0.2.0', author='Nicholas McKibben', author_email='nicholas.bgp@gmail.com', url='https://github.com/mckib2/scikit-glpk', From 8936c3779828cd629a3af8ec34c71b90e1b814bb Mon Sep 17 00:00:00 2001 From: mckib2 Date: Mon, 17 Aug 2020 16:33:13 -0700 Subject: [PATCH 2/2] Add tolerance options for simplex and mip --- README.rst | 5 ----- glpk/_glpk.py | 62 ++++++++++++++++++++++++++++++++++++++++++++++----- setup.py | 2 +- 3 files changed, 58 insertions(+), 11 deletions(-) diff --git a/README.rst b/README.rst index d3a2116..a1e5195 100644 --- a/README.rst +++ b/README.rst @@ -84,8 +84,3 @@ Since the underlying API is quite simple and written in C and only C, `ctypes` i GLPK is packaged but I may want to make it so the user can optionally specify where the installation is on a user's computer (i.e., path to the shared library) so GLPK is not packaged with `scikit-glpk` and/or scipy. `linprog` could then presumably route the problem to the GLPK backend instead of HiGHS or the existing native python solvers. The `ctypes` wrapper is required for integrating GLPK into the Python runtime. Instead of using MPS files to communicate problems and reading solutions from files, `scipy.sparse.coo_matrix` and `numpy` arrays can be passed directly to the library. More information can be extracted from GLPK this way as well (For example, there is no way to get iteration count except by reading directly from the underlying structs. It is only ever printed to stdout, no other way to get it). - -TODO ----- - -- Several GLPK solver options (notably tolerances) not wrapped yet diff --git a/glpk/_glpk.py b/glpk/_glpk.py index c775016..3fe3137 100644 --- a/glpk/_glpk.py +++ b/glpk/_glpk.py @@ -116,13 +116,38 @@ def glpk( - ``norelax`` : standard "textbook" ratio test - ``flip`` : long-step ratio test + - tol_bnd : double + Tolerance used to check if the basic solution is primal + feasible. (Default: 1e-7). + + - tol_dj : double + Tolerance used to check if the basic solution is dual + feasible. (Default: 1e-7). + + - tol_piv : double + Tolerance used to choose eligble pivotal elements of + the simplex table. (Default: 1e-10). + + - obj_ll : double + Lower limit of the objective function. If the objective + function reaches this limit and continues decreasing, + the solver terminates the search. Used in the dual simplex + only. (Default: -DBL_MAX -- the largest finite float64). + + - obj_ul : double + Upper limit of the objective function. If the objective + function reaches this limit and continues increasing, + the solver terminates the search. Used in the dual simplex + only. (Default: +DBL_MAX -- the largest finite float64). + - presolve : bool Use presolver (assumes ``scale=True`` and ``init_basis='adv'``. Default is ``True``. - exact : bool Use simplex method based on exact arithmetic. - Default is ``False``. + Default is ``False``. If ``True``, all other + ``simplex_option`` fields are ignored. ip_options : dict Options specific to interior-pooint solver. @@ -208,13 +233,31 @@ def glpk( - ``clique`` : clique cuts - ``all`` : generate all cuts above - - gap_tol : float - Relative mip gap tolerance. - + - tol_int : float + Absolute tolerance used to check if optimal solution to the + current LP relaxation is integer feasible. + (Default: 1e-5). + - tol_obj : float + Relative tolerance used to check if the objective value in + optimal solution to the current LP relaxation is not better + than in the best known integer feasible solution. + (Default: 1e-7). + - mip_gap : float + Relative mip gap tolerance. If the relative mip gap for + currently known best integer feasiblesolution falls below + this tolerance, the solver terminates the search. This allows + obtaining suboptimal integer feasible solutions if solving the + problem to optimality takes too long time. + (Default: 0.0). - bound : float add inequality obj <= bound (minimization) or obj >= bound (maximization) to integer feasibility problem (assumes ``minisat=True``). + + Notes + ----- + In general, don't change tolerances without a detailed understanding + of their purposes. ''' # Housekeeping @@ -375,6 +418,13 @@ def glpk( 'norelax': GLPK.GLP_RT_STD, 'flip': GLPK.GLP_RT_FLIP, }[simplex_options.get('ratio', 'relax')] + smcp.tol_bnd = simplex_options.get('tol_bnd', 1e-7) + smcp.tol_dj = simplex_options.get('tol_dj', 1e-7) + smcp.tol_piv = simplex_options.get('tol_piv', 1e-10) + if simplex_options.get('obj_ll', False): + smcp.obj_ll = simplex_options['obj_ll'] + if simplex_options.get('obj_ul', False): + smcp.obj_ul = simplex_options['obj_ul'] smcp.it_lim = maxit smcp.tm_lim = timeout smcp.presolve = { @@ -528,7 +578,9 @@ def glpk( if 'clique' in cuts: iocp.clq_cuts = GLPK.GLP_ON - iocp.mip_gap = mip_options.get('gap_tol', 0.0) + iocp.tol_int = mip_options.get('tol_int', 1e-5) + iocp.tol_obj = mip_options.get('tol_obj', 1e-7) + iocp.mip_gap = mip_options.get('mip_gap', 0.0) iocp.tm_lim = timeout iocp.presolve = { True: GLPK.GLP_ON, diff --git a/setup.py b/setup.py index 840d64f..fbf51ae 100644 --- a/setup.py +++ b/setup.py @@ -41,7 +41,7 @@ def get_export_symbols(self, ext): setup( name='scikit-glpk', - version='0.2.0', + version='0.2.1', author='Nicholas McKibben', author_email='nicholas.bgp@gmail.com', url='https://github.com/mckib2/scikit-glpk',