Return the infeasible solution from every solve instead of raising

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.
+* 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)

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

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,
@@ -227,11 +228,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."
         )
@@ -258,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 "

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)

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,
@@ -234,16 +235,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
@@ -316,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(
@@ -356,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:
@@ -372,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
@@ -427,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,
@@ -492,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 = []
@@ -518,15 +520,15 @@ 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)
         for k, reactions, _ in candidates:
             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

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

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:

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
 

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"

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))

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):

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 == []