Skip to content

Commit

Permalink
fix: catch cplex error for milp (#468)
Browse files Browse the repository at this point in the history
milp don't have reduced costs / shadow prices and cplex throws error
when accessing these attributes. Use ugly catch-all to work-around
  • Loading branch information
hredestig committed Mar 27, 2017
1 parent 5f01ece commit d88e595
Show file tree
Hide file tree
Showing 2 changed files with 55 additions and 42 deletions.
12 changes: 10 additions & 2 deletions cobra/core/solution.py
Original file line number Diff line number Diff line change
Expand Up @@ -291,7 +291,11 @@ def get_solution(model, reactions=None, metabolites=None):
fluxes = zeros(len(reactions))
var_primals = model.solver.primal_values
reduced = zeros(len(reactions))
var_duals = model.solver.reduced_costs
try:
# cplex throws error trying if accessing this attribute and not defined
var_duals = model.solver.reduced_costs
except Exception: # CplexSolverError
var_duals = Series([None] * len(reactions), index=rxn_index)
# reduced costs are not always defined, e.g. for integer problems
if var_duals[rxn_index[0]] is None:
reduced.fill(nan)
Expand All @@ -304,7 +308,11 @@ def get_solution(model, reactions=None, metabolites=None):
fluxes[i] = var_primals[forward] - var_primals[reverse]
reduced[i] = var_duals[forward] - var_duals[reverse]
met_index = [met.id for met in metabolites]
constr_duals = model.solver.shadow_prices
try:
# cplex throws error trying if accessing this attribute and not defined
constr_duals = model.solver.shadow_prices
except Exception: # CplexSolverError
constr_duals = Series([None] * len(metabolites), index=met_index)
shadow = asarray([constr_duals[met.id] for met in metabolites])
return Solution(model.solver.objective.value, model.solver.status,
reactions,
Expand Down
85 changes: 45 additions & 40 deletions cobra/test/test_flux_analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,36 @@
all_solvers = optlang_solvers + list(solver_dict)


def construct_ll_test_model():
test_model = Model()
test_model.add_metabolites(Metabolite("A"))
test_model.add_metabolites(Metabolite("B"))
test_model.add_metabolites(Metabolite("C"))
EX_A = Reaction("EX_A")
EX_A.add_metabolites({test_model.metabolites.A: 1})
DM_C = Reaction("DM_C")
DM_C.add_metabolites({test_model.metabolites.C: -1})
v1 = Reaction("v1")
v1.add_metabolites({test_model.metabolites.A: -1,
test_model.metabolites.B: 1})
v2 = Reaction("v2")
v2.add_metabolites({test_model.metabolites.B: -1,
test_model.metabolites.C: 1})
v3 = Reaction("v3")
v3.add_metabolites({test_model.metabolites.C: -1,
test_model.metabolites.A: 1})
test_model.add_reactions([EX_A, DM_C, v1, v2, v3])
DM_C.objective_coefficient = 1
return test_model


@pytest.fixture(scope="function", params=optlang_solvers)
def ll_test_model(request):
test_model = construct_ll_test_model()
test_model.solver = request.param
return test_model


@contextmanager
def captured_output():
""" A context manager to test the IO summary methods """
Expand Down Expand Up @@ -295,36 +325,13 @@ def test_find_blocked_reactions(self, model, solver):
open_exchanges=True)
assert result == []

@classmethod
def construct_ll_test_model(cls):
test_model = Model()
test_model.add_metabolites(Metabolite("A"))
test_model.add_metabolites(Metabolite("B"))
test_model.add_metabolites(Metabolite("C"))
EX_A = Reaction("EX_A")
EX_A.add_metabolites({test_model.metabolites.A: 1})
DM_C = Reaction("DM_C")
DM_C.add_metabolites({test_model.metabolites.C: -1})
v1 = Reaction("v1")
v1.add_metabolites({test_model.metabolites.A: -1,
test_model.metabolites.B: 1})
v2 = Reaction("v2")
v2.add_metabolites({test_model.metabolites.B: -1,
test_model.metabolites.C: 1})
v3 = Reaction("v3")
v3.add_metabolites({test_model.metabolites.C: -1,
test_model.metabolites.A: 1})
test_model.add_reactions([EX_A, DM_C, v1, v2, v3])
DM_C.objective_coefficient = 1
return test_model

def test_legacy_loopless_benchmark(self, benchmark):
test_model = self.construct_ll_test_model()
test_model = construct_ll_test_model()
benchmark(lambda: construct_loopless_model(test_model).optimize(
solver="cglpk"))

def test_loopless_benchmark_before(self, benchmark):
test_model = self.construct_ll_test_model()
test_model = construct_ll_test_model()

def _():
with test_model:
Expand All @@ -334,15 +341,15 @@ def _():
benchmark(_)

def test_loopless_benchmark_after(self, benchmark):
test_model = self.construct_ll_test_model()
test_model = construct_ll_test_model()
benchmark(loopless_solution, test_model)

def test_legacy_loopless(self):
try:
get_solver_name(mip=True)
except SolverNotFound:
pytest.skip("no MILP solver found")
test_model = self.construct_ll_test_model()
test_model = construct_ll_test_model()
feasible_sol = construct_loopless_model(test_model).optimize(
solver="cglpk")
assert feasible_sol.status == "optimal"
Expand All @@ -354,12 +361,11 @@ def test_legacy_loopless(self):
with pytest.raises(UserWarning):
infeasible_mod.optimize(solver="cglpk")

def test_loopless_solution(self):
test_model = self.construct_ll_test_model()
solution_feasible = loopless_solution(test_model)
test_model.reactions.v3.lower_bound = 1
test_model.optimize()
solution_infeasible = loopless_solution(test_model)
def test_loopless_solution(self, ll_test_model):
solution_feasible = loopless_solution(ll_test_model)
ll_test_model.reactions.v3.lower_bound = 1
ll_test_model.optimize()
solution_infeasible = loopless_solution(ll_test_model)
assert solution_feasible.fluxes["v3"] == 0.0
assert solution_infeasible.fluxes["v3"] == 1.0

Expand All @@ -373,13 +379,12 @@ def test_loopless_solution_fluxes(self, model):
with pytest.raises(UserWarning):
loopless_solution(model, fluxes=fluxes)

def test_add_loopless(self):
test_model = self.construct_ll_test_model()
add_loopless(test_model)
feasible_status = test_model.optimize().status
test_model.reactions.v3.lower_bound = 1
test_model.solver.optimize()
infeasible_status = test_model.solver.status
def test_add_loopless(self, ll_test_model):
add_loopless(ll_test_model)
feasible_status = ll_test_model.optimize().status
ll_test_model.reactions.v3.lower_bound = 1
ll_test_model.solver.optimize()
infeasible_status = ll_test_model.solver.status
assert feasible_status == "optimal"
assert infeasible_status == "infeasible"

Expand Down

0 comments on commit d88e595

Please sign in to comment.