Merge pull request #9 from mckib2/mps-write

MPS/LP writer

README.rst

@@ -23,8 +23,8 @@ There are a few things in this package:
 
 - `glpk()` : the wrappers over the solvers (basically acts like Python-friendly `glpsol`)
 - `mpsread()` : convert an MPS file to some matrices
-
-
+- `mpswrite()` : convert matrices to MPS file
+- `lpwrite()` : convert matrices to CPLEX LP file
 
 .. code-block::
 
@@ -33,10 +33,16 @@ There are a few things in this package:
    res = glpk(
        c, A_ub, b_ub, A_eq, b_eq, bounds, solver, sense, maxit, timeout,
        basis_fac, message_level, disp, simplex_options, ip_options,
-       mip_options, libpath)
+       mip_options)
+
+   from glpk import mpsread, mpswrite, lpwrite
 
    c, A_ub, b_ub, A_eq, b_eq, bounds = mpsread(
-       filename, fmt=GLPK.GLP_MPS_FILE, ret_glp_prob=False, libpath=None)
+       filename, fmt=GLPK.GLP_MPS_FILE, ret_glp_prob=False)
+
+   success = mpswrite(c, A_ub, b_ub, A_eq, b_eq, bounds, sense, filename, fmt)
+
+   success = lpwrite(c, A_ub, b_ub, A_eq, b_eq, bounds, sense, filename)
 
 There's lots of information in the docstrings for these functions, please check there for a complete listing and explanation.
 
@@ -73,7 +79,7 @@ Note that there are several projects that aim for something like this, but which
 Why do we want this?
 --------------------
 
-GLPK has a lot of options that the current scipy solvers lack as well as robust MIP support (only basic in HiGHS).  It is also a standard, well known solver in the optimization community.  The only thing that I want that it lacks on an API level is robust support for column generation.  Easy access to GLPK as a backend to `linprog` would be very welcome (to me at least).
+GLPK has a lot of options that the current scipy solvers lack as well as robust MIP support (only basic in HiGHS).  It is also a standard, well known solver in the optimization community.  The only thing that I want that it lacks on an API level is robust support for column generation.  Easy access to GLPK as a backend to `linprog` would be very welcome (to me at least).  I also find access to a linprog LP description (c, A_ub, etc.) to MPS/LP file format convienent for interacting with other solvers such as `HiGHS <https://github.com/ERGO-Code/HiGHS>`_.
 
 Approach
 --------

glpk/__init__.py

@@ -1,4 +1,4 @@
 
 from ._glpk import glpk
-from ._utils import mpsread, mpswrite
+from ._fileio import mpsread, mpswrite, lpwrite
 from ._glpk_defines import GLPK

glpk/_fileio.py

@@ -0,0 +1,159 @@
+import ctypes
+import pathlib
+import os
+
+import numpy as np
+from scipy.sparse import coo_matrix
+
+from ._glpk_defines import GLPK
+from ._utils import _fill_prob
+
+
+def mpsread(filename, fmt=GLPK.GLP_MPS_FILE, ret_glp_prob=False):
+    '''Read an MPS file.
+
+    Parameters
+    ----------
+    filename : str
+        Path to MPS file.
+    fmt : { GLP_MPS_DECK, GLP_MPS_FILE }
+        Type of MPS file (original or free-format). Default is
+        free-format (``GLP_MPS_FILE``).
+    ret_glp_prob : bool
+        Return the ``glp_prob`` structure. If ``False``, `linprog` style
+        matrices and bounds are returned, i.e., ``A_ub``, ``b_ub``, etc.
+        Default is ``False``.
+    '''
+
+    _lib = GLPK()._lib
+
+    # Populate a glp_prob structure
+    prob = _lib.glp_create_prob()
+    filename = str(pathlib.Path(filename).expanduser().resolve())
+    _lib.glp_read_mps(prob, fmt, None, filename.encode())
+
+    # Return the GLPK object if requested
+    if ret_glp_prob:
+        return prob
+
+    # Otherwise read the matrices
+    n = _lib.glp_get_num_cols(prob)
+    c = np.empty((n,))
+    bounds = []
+    for ii in range(1, n+1):
+        c[ii-1] = _lib.glp_get_obj_coef(prob, ii)
+
+        # get lb and ub; have to do it like this for some reason...
+        tipe = _lib.glp_get_col_type(prob, ii)
+        if tipe == GLPK.GLP_FR:
+            bnd = (-np.inf, np.inf)
+        elif tipe == GLPK.GLP_LO:
+            bnd = (_lib.glp_get_col_lb(prob, ii), np.inf)
+        elif tipe == GLPK.GLP_UP:
+            bnd = (-np.inf, _lib.glp_get_col_ub(prob, ii))
+        elif tipe == GLPK.GLP_DB:
+            bnd = (_lib.glp_get_col_lb(prob, ii), _lib.glp_get_col_ub(prob, ii))
+        else:
+            bnd = (_lib.glp_get_col_lb(prob, ii), _lib.glp_get_col_lb(prob, ii))
+        bounds.append(bnd)
+
+    m = _lib.glp_get_num_rows(prob)
+    nnz = _lib.glp_get_num_nz(prob)
+    ub_cols, ub_rows, ub_vals = [], [], []
+    eq_cols, eq_rows, eq_vals = [], [], []
+    b_ub, b_eq = [], []
+    indarr = ctypes.c_int*(n+1)
+    valarr = ctypes.c_double*(n+1)
+    col_ind = indarr()
+    row_val = valarr()
+    for ii in range(1, m+1):
+        row_type = _lib.glp_get_row_type(prob, ii)
+        l = _lib.glp_get_mat_row(prob, ii, col_ind, row_val)
+        if l == 0:
+            # if there are no elements in the row, skip,
+            # otherwise we would add an extra element to RHS
+            continue
+
+        if row_type == GLPK.GLP_UP:
+            b_ub.append(_lib.glp_get_row_ub(prob, ii))
+            for jj in range(1, l+1):
+                ub_cols.append(col_ind[jj] - 1)
+                ub_rows.append(ii - 1)
+                ub_vals.append(row_val[jj])
+        elif row_type == GLPK.GLP_FX:
+            b_eq.append(_lib.glp_get_row_lb(prob, ii))
+            for jj in range(1, l+1):
+                eq_cols.append(col_ind[jj] - 1)
+                eq_rows.append(ii - 1)
+                eq_vals.append(row_val[jj])
+        elif row_type == GLPK.GLP_LO:
+            b_ub.append(-1*_lib.glp_get_row_lb(prob, ii))
+            for jj in range(1, l+1):
+                ub_cols.append(col_ind[jj] - 1)
+                ub_rows.append(ii - 1)
+                ub_vals.append(-1*row_val[jj])
+        else:
+            raise NotImplementedError()
+
+    if ub_vals:
+        # Converting to csc_matrix gets rid of all-zero rows
+        A_ub = coo_matrix((ub_vals, (ub_rows, ub_cols)), shape=(max(ub_rows)+1, n))
+        A_ub = A_ub.tocsc()
+        A_ub = A_ub[A_ub.getnnz(axis=1) > 0]
+    else:
+        A_ub = None
+        b_ub = None
+
+    if eq_vals:
+        A_eq = coo_matrix((eq_vals, (eq_rows, eq_cols)), shape=(max(eq_rows)+1, n))
+        A_eq = A_eq.tocsc()
+        A_eq = A_eq[A_eq.getnnz(axis=1) > 0]
+    else:
+        A_eq = None
+        b_eq = None
+
+    return(c, A_ub, b_ub, A_eq, b_eq, bounds)
+
+
+def mpswrite(
+        c,
+        A_ub=None,
+        b_ub=None,
+        A_eq=None,
+        b_eq=None,
+        bounds=None,
+        sense=GLPK.GLP_MIN,
+        filename='prob.mps',
+        fmt=GLPK.GLP_MPS_FILE):
+    '''Write an MPS file.'''
+
+    filename = pathlib.Path(filename)
+    prob, _lp = _fill_prob(c, A_ub, b_ub, A_eq, b_eq, bounds, sense, filename.stem)
+
+    # Get the library
+    _lib = GLPK()._lib
+
+    # Call write
+    res = _lib.glp_write_mps(prob, fmt, None, str(filename).encode())
+    return res == 0
+
+
+def lpwrite(c,
+        A_ub=None,
+        b_ub=None,
+        A_eq=None,
+        b_eq=None,
+        bounds=None,
+        sense=GLPK.GLP_MIN,
+        filename='prob.lp'):
+    '''Write a CPLEX LP file.'''
+
+    filename = pathlib.Path(filename)
+    prob, _lp = _fill_prob(c, A_ub, b_ub, A_eq, b_eq, bounds, sense, filename.stem)
+
+    # Get the library
+    _lib = GLPK()._lib
+
+    # Call write
+    res = _lib.glp_write_lp(prob, None, str(filename).encode())
+    return res == 0

glpk/_glpk.py

@@ -6,8 +6,10 @@
 import numpy as np
 from scipy.sparse import coo_matrix
 from scipy.optimize import OptimizeWarning, OptimizeResult
+from scipy.optimize._linprog_util import _LPProblem, _clean_inputs
 
 from ._glpk_defines import GLPK, glp_smcp, glp_iptcp, glp_bfcp, glp_iocp
+from ._utils import _fill_prob
 
 
 def glpk(
@@ -261,112 +263,20 @@ def glpk(
     '''
 
     # Housekeeping
-    c = np.array(c).flatten()
-    if bounds is None:
-        # Default bound is [0, inf)
-        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)))
-    if b_ub is None:
-        # just need something iterable
-        b_ub = []
-    if A_eq is None:
-        A_eq = np.empty((0, len(c)))
-    if b_eq is None:
-        b_eq = []
     if simplex_options is None:
         simplex_options = {}
     if ip_options is None:
         ip_options = {}
     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)
+    # Create and fill the GLPK problem struct
+    prob, lp = _fill_prob(c, A_ub, b_ub, A_eq, b_eq, bounds, sense, 'problem-name')
+    c, A_ub, b_ub, A_eq, b_eq, bounds, _x0 = lp
 
     # Get the library
     _lib = GLPK()._lib
 
-    # Create problem instance
-    prob = _lib.glp_create_prob()
-
-    # Give problem a name
-    _lib.glp_set_prob_name(prob, b'problem-name')
-
-    # Set objective name
-    _lib.glp_set_obj_name(prob, b'obj-name')
-
-    # Set objective sense
-    _lib.glp_set_obj_dir(prob, sense)
-
-    # Set objective coefficients and column bounds
-    first_col = _lib.glp_add_cols(prob, len(c))
-    for ii, (c0, bnd) in enumerate(zip(c, bounds)):
-        _lib.glp_set_obj_coef(prob, ii + first_col, c0)
-        _lib.glp_set_col_name(prob, ii + first_col, b'c%d' % ii) # name is c[idx], idx is 0-based index
-
-        if bnd is not None:
-            _lib.glp_set_col_bnds(prob, ii + first_col, bnd[0], bnd[1], bnd[2])
-        # else: default is GLP_FX with lb=0, ub=0
-
-    # Need to load both matrices at the same time
-    first_row = _lib.glp_add_rows(prob, A.shape[0])
-
-    # prepend an element and make 1-based index
-    # b/c GLPK expects indices starting at 1
-    nnz = A.nnz
-    rows = np.concatenate(([-1], A.row + first_row)).astype(ctypes.c_int)
-    cols = np.concatenate(([-1], A.col + first_col)).astype(ctypes.c_int)
-    values = np.concatenate(([0], A.data)).astype(ctypes.c_double)
-    _lib.glp_load_matrix(
-        prob,
-        nnz,
-        rows,
-        cols,
-        values,
-    )
-
-    # Set row bounds
-    # Upper bounds (b_ub):
-    for ii, b0 in enumerate(b_ub):
-        # lb is ignored for upper bounds
-        _lib.glp_set_row_bnds(prob, ii + first_row, GLPK.GLP_UP, 0, b0)
-    # Equalities (b_eq)
-    for ii, b0 in enumerate(b_eq):
-        _lib.glp_set_row_bnds(prob, ii + first_row + len(b_ub), GLPK.GLP_FX, b0, b0)
-
     # Scale the problem
     _lib.glp_scale_prob(prob, GLPK.GLP_SF_AUTO) # do auto scaling for now
 

glpk/_glpk_defines.py

@@ -528,6 +528,13 @@ def __init__(self):
             ctypes.c_char_p,           # filename
         ]
 
+        _lib.glp_write_lp.restype = ctypes.c_int
+        _lib.glp_write_lp.argtypes = [
+            ctypes.POINTER(glp_prob),
+            ctypes.POINTER(glp_mpscp), # should be NULL
+            ctypes.c_char_p,           # filename
+        ]
+
         # LP Basis Factorization
         _lib.glp_get_bfcp.restype = None
         _lib.glp_get_bfcp.argtypes = [ctypes.POINTER(glp_prob), ctypes.POINTER(glp_bfcp)]