Fix GLPK tolerance issue

This fixes issue #22 by (1) making sure that plot_flux_space uses one solver consistently throughout all optimizations (2) GLPK tol_bnd was reset from 1e-9 to the default (1e-7). Hopefully that change won't have to be reversed, later.

Author: VonAlphaBisZulu

straindesign/glpk_interface.py

@@ -183,7 +183,7 @@ def __init__(self, c, A_ineq, b_ineq, A_eq, b_eq, lb, ub, vtype, indic_constr, M
             glp_add_rows(self.glpk, A_ineq.shape[0] + A_eq.shape[0] + A_indic.shape[0])
             eq_type = [GLP_UP] * len(b_ineq) + [GLP_FX] * len(b_eq) + [GLP_UP] * len(b_indic)
             for i, t, b in zip(range(len(b_ineq + b_eq + b_indic)), eq_type, b_ineq + b_eq + b_indic):
-                glp_set_row_bnds(self.glpk, i + 1, t, b, b)
+                glp_set_row_bnds(self.glpk, i + 1, t, float(b), float(b))
 
             A = sparse.vstack((A_ineq, A_eq, A_indic), 'coo')
             ia = intArray(A.nnz + 1)
@@ -201,7 +201,7 @@ def __init__(self, c, A_ineq, b_ineq, A_eq, b_eq, lb, ub, vtype, indic_constr, M
         self.lp_params = glp_smcp()
         glp_init_smcp(self.lp_params)
         self.max_tlim = self.lp_params.tm_lim
-        self.lp_params.tol_bnd = 1e-9
+        # self.lp_params.tol_bnd = 1e-9
         self.lp_params.msg_lev = 0
         # MILP parameters
         if self.ismilp:
@@ -368,13 +368,13 @@ def set_ub(self, ub):
         type = [glp_get_col_type(self.glpk, i + 1) for i in setvars]
         for i, l, u, t in zip(setvars, lb, ub, type):
             if t in [GLP_FR, GLP_LO] and isinf(u):
-                glp_set_col_bnds(self.glpk, i + 1, t, l, u)
+                glp_set_col_bnds(self.glpk, i + 1, t, float(l), float(u))
             elif t == GLP_UP and isinf(u):
-                glp_set_col_bnds(self.glpk, i + 1, GLP_FR, l, u)
+                glp_set_col_bnds(self.glpk, i + 1, GLP_FR, float(l), float(u))
             elif t in [GLP_LO, GLP_DB, GLP_FX] and not isinf(u) and l < u:
-                glp_set_col_bnds(self.glpk, i + 1, GLP_DB, l, u)
+                glp_set_col_bnds(self.glpk, i + 1, GLP_DB, float(l), float(u))
             elif t in [GLP_LO, GLP_DB, GLP_FX] and not isinf(u) and l == u:
-                glp_set_col_bnds(self.glpk, i + 1, GLP_FX, l, u)
+                glp_set_col_bnds(self.glpk, i + 1, GLP_FX, float(l), float(u))
 
     def set_time_limit(self, t):
         """Set the computation time limit (in seconds)"""
@@ -413,9 +413,9 @@ def add_ineq_constraints(self, A_ineq, b_ineq):
                 val[i + 1] = float(v)
             glp_set_mat_row(self.glpk, numrows + j + 1, numvars, col, val)
             if isinf(b_ineq[j]):
-                glp_set_row_bnds(self.glpk, numrows + j + 1, GLP_FR, -inf, b_ineq[j])
+                glp_set_row_bnds(self.glpk, numrows + j + 1, GLP_FR, -inf, float(b_ineq[j]))
             else:
-                glp_set_row_bnds(self.glpk, numrows + j + 1, GLP_UP, -inf, b_ineq[j])
+                glp_set_row_bnds(self.glpk, numrows + j + 1, GLP_UP, -inf, float(b_ineq[j]))
 
     def add_eq_constraints(self, A_eq, b_eq):
         """Add equality constraints to the model
@@ -442,7 +442,7 @@ def add_eq_constraints(self, A_eq, b_eq):
                 col[i + 1] = i + 1
                 val[i + 1] = float(v)
             glp_set_mat_row(self.glpk, numrows + j + 1, numvars, col, val)
-            glp_set_row_bnds(self.glpk, numrows + j + 1, GLP_FX, b_eq[j], b_eq[j])
+            glp_set_row_bnds(self.glpk, numrows + j + 1, GLP_FX, float(b_eq[j]), float(b_eq[j]))
 
     def set_ineq_constraint(self, idx, a_ineq, b_ineq):
         """Replace a specific inequality constraint
@@ -467,9 +467,9 @@ def set_ineq_constraint(self, idx, a_ineq, b_ineq):
             val[i + 1] = float(v)
         glp_set_mat_row(self.glpk, idx + 1, numvars, col, val)
         if isinf(b_ineq):
-            glp_set_row_bnds(self.glpk, idx + 1, GLP_FR, -inf, b_ineq)
+            glp_set_row_bnds(self.glpk, idx + 1, GLP_FR, -inf, float(b_ineq))
         else:
-            glp_set_row_bnds(self.glpk, idx + 1, GLP_UP, -inf, b_ineq)
+            glp_set_row_bnds(self.glpk, idx + 1, GLP_UP, -inf, float(b_ineq))
 
     def getSolution(self, status) -> list:
         """Retrieve solution from GLPK backend"""

straindesign/lptools.py

@@ -855,11 +855,11 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]:
         elif ax_type[0] == 'yield':
             constr += [[{**axes[0][0], **{k: -v * x for k, v in axes[0][1].items()}}, '=', 0]]
         if ax_type[1] == 'rate':
-            sol_vmin = fba(model, constraints=constr, obj=axes[1][0], obj_sense='minimize')
-            sol_vmax = fba(model, constraints=constr, obj=axes[1][0], obj_sense='maximize')
+            sol_vmin = fba(model, constraints=constr, obj=axes[1][0], solver=solver, obj_sense='minimize')
+            sol_vmax = fba(model, constraints=constr, obj=axes[1][0], solver=solver, obj_sense='maximize')
         elif ax_type[1] == 'yield':
-            sol_vmin = yopt(model, constraints=constr, obj_num=axes[1][0], obj_den=axes[1][1], obj_sense='minimize')
-            sol_vmax = yopt(model, constraints=constr, obj_num=axes[1][0], obj_den=axes[1][1], obj_sense='maximize')
+            sol_vmin = yopt(model, constraints=constr, obj_num=axes[1][0], obj_den=axes[1][1], solver=solver, obj_sense='minimize')
+            sol_vmax = yopt(model, constraints=constr, obj_num=axes[1][0], obj_den=axes[1][1], solver=solver, obj_sense='maximize')
         lb[i] = ceil_dec(sol_vmin.objective_value, 9)
         ub[i] = floor_dec(sol_vmax.objective_value, 9)