From 215c7ba9ee1947a60c1f9e5319ab7e9d4d8ed625 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 30 Oct 2026 23:30:31 +0100 Subject: [PATCH 1/2] Check the starting LP of flexibilize_enz_concs and sensitivity_tuning from the solver status after slim_optimize --- CHANGELOG.md | 1 + src/geckopy/kcat_tuning/sensitivity_tuning.py | 8 +++++--- .../limit_proteins/flexibilize_enz_concs.py | 9 +++++---- tests/test_flexibilize_enz_concs.py | 16 +++++++++++++++- tests/test_sensitivity_tuning.py | 17 ++++++++++++++++- 5 files changed, 42 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 41bc62a..7bb4cf2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,7 @@ * With scipy as cobra's default solver, `load_conventional_gem` reads the model on glpk and `make_ec_model` builds the ecModel on glpk, and both convert the result to scipy (`geckopy.solver_support.switch_to_scipy`). Building the yeast-GEM ecModel on the scipy interface took more than 40 minutes in optlang's dense matrix operations without finishing. * `Enzyme.shadow_price`, `get_enzyme_bottlenecks` and `relax_proteomics_greedy` work on the scipy solver, which has no constraint duals in optlang: `geckopy.solver_support.constraint_shadow_prices` reads the HiGHS row marginals of the last solve and uses them only if they are an optimum of the problem as it is now. `relax_proteomics_greedy` computes the shadow prices of all candidates once per iteration and takes the lowest ID among enzymes with equal absolute shadow price. * `flexibilize_enz_concs` and `sensitivity_tuning` decide that the starting LP is infeasible from the solver status, not from the objective value, so they behave the same on every solver. Before, an infeasible start raised on Gurobi (no objective value) but not on GLPK (a finite objective value). `flexibilize_enz_concs` now returns an empty result for it on all solvers, as MATLAB's `flexibilizeEnzConcs` does, instead of raising on Gurobi. `sensitivity_tuning` raises on all solvers and stops with a warning when a later solve is not optimal. +* The starting-LP check of `flexibilize_enz_concs` and `sensitivity_tuning` reads the solver status after `slim_optimize`, not through the `Solution` of `optimize`. Gurobi holds no primal values for an infeasible LP that it presolved, so building the solution raised `GurobiError: Unable to retrieve attribute 'X'` instead of reaching the infeasible-start handling. * `apply_kcat_constraints` treats an empty `update_rxns` as all reactions, as MATLAB's `applyKcatConstraints` does with an empty `updateRxns`; it used to do nothing. The callers that build the list from a filter (`set_kcat_for_reactions`, the evotune scoring and reset) diff --git a/src/geckopy/kcat_tuning/sensitivity_tuning.py b/src/geckopy/kcat_tuning/sensitivity_tuning.py index a0a67f0..3956098 100644 --- a/src/geckopy/kcat_tuning/sensitivity_tuning.py +++ b/src/geckopy/kcat_tuning/sensitivity_tuning.py @@ -227,11 +227,13 @@ def sensitivity_tuning( if bio_rxn_id: model.objective = bio_rxn_id - initial_sol = model.optimize() - if initial_sol.status != "optimal": + # slim_optimize does not read the primal values, which Gurobi does not hold + # for an infeasible LP that it presolved (reading them raises GurobiError). + model.slim_optimize(error_value=float("nan")) + if model.solver.status != "optimal": raise ValueError( "FBA of the input model is not optimal (status " - f"{initial_sol.status!r}). Reduce the protein " + f"{model.solver.status!r}). Reduce the protein " "pool constraint with set_prot_pool_size and/or check the " "exchange constraints." ) diff --git a/src/geckopy/limit_proteins/flexibilize_enz_concs.py b/src/geckopy/limit_proteins/flexibilize_enz_concs.py index c2f9f81..aab451e 100644 --- a/src/geckopy/limit_proteins/flexibilize_enz_concs.py +++ b/src/geckopy/limit_proteins/flexibilize_enz_concs.py @@ -234,16 +234,17 @@ def flexibilize_enz_concs( units.append(members) unit_frequence = np.zeros(len(units), dtype=int) - initial_sol = model.optimize() - if initial_sol.status != "optimal": + # slim_optimize does not read the primal values, which Gurobi does not hold + # for an infeasible LP that it presolved (reading them raises GurobiError). + pred_growth = float(model.slim_optimize(error_value=float("nan"))) + if model.solver.status != "optimal": logger.warning( "flexibilize_enz_concs: FBA of the input model is %s, so there is " "no growth to flexibilize towards; no enzyme is flexibilized. " "Check the enzyme-concentration and exchange constraints.", - initial_sol.status, + model.solver.status, ) return FlexEnzResult() - pred_growth = float(initial_sol.objective_value) flex_break = False pool_old: Optional[float] = None diff --git a/tests/test_flexibilize_enz_concs.py b/tests/test_flexibilize_enz_concs.py index 935bf58..0282e3d 100644 --- a/tests/test_flexibilize_enz_concs.py +++ b/tests/test_flexibilize_enz_concs.py @@ -278,8 +278,22 @@ def test_infeasible_start_returns_empty_result_on_every_solver(tmp_path, solver) model = _build_toy(adapter, enzyme_ub=1.0, measured_conc=1.0) model.solver = solver model.reactions.get_by_id("R_AB").lower_bound = 50.0 - assert model.optimize().status == "infeasible" + model.slim_optimize(error_value=float("nan")) + assert model.solver.status == "infeasible" result = flexibilize_enz_concs(model) assert result.uniprot_ids == [] assert model.reactions.get_by_id("usage_prot_E").upper_bound == 1.0 assert model.reactions.get_by_id("R_AB").lower_bound == 50.0 + + +def test_infeasible_start_returns_empty_result_when_gurobi_presolves_it(tmp_path): + """Gurobi holds no primal values for an infeasible LP it presolved; the start check must not read them.""" + adapter = _adapter(tmp_path, gr_exp=500.0, bio_rxn="biomass") + model = _build_toy(adapter, enzyme_ub=1.0, measured_conc=1.0) + try: + model.solver = "gurobi" + except Exception as exc: + pytest.skip(f"gurobi unusable: {exc}") + model.solver.problem.Params.Presolve = 2 + model.reactions.get_by_id("R_AB").lower_bound = 50.0 + assert flexibilize_enz_concs(model).uniprot_ids == [] diff --git a/tests/test_sensitivity_tuning.py b/tests/test_sensitivity_tuning.py index 5d8da43..6dcabf3 100644 --- a/tests/test_sensitivity_tuning.py +++ b/tests/test_sensitivity_tuning.py @@ -286,6 +286,21 @@ def test_infeasible_start_raises_on_every_solver(tmp_path, solver): model = _build_toy(adapter, initial_kcat=1.0, enzyme_mw=100.0) model.solver = solver model.reactions.get_by_id("biomass").lower_bound = 500.0 - assert model.optimize().status == "infeasible" + model.slim_optimize(error_value=float("nan")) + assert model.solver.status == "infeasible" + with pytest.raises(ValueError, match="not optimal"): + sensitivity_tuning(model) + + +def test_infeasible_start_raises_when_gurobi_presolves_it(tmp_path): + """Gurobi holds no primal values for an infeasible LP it presolved; the start check must not read them.""" + adapter = _adapter(tmp_path, gr_exp=300.0) + model = _build_toy(adapter, initial_kcat=1.0, enzyme_mw=100.0) + try: + model.solver = "gurobi" + except Exception as exc: + pytest.skip(f"gurobi unusable: {exc}") + model.solver.problem.Params.Presolve = 2 + model.reactions.get_by_id("biomass").lower_bound = 500.0 with pytest.raises(ValueError, match="not optimal"): sensitivity_tuning(model) From b2bc4c113f7b013bc9a1f6c058f8cf76ea7fbf1b Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 31 Oct 2026 00:42:28 +0100 Subject: [PATCH 2/2] optimize_model returns the infeasible solution on every solver; the limit-protein, kcat-tuning and analysis solves use it --- CHANGELOG.md | 2 +- src/geckopy/kcat_tuning/evotune/simulate.py | 3 +- src/geckopy/kcat_tuning/sensitivity_tuning.py | 3 +- src/geckopy/kcat_tuning/sigma_fitter.py | 3 +- .../limit_proteins/flexibilize_enz_concs.py | 15 ++--- .../limit_proteins/get_conc_control_coeffs.py | 9 +-- src/geckopy/limit_proteins/relax_enz_concs.py | 3 +- src/geckopy/solver_support.py | 25 ++++++++ src/geckopy/utilities/bottlenecks.py | 4 +- src/geckopy/utilities/ec_fseof.py | 5 +- src/geckopy/utilities/ec_fva.py | 3 +- tests/test_optimize_model.py | 62 +++++++++++++++++++ 12 files changed, 116 insertions(+), 21 deletions(-) create mode 100644 tests/test_optimize_model.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 7bb4cf2..6e50d57 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,7 +18,7 @@ * With scipy as cobra's default solver, `load_conventional_gem` reads the model on glpk and `make_ec_model` builds the ecModel on glpk, and both convert the result to scipy (`geckopy.solver_support.switch_to_scipy`). Building the yeast-GEM ecModel on the scipy interface took more than 40 minutes in optlang's dense matrix operations without finishing. * `Enzyme.shadow_price`, `get_enzyme_bottlenecks` and `relax_proteomics_greedy` work on the scipy solver, which has no constraint duals in optlang: `geckopy.solver_support.constraint_shadow_prices` reads the HiGHS row marginals of the last solve and uses them only if they are an optimum of the problem as it is now. `relax_proteomics_greedy` computes the shadow prices of all candidates once per iteration and takes the lowest ID among enzymes with equal absolute shadow price. * `flexibilize_enz_concs` and `sensitivity_tuning` decide that the starting LP is infeasible from the solver status, not from the objective value, so they behave the same on every solver. Before, an infeasible start raised on Gurobi (no objective value) but not on GLPK (a finite objective value). `flexibilize_enz_concs` now returns an empty result for it on all solvers, as MATLAB's `flexibilizeEnzConcs` does, instead of raising on Gurobi. `sensitivity_tuning` raises on all solvers and stops with a warning when a later solve is not optimal. -* The starting-LP check of `flexibilize_enz_concs` and `sensitivity_tuning` reads the solver status after `slim_optimize`, not through the `Solution` of `optimize`. Gurobi holds no primal values for an infeasible LP that it presolved, so building the solution raised `GurobiError: Unable to retrieve attribute 'X'` instead of reaching the infeasible-start handling. +* Gurobi holds no primal values for an infeasible LP that it presolved and the scipy interface returns `None` for them, so `model.optimize()` raised `GurobiError: Unable to retrieve attribute 'X'` or `TypeError` there, and the checks that an LP is infeasible failed in a few percent of full test runs. The starting-LP check of `flexibilize_enz_concs` and `sensitivity_tuning` reads the solver status after `slim_optimize`, and `geckopy.solver_support.optimize_model` returns the infeasible solution (status `infeasible`, NaN objective and fluxes) in the other solves of the limit-protein, kcat-tuning and analysis functions. * `apply_kcat_constraints` treats an empty `update_rxns` as all reactions, as MATLAB's `applyKcatConstraints` does with an empty `updateRxns`; it used to do nothing. The callers that build the list from a filter (`set_kcat_for_reactions`, the evotune scoring and reset) diff --git a/src/geckopy/kcat_tuning/evotune/simulate.py b/src/geckopy/kcat_tuning/evotune/simulate.py index d9b40bd..03c43d6 100644 --- a/src/geckopy/kcat_tuning/evotune/simulate.py +++ b/src/geckopy/kcat_tuning/evotune/simulate.py @@ -14,6 +14,7 @@ from typing import TYPE_CHECKING, Callable, Optional import numpy as np +from ...solver_support import optimize_model if TYPE_CHECKING: from ...databases.flux_data import FluxData @@ -146,7 +147,7 @@ def simulate_evotune_dataset( if change_protein_biomass is not None: change_protein_biomass(model, ptot) - sol = model.optimize() + sol = optimize_model(model) if ( sol.status != "optimal" or sol.objective_value is None diff --git a/src/geckopy/kcat_tuning/sensitivity_tuning.py b/src/geckopy/kcat_tuning/sensitivity_tuning.py index 3956098..5443119 100644 --- a/src/geckopy/kcat_tuning/sensitivity_tuning.py +++ b/src/geckopy/kcat_tuning/sensitivity_tuning.py @@ -31,6 +31,7 @@ import numpy as np +from ..solver_support import optimize_model from ..ec_model.constants import ( POOL_ID, PROT_PREFIX, @@ -260,7 +261,7 @@ def sensitivity_tuning( iteration = 1 while iteration <= max_iterations: - sol = model.optimize() + sol = optimize_model(model) if sol.status != "optimal": logger.warning( "sensitivity_tuning: the model has no optimal solution with the " diff --git a/src/geckopy/kcat_tuning/sigma_fitter.py b/src/geckopy/kcat_tuning/sigma_fitter.py index d7fe490..44ca07b 100644 --- a/src/geckopy/kcat_tuning/sigma_fitter.py +++ b/src/geckopy/kcat_tuning/sigma_fitter.py @@ -23,6 +23,7 @@ import numpy as np +from ..solver_support import optimize_model from ..ec_model.pipeline.protein_pool import set_prot_pool_size if TYPE_CHECKING: @@ -75,7 +76,7 @@ def _growth_at_sigma( """ with model: set_prot_pool_size(model, p_tot=p_tot, f=f, sigma=float(sigma)) - sol = model.optimize() + sol = optimize_model(model) if sol.status != "optimal" or sol.objective_value is None: return 0.0 value = float(sol.objective_value) diff --git a/src/geckopy/limit_proteins/flexibilize_enz_concs.py b/src/geckopy/limit_proteins/flexibilize_enz_concs.py index aab451e..79fb0b2 100644 --- a/src/geckopy/limit_proteins/flexibilize_enz_concs.py +++ b/src/geckopy/limit_proteins/flexibilize_enz_concs.py @@ -32,6 +32,7 @@ import numpy as np +from ..solver_support import optimize_model from ..ec_model.constants import ( POOL_EXCHANGE_ID, USAGE_PREFIX, @@ -317,7 +318,7 @@ def flexibilize_enz_concs( float(model.ec.concs[measured_idx[i]]) * (1.0 + increase), ) - new_sol = model.optimize() + new_sol = optimize_model(model) pred_growth = _optimal_growth(new_sol) if verbose: logger.info( @@ -357,7 +358,7 @@ def _try_relax_pool( with model: pool_rxn.upper_bound = 1000.0 bio_rxn.upper_bound = exp_growth - sol = model.optimize() + sol = optimize_model(model) candidate_growth = _optimal_growth(sol) if candidate_growth - pred_growth <= _GROWTH_DELTA_THRESHOLD: @@ -373,7 +374,7 @@ def _try_relax_pool( bio_rxn.lower_bound = candidate_growth model.objective = POOL_EXCHANGE_ID model.objective_direction = "min" - sol = model.optimize() + sol = optimize_model(model) new_pool = _optimal_growth(sol) pool_rxn.upper_bound = new_pool @@ -428,7 +429,7 @@ def _build_result( rxn.upper_bound = exp_growth * (1 + band) model.objective = POOL_EXCHANGE_ID model.objective_direction = "min" - sol = model.optimize() + sol = optimize_model(model) new_concs = np.array( [float(sol.fluxes[rid]) for rid in usage_rxn_ids], dtype=float, @@ -493,7 +494,7 @@ def _unit_scores( from .get_conc_control_coeffs import _solution_is_optimal, _solver_supports_duals scores = np.zeros(len(units)) - sol = model.optimize() + sol = optimize_model(model) if not _solution_is_optimal(sol): return scores candidates = [] @@ -519,7 +520,7 @@ def _unit_scores( pool_rxn = model.reactions.get_by_id(POOL_EXCHANGE_ID) with model: pool_rxn.upper_bound = max(pool_rxn.upper_bound, 1000.0) - base = model.optimize() + base = optimize_model(model) if not _solution_is_optimal(base): return scores baseline = float(base.objective_value) @@ -527,7 +528,7 @@ def _unit_scores( with model: for r in reactions: r.upper_bound = r.upper_bound * fold_change - probe = model.optimize() + probe = optimize_model(model) if _solution_is_optimal(probe): gain = float(probe.objective_value) - baseline scores[k] = gain / (fold_change - 1.0) if gain > _GROWTH_DELTA_THRESHOLD else 0.0 diff --git a/src/geckopy/limit_proteins/get_conc_control_coeffs.py b/src/geckopy/limit_proteins/get_conc_control_coeffs.py index 012ac44..280634b 100644 --- a/src/geckopy/limit_proteins/get_conc_control_coeffs.py +++ b/src/geckopy/limit_proteins/get_conc_control_coeffs.py @@ -14,6 +14,7 @@ from ..ec_model.ec_model import EcModel +from ..solver_support import optimize_model from ..ec_model.constants import USAGE_PREFIX logger = logging.getLogger(__name__) @@ -130,7 +131,7 @@ def _shadow_price_coeffs( enz_mask = np.zeros(n, dtype=bool) coeffs = np.zeros(n, dtype=float) - sol = model.optimize() + sol = optimize_model(model) if not _solution_is_optimal(sol): return enz_mask, coeffs initial_growth = sol.objective_value @@ -167,7 +168,7 @@ def _shadow_price_coeffs( predicted_gain = rc * rxn.upper_bound * _VERIFY_STEP with model: rxn.upper_bound = rxn.upper_bound * (1.0 + _VERIFY_STEP) - new_sol = model.optimize() + new_sol = optimize_model(model) if not _solution_is_optimal(new_sol): continue gain = new_sol.objective_value - initial_growth @@ -248,7 +249,7 @@ def _finite_difference_coeffs( enz_mask = np.zeros(n, dtype=bool) coeffs = np.zeros(n, dtype=float) - initial_sol = model.optimize() + initial_sol = optimize_model(model) if not _solution_is_optimal(initial_sol): return enz_mask, coeffs initial_growth = initial_sol.objective_value @@ -270,7 +271,7 @@ def _finite_difference_coeffs( new_ub = prev_ub * fold_change with model: rxn.upper_bound = new_ub - new_sol = model.optimize() + new_sol = optimize_model(model) if not _solution_is_optimal(new_sol): continue new_growth = new_sol.objective_value diff --git a/src/geckopy/limit_proteins/relax_enz_concs.py b/src/geckopy/limit_proteins/relax_enz_concs.py index be3b8ba..a700fbd 100644 --- a/src/geckopy/limit_proteins/relax_enz_concs.py +++ b/src/geckopy/limit_proteins/relax_enz_concs.py @@ -24,6 +24,7 @@ import numpy as np +from ..solver_support import optimize_model from ..ec_model.constants import POOL_EXCHANGE_ID, USAGE_PREFIX from .flexibilize_enz_concs import ( FlexEnzResult, @@ -161,7 +162,7 @@ def relax_enz_concs( terms = [(1.0 if weight == "relative" else sum(max(base[i], floor) for i in unit)) * x for x, unit in folds] terms.append(pool_weight * (1.0 if weight == "relative" else pool_cap) * pool_x) model.objective = model.problem.Objective(add(terms), direction="min") - solution = model.optimize() + solution = optimize_model(model) feasible = solution.status == "optimal" values = model.solver.primal_values if feasible else {} finally: diff --git a/src/geckopy/solver_support.py b/src/geckopy/solver_support.py index 01e4c7e..88d7707 100644 --- a/src/geckopy/solver_support.py +++ b/src/geckopy/solver_support.py @@ -74,8 +74,33 @@ "edit_on_fast_solver", "is_scipy_interface", "constraint_shadow_prices", + "optimize_model", ] +def optimize_model(model, **kwargs): + """``model.optimize()`` that returns the solution of an infeasible LP on every solver. + + Gurobi holds no primal values for an LP that it presolved as infeasible and optlang's + scipy interface returns ``None`` for them, and cobra reads them while it builds the + solution, so ``model.optimize()`` raises ``GurobiError`` or ``TypeError`` there. When + the solver status is not optimal this function returns a ``Solution`` with that status, + a NaN objective and NaN fluxes instead; with an optimal status every error is raised. + """ + try: + return model.optimize(**kwargs) + except Exception: + if model.solver.status == "optimal": + raise + import pandas as pd + from cobra.core import Solution + + return Solution( + objective_value=float("nan"), + status=model.solver.status, + fluxes=pd.Series(float("nan"), index=[r.id for r in model.reactions]), + ) + + # Relative tolerance of the optimality certificate, scaled by the size of the quantity checked. _CERTIFICATE_TOL = 1e-6 diff --git a/src/geckopy/utilities/bottlenecks.py b/src/geckopy/utilities/bottlenecks.py index 39afee3..ffdfdef 100644 --- a/src/geckopy/utilities/bottlenecks.py +++ b/src/geckopy/utilities/bottlenecks.py @@ -22,7 +22,7 @@ import pandas as pd -from ..solver_support import constraint_shadow_prices +from ..solver_support import constraint_shadow_prices, optimize_model if TYPE_CHECKING: from ..ec_model.ec_model import EcModel @@ -59,7 +59,7 @@ def get_enzyme_bottlenecks(model: "EcModel", *, top: int = 10) -> pd.DataFrame: RuntimeError If the solve status is not optimal. """ - sol = model.optimize() + sol = optimize_model(model) if sol.status != "optimal": raise RuntimeError( f"Solver status {sol.status!r}; cannot rank bottlenecks" diff --git a/src/geckopy/utilities/ec_fseof.py b/src/geckopy/utilities/ec_fseof.py index c41b79e..9bebcee 100644 --- a/src/geckopy/utilities/ec_fseof.py +++ b/src/geckopy/utilities/ec_fseof.py @@ -47,6 +47,7 @@ from typing import TYPE_CHECKING, Optional from raven_toolbox.analysis.fseof import FSEOFResult, fseof +from ..solver_support import optimize_model if TYPE_CHECKING: from ..ec_model.ec_model import EcModel @@ -162,7 +163,7 @@ def _check_cs_consistency( cs_rxn_obj = model.reactions.get_by_id(cs_rxn_id) with model: model.objective = bio_rxn_id - sol = model.optimize() + sol = optimize_model(model) if sol.status != "optimal": return cs_flux = float(sol.fluxes.get(cs_rxn_id, 0.0)) @@ -187,7 +188,7 @@ def _biomass_optimal_target_flux( """ with model: model.objective = bio_rxn_id - sol = model.optimize() + sol = optimize_model(model) if sol.status != "optimal": return 0.0 return float(sol.fluxes.get(prod_target_rxn_id, 0.0)) diff --git a/src/geckopy/utilities/ec_fva.py b/src/geckopy/utilities/ec_fva.py index 391bb76..f51e520 100644 --- a/src/geckopy/utilities/ec_fva.py +++ b/src/geckopy/utilities/ec_fva.py @@ -31,6 +31,7 @@ import pandas as pd from cobra.util import ProcessPool +from ..solver_support import optimize_model from ..ec_model.constants import canonicalize_rxn_id try: @@ -258,7 +259,7 @@ def _solve_for_conv_rxn( objective.update({rxns[i]: -1.0 for i in reverse_neg}) ec_model.objective = objective ec_model.objective_direction = sense - sol = ec_model.optimize() + sol = optimize_model(ec_model) if sol.status != "optimal": return out for i, rxn in enumerate(rxns): diff --git a/tests/test_optimize_model.py b/tests/test_optimize_model.py new file mode 100644 index 0000000..b09ed42 --- /dev/null +++ b/tests/test_optimize_model.py @@ -0,0 +1,62 @@ +"""Tests for optimize_model.""" +import math + +import cobra +import pytest + +from geckopy.solver_support import optimize_model + + +def _infeasible_model(solver: str) -> cobra.Model: + model = cobra.Model("m") + a = cobra.Metabolite("A", compartment="c") + rxn = cobra.Reaction("R1") + rxn.add_metabolites({a: 1.0}) + rxn.bounds = (0.0, 10.0) + sink = cobra.Reaction("S1") + sink.add_metabolites({a: -1.0}) + sink.bounds = (20.0, 30.0) + model.add_reactions([rxn, sink]) + model.objective = "R1" + model.solver = solver + return model + + +@pytest.mark.parametrize("solver", ["glpk", "scipy"]) +def test_infeasible_lp_gives_an_infeasible_solution(solver): + solution = optimize_model(_infeasible_model(solver)) + assert solution.status == "infeasible" + + +def test_feasible_lp_matches_model_optimize(): + model = _infeasible_model("glpk") + model.reactions.get_by_id("S1").bounds = (0.0, 5.0) + solution = optimize_model(model) + assert solution.status == "optimal" + assert solution.objective_value == pytest.approx(5.0) + + +def test_gurobi_presolved_infeasible_lp_gives_an_infeasible_solution(): + """Gurobi holds no primal values for an infeasible LP that it presolved, which makes + ``model.optimize()`` raise ``GurobiError``.""" + try: + model = _infeasible_model("gurobi") + except Exception as exc: + pytest.skip(f"gurobi unusable: {exc}") + model.solver.problem.Params.Presolve = 2 + solution = optimize_model(model) + assert solution.status == "infeasible" + assert math.isnan(solution.objective_value) + assert set(solution.fluxes.index) == {"R1", "S1"} + + +def test_other_errors_are_raised(): + class Broken: + class solver: + status = "optimal" + + def optimize(self, **kwargs): + raise RuntimeError("boom") + + with pytest.raises(RuntimeError, match="boom"): + optimize_model(Broken())