glpk backend crashing

[Theorell]

Problem description

I'm trying to calculate the chebyshev center of a polytope using optlang with the backends gurobi and glpk. In an effort to find the problem, I have created a minimal example which first finds the minimum of a 2d square (works with both backends) and then extends the problem to find the center of the same square (fails for glpk but not gurobi).

It is of course possibly that this is a general glpk problem, however, that is hard for me to believe, since it fails already for a very trivial example which glpk should be able to solve. Instead, I believe that the problem lies in how the problem is served to glpk.

Code Sample

from cobra.util.solver import solvers
from scipy import sparse as sp
from optlang.symbolics import Zero
import pandas as pd
specific_interface = "glpk"
interface = solvers[specific_interface]
backend = specific_interface

def make_sparse_constraint_system(m, A, b, x):
    A = sp.csr_matrix(A)
    constr_dict = dict()
    for row_ind, bound in enumerate(b):
        constr_dict[interface.Constraint(Zero, ub=bound)] = {
            x[var_ind]: A[row_ind, var_ind]
            for var_ind in A[row_ind, :].nonzero()[1]
        }
    m.add(constr_dict.keys())
    m.update()
    for constr, terms in constr_dict.items():
        constr.set_linear_coefficients(terms)

def gurobi_solve(obj, A, b):
    m = interface.Model()
    # m.configuration.presolve = True
    r_names = [str(r) for r in range(A.shape[1])]
    x = np.array([interface.Variable(name, lb=-np.inf) for name in r_names])
    m.add(x)
    m.objective = interface.Objective(obj @ x, direction="min")
    make_sparse_constraint_system(m, A, b, x)
    m.update()
    m.optimize()
    return m

def main():

    # min square 2d
    b = np.array([1, 1, 1, 1])
    A_ext = np.array([[1, 0], [0, 1], [-1, 0], [0, -1]])
    obj = np.array([1, 1])
    m = gurobi_solve(obj, A_ext, b)
    values = pd.Series(m.primal_values).values
    print('does not crash')

    # chebyshev center of 2d square
    b = np.array([1, 1, 1, 1])
    A_ext = np.array([[1, 0, 1], [0, 1, 1], [-1, 0, 1], [0, -1, 1]])
    obj = np.array([0, 0, -1])
    m = gurobi_solve(obj, A_ext, b)
    values = pd.Series(m.primal_values).values
    print('crashes before this print if backend is glpk')

if __name__ == "__main__":
    main()

does not crash
Assertion failed: teta >= 0.0
Error detected in file simplex/spxchuzr.c at line 294

Context

``` python 3.8 macOS 10.15.7 pip freeze appdirs==1.4.4 certifi==2020.11.8 cobra==0.20.0 colorama==0.4.4 commonmark==0.9.1 depinfo==1.6.0 diskcache==5.1.0 future==0.18.2 h11==0.11.0 h5py==2.10.0 httpcore==0.12.2 httpx==0.16.1 idna==2.10 importlib-resources==3.3.0 mpmath==1.1.0 numexpr==2.7.1 numpy==1.18.3 optlang==1.4.4 pandas==1.0.1 pydantic==1.7.3 Pygments==2.7.2 python-dateutil==2.8.1 python-libsbml==5.18.0 python-libsbml-experimental==5.18.3 pytz==2020.4 rfc3986==1.4.0 rich==6.2.0 ruamel.yaml==0.16.12 ruamel.yaml.clib==0.2.2 scipy==1.4.1 six==1.15.0 sniffio==1.2.0 swiglpk==4.65.1 sympy==1.7 tables==3.6.1 typing-extensions==3.7.4.3 ```

[KristianJensen]

Hi

Thanks for the report. I think the problem might be that you set the lower bound of variables using lb=-np.inf in line 4 of gurobi_solve. Could you try changing it to lb=None?

[cdiener]

Can confirm that it works well with lb=None. That is the output I get with GLPK:

Scaling...
 A: min|aij| =  1.000e+00  max|aij| =  1.000e+00  ratio =  1.000e+00
Problem data seem to be well scaled
does not crash
crashes before this print if backend is glpk

[Theorell]

Yes, it works for me too, thanks! Maybe it would be helpful with a warning or a check for this, since it crashes without stack trace.