From 17720c9c92e624eee03b43e5ad72aaa214b26b45 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:39:04 +0200 Subject: [PATCH 01/17] fix: read the SBO term from YAML as a string read_yaml_model wrapped every annotation value in a list, so annotation["sbo"] came back as ['SBO:...'] while add_sbo_terms and cobra's SBML reader give a bare string. sbo is single-valued and is now kept as a string, also when the file spells it as a one-item list. --- src/raven_toolbox/io/yaml.py | 13 +++++++---- tests/test_io_yaml.py | 45 ++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+), 4 deletions(-) diff --git a/src/raven_toolbox/io/yaml.py b/src/raven_toolbox/io/yaml.py index 09977b7..2be733a 100644 --- a/src/raven_toolbox/io/yaml.py +++ b/src/raven_toolbox/io/yaml.py @@ -239,7 +239,7 @@ def model_from_yaml_data(raw: dict) -> cobra.Model: # single-item list here, before model_from_dict, so cobra/geckopy always # see annotation[key] as list[str] regardless of whether the source file # used a scalar or a list — cobra's own model_from_dict does not do this - # normalisation itself. + # normalisation itself. ``sbo`` stays a bare string. for section in ("metabolites", "reactions", "genes"): _normalize_annotation_values(raw.get(section)) @@ -334,8 +334,10 @@ def _normalize_annotation_values(entries) -> None: store whatever shape the source used — means ``annotation[key]`` is always ``list[str]`` however the file spelled it, matching what cobra and geckopy expect and keeping a scalar-collapsed value from coming - back as a bare string after a write/read round trip. Normalises in - place; a no-op when ``entries`` is falsy. + back as a bare string after a write/read round trip. The exception is + ``sbo``, which is single-valued: it is kept as a bare string, the type + cobra's SBML reader and :func:`~raven_toolbox.annotation.add_sbo_terms` + give it. Normalises in place; a no-op when ``entries`` is falsy. """ if not entries: return @@ -346,7 +348,10 @@ def _normalize_annotation_values(entries) -> None: if not isinstance(annotation, dict): continue for key, value in annotation.items(): - if not isinstance(value, list): + if key == "sbo": + if isinstance(value, list) and len(value) == 1: + annotation[key] = value[0] + elif not isinstance(value, list): annotation[key] = [value] diff --git a/tests/test_io_yaml.py b/tests/test_io_yaml.py index a925497..54e4f74 100644 --- a/tests/test_io_yaml.py +++ b/tests/test_io_yaml.py @@ -345,3 +345,48 @@ def test_read_yaml_model_gives_integer_charges(tmp_path): sbml = tmp_path / "float_charges.xml" cobra.io.write_sbml_model(model, str(sbml)) assert cobra.io.read_sbml_model(str(sbml)).metabolites.get_by_id("a_c").charge == -1 + + +def test_sbo_term_is_a_string_after_yaml_round_trip(tmp_path): + """``annotation['sbo']`` is a bare string after a YAML round trip, as after add_sbo_terms or SBML read.""" + m = cobra.Model("sbo") + a = cobra.Metabolite("a_c", compartment="c") + a.annotation["sbo"] = "SBO:0000247" + a.annotation["kegg.compound"] = ["C00031"] + r = cobra.Reaction("r1", lower_bound=0, upper_bound=1000) + r.add_metabolites({a: -1}) + r.annotation["sbo"] = "SBO:0000627" + m.add_reactions([r]) + out = tmp_path / "sbo.yml" + write_yaml_model(m, out) + back = read_yaml_model(out) + assert back.reactions.get_by_id("r1").annotation["sbo"] == "SBO:0000627" + assert back.metabolites.get_by_id("a_c").annotation["sbo"] == "SBO:0000247" + assert back.metabolites.get_by_id("a_c").annotation["kegg.compound"] == ["C00031"] + + +def test_sbo_term_written_as_a_list_is_read_as_a_string(tmp_path): + path = tmp_path / "sbo_list.yml" + path.write_text( + "!!omap\n" + "- metabolites:\n" + " - !!omap\n" + " - id: a_c\n" + " - compartment: c\n" + "- reactions:\n" + " - !!omap\n" + " - id: r1\n" + " - metabolites: !!omap\n" + " - a_c: -1.0\n" + " - lower_bound: 0.0\n" + " - upper_bound: 1000.0\n" + " - annotation: !!omap\n" + " - sbo:\n" + " - SBO:0000627\n" + "- genes: []\n" + "- compartments: !!omap\n" + " - c: cytosol\n", + encoding="utf-8", + ) + model = read_yaml_model(path) + assert model.reactions.get_by_id("r1").annotation["sbo"] == "SBO:0000627" From 4ec79977b53ee4e12ae8c3a6dfd32b32e6358c23 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:39:44 +0200 Subject: [PATCH 02/17] fix: store batch-curated notes under the keys the YAML layer uses batch_curate kept metNotes, rxnNotes, rxnReferences and rxnConfidenceScores under their table names in .notes, and the InChI in annotation["inchi"]. read_yaml_model and the YAML/MAT writers use note, references, confidence_score and inchis, so write_yaml_model emitted a nested mapping that MATLAB read back as "metNotes: ...", and a model read from YAML ended up with two different InChIs. The values now go under those keys, the confidence score as a number, and an InChI already held as an annotation is updated too. --- src/raven_toolbox/curation/__init__.py | 7 +++ src/raven_toolbox/curation/batch.py | 50 ++++++++++++-------- tests/test_curation.py | 64 ++++++++++++++++++++++++-- 3 files changed, 98 insertions(+), 23 deletions(-) diff --git a/src/raven_toolbox/curation/__init__.py b/src/raven_toolbox/curation/__init__.py index 131f6e0..3ed6d35 100644 --- a/src/raven_toolbox/curation/__init__.py +++ b/src/raven_toolbox/curation/__init__.py @@ -37,6 +37,13 @@ Everything after the core columns in any of the four tables is interpreted as a MIRIAM annotation: the column header becomes the namespace key (``met.annotation[
] = ``). + +``inchi``, ``metNotes``, ``rxnNotes``, ``rxnReferences`` and +``rxnConfidenceScores`` are RAVEN fields and go to ``.notes`` under the keys +:func:`~raven_toolbox.io.read_yaml_model` uses for them (``inchis``, +``note``, ``note``, ``references`` and ``confidence_score``, the last as a +number), so the YAML and MAT writers emit them as those fields. An InChI the +metabolite already carries as ``annotation["inchi"]`` is updated as well. """ from raven_toolbox.curation.batch import ( DEFAULT_CORE_GENE_COLUMNS, diff --git a/src/raven_toolbox/curation/batch.py b/src/raven_toolbox/curation/batch.py index 5004349..7043931 100644 --- a/src/raven_toolbox/curation/batch.py +++ b/src/raven_toolbox/curation/batch.py @@ -184,10 +184,7 @@ def _update_metabolite(met: cobra.Metabolite, row: pd.Series, miriam_cols: list[ met.formula = str(row["formula"]) if _has_value(row.get("charge")): met.charge = _coerce_int(row["charge"]) - if _has_value(row.get("inchi")): - met.annotation["inchi"] = str(row["inchi"]) - if _has_value(row.get("metNotes")): - met.notes["metNotes"] = str(row["metNotes"]) + _apply_met_notes(met, row) _apply_miriam(met, row, miriam_cols) @@ -209,10 +206,7 @@ def _build_new_metabolites( met.formula = str(row["formula"]) if _has_value(row.get("charge")): met.charge = _coerce_int(row["charge"]) - if _has_value(row.get("inchi")): - met.annotation["inchi"] = str(row["inchi"]) - if _has_value(row.get("metNotes")): - met.notes["metNotes"] = str(row["metNotes"]) + _apply_met_notes(met, row) _apply_miriam(met, row, miriam_cols) new_mets.append(met) return new_mets @@ -374,12 +368,7 @@ def _update_reaction(rxn: cobra.Reaction, row: pd.Series, miriam_cols: list[str] rxn.subsystem = subsystem_to_str(row["subSystems"]) if _has_value(row.get("eccodes")): rxn.annotation["ec-code"] = str(row["eccodes"]) - if _has_value(row.get("rxnNotes")): - rxn.notes["rxnNotes"] = str(row["rxnNotes"]) - if _has_value(row.get("rxnReferences")): - rxn.notes["rxnReferences"] = str(row["rxnReferences"]) - if _has_value(row.get("rxnConfidenceScores")): - rxn.notes["rxnConfidenceScores"] = str(row["rxnConfidenceScores"]) + _apply_rxn_notes(rxn, row) _apply_miriam(rxn, row, miriam_cols) @@ -405,12 +394,7 @@ def _build_new_reactions( rxn.subsystem = subsystem_to_str(row["subSystems"]) if _has_value(row.get("eccodes")): rxn.annotation["ec-code"] = str(row["eccodes"]) - if _has_value(row.get("rxnNotes")): - rxn.notes["rxnNotes"] = str(row["rxnNotes"]) - if _has_value(row.get("rxnReferences")): - rxn.notes["rxnReferences"] = str(row["rxnReferences"]) - if _has_value(row.get("rxnConfidenceScores")): - rxn.notes["rxnConfidenceScores"] = str(row["rxnConfidenceScores"]) + _apply_rxn_notes(rxn, row) _apply_miriam(rxn, row, miriam_cols) new_rxns.append(rxn) return new_rxns @@ -418,6 +402,32 @@ def _build_new_reactions( # --- shared helpers --------------------------------------------------- +# Table columns stored in ``.notes``, mapped to the keys read_yaml_model and +# read_mat_model use for the same RAVEN fields (see _MET_FIELDS/_RXN_FIELDS in +# raven_toolbox.io.yaml), so the writers emit them as RAVEN fields. +_MET_NOTE_COLUMNS = (("inchi", "inchis"), ("metNotes", "note")) +_RXN_NOTE_COLUMNS = (("rxnNotes", "note"), ("rxnReferences", "references")) + + +def _apply_met_notes(met: cobra.Metabolite, row: pd.Series) -> None: + for col, key in _MET_NOTE_COLUMNS: + if _has_value(row.get(col)): + met.notes[key] = str(row[col]) + # An InChI already kept as a MIRIAM annotation is updated as well, so the + # metabolite does not carry two different structures. + if _has_value(row.get("inchi")) and "inchi" in met.annotation: + met.annotation["inchi"] = str(row["inchi"]) + + +def _apply_rxn_notes(rxn: cobra.Reaction, row: pd.Series) -> None: + for col, key in _RXN_NOTE_COLUMNS: + if _has_value(row.get(col)): + rxn.notes[key] = str(row[col]) + if _has_value(row.get("rxnConfidenceScores")): + score = float(row["rxnConfidenceScores"]) + rxn.notes["confidence_score"] = int(score) if score.is_integer() else score + + def _apply_miriam(entity, row: pd.Series, miriam_cols: list[str]) -> None: for col in miriam_cols: value = row.get(col) diff --git a/tests/test_curation.py b/tests/test_curation.py index e7f9bd0..8f0c4fa 100644 --- a/tests/test_curation.py +++ b/tests/test_curation.py @@ -50,7 +50,7 @@ def test_add_new_metabolite(): assert new.compartment == "c" assert new.formula == "C21H27N7O14P2" assert new.charge == -2 - assert new.notes.get("metNotes") == "added by curation" + assert new.notes.get("note") == "added by curation" def test_update_existing_metabolite(): @@ -66,9 +66,67 @@ def test_update_existing_metabolite(): atp = m.metabolites.get_by_id("s_0001") assert atp.formula == "C10H16N5O13P3" # overwritten assert atp.charge == -3 + assert atp.notes["inchis"] == "InChI=ATP" + + +def test_update_keeps_an_annotated_inchi_in_step(): + m = _base_model() + atp = m.metabolites.get_by_id("s_0001") + atp.annotation["inchi"] = "InChI=old" + df = pd.DataFrame([ + {"metNames": "ATP", "comps": "c", "formula": "", "charge": "", + "inchi": "InChI=ATP", "metNotes": ""}, + ]) + with pytest.warns(UserWarning, match="overwritten"): + batch_curate(m, mets_df=df, met_id_prefix="s_") + assert atp.notes["inchis"] == "InChI=ATP" assert atp.annotation["inchi"] == "InChI=ATP" +def test_curated_notes_are_written_as_raven_yaml_fields(tmp_path): + """Notes, references, confidence scores and InChIs reach YAML as RAVEN's own fields.""" + from cobra.io.yaml import yaml as cobra_yaml + + from raven_toolbox.io import read_yaml_model, write_yaml_model + + m = _base_model() + mets_df = pd.DataFrame([ + {"metNames": "ATP", "comps": "c", "formula": "", "charge": "", + "inchi": "InChI=ATP", "metNotes": "charge assigned at pH 7.3"}, + ]) + rxns_df = pd.DataFrame([ + {"rxnNames": "ATP hydrolysis", "grRules": "", "lb": "", "ub": "", + "rev": "", "subSystems": "", "eccodes": "", "rxnNotes": "curated", + "rxnReferences": "PMID:1", "rxnConfidenceScores": 2}, + ]) + coeffs_df = pd.DataFrame([ + {"rxnNames": "ATP hydrolysis", "metNames": "ATP", "comps": "c", + "coefficient": -1.0}, + {"rxnNames": "ATP hydrolysis", "metNames": "ADP", "comps": "c", + "coefficient": 1.0}, + ]) + with pytest.warns(UserWarning, match="overwritten"): + batch_curate(m, mets_df=mets_df, rxns_df=rxns_df, + rxns_coeffs_df=coeffs_df, met_id_prefix="s_", + rxn_id_prefix="r_") + out = tmp_path / "curated.yml" + write_yaml_model(m, out) + + with open(out, encoding="utf-8") as fh: + raw = cobra_yaml.load(fh) + met = next(x for x in raw["metabolites"] if x["id"] == "s_0001") + assert met["notes"] == "charge assigned at pH 7.3" + assert met["inchis"] == "InChI=ATP" + rxn = next(x for x in raw["reactions"] if x["id"] == "r_0001") + assert rxn["notes"] == "curated" + assert rxn["references"] == "PMID:1" + assert rxn["confidence_score"] == 2 + + back = read_yaml_model(out) + assert back.metabolites.get_by_id("s_0001").notes["note"] == "charge assigned at pH 7.3" + assert back.reactions.get_by_id("r_0001").notes["confidence_score"] == 2 + + def test_miriam_columns_auto_detected(): m = _base_model() df = pd.DataFrame([ @@ -144,7 +202,7 @@ def test_add_new_reaction_with_existing_mets(): assert new.gene_reaction_rule == "YBR456W" assert new.annotation["ec-code"] == "2.7.4.6" assert new.annotation["kegg.reaction"] == "R00187" - assert new.notes.get("rxnConfidenceScores") == "3" + assert new.notes.get("confidence_score") == 3 def test_add_new_reaction_joins_list_subsystems(): @@ -192,7 +250,7 @@ def test_update_existing_reaction_by_stoichiometry(): assert rxn.name == "ATP hydrolysis renamed" assert rxn.gene_reaction_rule == "YAL999W" assert rxn.lower_bound == -100 - assert rxn.notes.get("rxnNotes") == "updated" + assert rxn.notes.get("note") == "updated" def test_rxn_with_unknown_met_raises(): From 6b690feb868bb114c299db504103891a464507c2 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:40:47 +0200 Subject: [PATCH 03/17] fix: raise in scale_biomass when the balancing component cannot absorb the change With balance_out, scale_biomass rescaled the balancing component so the total is 1 g/gDW without checking that its new fraction stays positive. A negative factor turned the pseudoreaction's substrates into products and the model was changed without an error; a balancing component with fraction 0 was skipped silently. Both now raise before the model is changed, and the docstring says the total becomes 1 g/gDW rather than staying there. --- src/raven_toolbox/biomass/scale.py | 25 ++++++++++++++++++++++--- tests/test_biomass.py | 22 ++++++++++++++++++++++ 2 files changed, 44 insertions(+), 3 deletions(-) diff --git a/src/raven_toolbox/biomass/scale.py b/src/raven_toolbox/biomass/scale.py index ea32d08..31c93ae 100644 --- a/src/raven_toolbox/biomass/scale.py +++ b/src/raven_toolbox/biomass/scale.py @@ -93,11 +93,16 @@ def scale_biomass( balance_out: str | None = None, ) -> None: """Rescale a component to a new mass fraction, optionally balancing - out a second component to keep the total at 1 g/gDW. + out a second component so the biomass total becomes 1 g/gDW. Ports ``scaleBioMass.m``. The rescaling factor is derived from :func:`sum_biomass` on the *current* model state, so call this with an in-place ``model`` mutation, not on a stale snapshot. + + With ``balance_out``, that component is rescaled so the components sum + to exactly 1 g/gDW, whatever the total was before. It must be large + enough to absorb the difference: a ``ValueError`` is raised, before the + model is changed, when its fraction is 0 or would become 0 or negative. """ fractions = sum_biomass(model, config) current = fractions[component_name] @@ -106,6 +111,22 @@ def scale_biomass( f"Cannot scale {component_name!r} to {new_value}: current " "fraction is 0 (pseudoreaction missing or empty)." ) + if balance_out is not None: + balance_current = fractions[balance_out] + if balance_current == 0: + raise ValueError( + f"Cannot balance with {balance_out!r}: current fraction is 0 " + "(pseudoreaction missing or empty)." + ) + # Fractions are linear in the substrate coefficients, so the total + # after the first rescaling is known before the model is changed. + balance_new = balance_current + 1.0 - (fractions["total"] - current + new_value) + if balance_new <= 0: + raise ValueError( + f"Cannot balance {component_name!r} = {new_value:g} with " + f"{balance_out!r}: it would have to go from {balance_current:.4g} " + f"to {balance_new:.4g} g/gDW to bring the total to 1." + ) factor = new_value / current rescale_pseudoreaction(model, config, component_name, factor) @@ -114,8 +135,6 @@ def scale_biomass( fractions = sum_biomass(model, config) total = fractions["total"] balance_current = fractions[balance_out] - if balance_current == 0: - return # nothing we can do balance_factor = (balance_current + (1.0 - total)) / balance_current rescale_pseudoreaction(model, config, balance_out, balance_factor) diff --git a/tests/test_biomass.py b/tests/test_biomass.py index da3e31a..288dc1a 100644 --- a/tests/test_biomass.py +++ b/tests/test_biomass.py @@ -207,6 +207,28 @@ def test_scale_biomass_zero_current_raises(): scale_biomass(m, cfg2, "DNA", 0.1) +def test_scale_biomass_balance_that_would_go_negative_raises(): + m, cfg = _toy_biomass_model() + before = {r.id: dict(r.metabolites) for r in m.reactions} + # Protein above 1 g/gDW leaves carbohydrate a negative share of the total, + # which would turn its substrates into products. + with pytest.raises(ValueError, match="carbohydrate"): + scale_biomass(m, cfg, "protein", 1.5, balance_out="carbohydrate") + assert {r.id: dict(r.metabolites) for r in m.reactions} == before + + +def test_scale_biomass_balance_on_empty_component_raises(): + m, cfg = _toy_biomass_model() + cfg2 = BiomassConfig.from_components( + cfg.biomass_rxn, cfg.proton_met, + [*cfg.components, BiomassComponent("DNA", "DNA pseudoreaction")], + ) + before = {r.id: dict(r.metabolites) for r in m.reactions} + with pytest.raises(ValueError, match="DNA"): + scale_biomass(m, cfg2, "protein", 0.5, balance_out="DNA") + assert {r.id: dict(r.metabolites) for r in m.reactions} == before + + # --- set_gam -------------------------------------------------------- def test_set_gam_scales_cofactor_coefficients(): From 2f75eed9dc83a9dca7fa4cb796fff9cfdae29958 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:43:44 +0200 Subject: [PATCH 04/17] fix: stop fill_gaps_kumar_milp reporting a repair that does not restore growth The indicator coupling uses big-M 1000, so with the solver's default integrality tolerance (about 1e-5) a "closed" candidate could carry about 0.01 flux, enough for the default 10 % growth floor. On smallYeast this returned "optimal" with nothing added for the RPI, DAR, CIT, ACO, IDP1, CO2TRANS and GROWTH removals on GLPK. The MILP now runs with an integrality tolerance of 1e-9, as RAVEN's optimizeProb does, and the rounded selection is checked on the filled model: when it does not reach min_growth the result has exit_status "not_restored". --- src/raven_toolbox/gapfilling/kumar_milp.py | 50 +++++++++++++++++++++- tests/test_gapfilling_kumar_milp.py | 34 +++++++++++++++ 2 files changed, 83 insertions(+), 1 deletion(-) diff --git a/src/raven_toolbox/gapfilling/kumar_milp.py b/src/raven_toolbox/gapfilling/kumar_milp.py index d2e5dfc..765daaf 100644 --- a/src/raven_toolbox/gapfilling/kumar_milp.py +++ b/src/raven_toolbox/gapfilling/kumar_milp.py @@ -24,6 +24,7 @@ """ from __future__ import annotations +import contextlib import math from collections.abc import Iterable from dataclasses import dataclass, field @@ -34,6 +35,9 @@ from raven_toolbox.manipulation.transfer import add_reactions_from_model +# Integrality tolerance for the MILP, as RAVEN's optimizeProb sets IntFeasTol. +_INTEGRALITY_TOL = 1e-9 + @dataclass class KumarGapFillResult: @@ -50,7 +54,10 @@ class KumarGapFillResult: Draft model with *added_reactions* incorporated and *reversed_reactions* bounds updated. exit_status: - Solver status string (e.g. ``"optimal"``, ``"infeasible"``). + Solver status string (e.g. ``"optimal"``, ``"infeasible"``), or + ``"not_restored"`` when the solver reported a solution but the + selected repairs do not bring the objective to *min_growth* in + *model*. The selection is still returned for inspection. """ added_reactions: list[str] = field(default_factory=list) @@ -250,6 +257,13 @@ def fill_gaps_kumar_milp( f"({len(rev_indicators)} reversal, {len(db_indicators)} database candidates)..." ) + # A fractional indicator within the solver's default integrality tolerance + # (about 1e-5 on GLPK and Gurobi) times big_m lets a "closed" candidate + # carry enough flux to meet a small growth floor. RAVEN's optimizeProb uses + # the same tight tolerance for MILPs. + with contextlib.suppress(Exception): + working.solver.configuration.tolerances.integrality = _INTEGRALITY_TOL + working.optimize() status = working.solver.status @@ -289,9 +303,43 @@ def fill_gaps_kumar_milp( add_reactions_from_model(filled, t, ids_from_t, genes=False, note=None) remaining -= set(ids_from_t) + # The rounded selection must meet the floor on its own; a solution that + # relied on fractional indicators does not. + if not _meets_floor(filled, obj_coeffs, min_growth): + if verbose: + print( + "fill_gaps_kumar_milp: the selected repairs do not reach " + f"min_growth {min_growth:.4g} in the filled model." + ) + status = "not_restored" + return KumarGapFillResult( added_reactions=sorted(added_rxn_ids), reversed_reactions=sorted(reversed_rxn_ids), model=filled, exit_status=status, ) + + +def _meets_floor(model: cobra.Model, obj_coeffs: dict, min_growth: float) -> bool: + """Whether ``model`` is feasible with the MILP's growth floor applied. + + The floor is the same bound the MILP puts on each objective reaction, + relaxed by a relative 1e-6 for the LP feasibility tolerance. + """ + floor = min_growth * (1 - 1e-6) + with model: + for rxn, coeff in obj_coeffs.items(): + if rxn.id not in model.reactions: + return False + target = model.reactions.get_by_id(rxn.id) + if coeff > 0: + if target.upper_bound < floor: + return False + target.lower_bound = max(target.lower_bound, floor) + elif coeff < 0: + if target.lower_bound > -floor: + return False + target.upper_bound = min(target.upper_bound, -floor) + model.slim_optimize(error_value=math.nan) + return model.solver.status == "optimal" diff --git a/tests/test_gapfilling_kumar_milp.py b/tests/test_gapfilling_kumar_milp.py index 6f1bd42..23a5533 100644 --- a/tests/test_gapfilling_kumar_milp.py +++ b/tests/test_gapfilling_kumar_milp.py @@ -129,3 +129,37 @@ def test_exit_status_on_infeasible(): template = cobra.Model("empty_template") result = fill_gaps_kumar_milp(m, template, verbose=False) assert result.exit_status != "optimal" + + +def _solvers(): + out = ["glpk"] + try: + import gurobipy # noqa: F401 + out.append("gurobi") + except ImportError: + pass + return out + + +@pytest.mark.parametrize("solver", _solvers()) +def test_small_floor_is_not_met_by_a_fractional_indicator(growth_gap, solver): + """With a floor below 1e-5 x big-M, a near-zero indicator must not pass as 'not added'.""" + draft, template = growth_gap + draft.solver = solver + result = fill_gaps_kumar_milp(draft, template, min_growth=0.005, verbose=False) + assert result.exit_status == "optimal" + assert result.added_reactions == ["r_missing"] + assert result.model.slim_optimize() >= 0.005 + + +def test_selection_that_does_not_restore_growth_is_reported(growth_gap, monkeypatch): + """A rounded MILP selection is checked on the filled model.""" + from raven_toolbox.gapfilling import kumar_milp + + draft, template = growth_gap + draft.solver = "glpk" + # A loose integrality tolerance lets the indicator stay fractional. + monkeypatch.setattr(kumar_milp, "_INTEGRALITY_TOL", 1e-4) + result = fill_gaps_kumar_milp(draft, template, min_growth=0.005, verbose=False) + assert result.added_reactions == [] + assert result.exit_status == "not_restored" From be4f4588b2956646e24c94fc16775caebb9f6df5 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:47:22 +0200 Subject: [PATCH 05/17] docs: describe the reporter_metabolites background correction and RAVEN parity as implemented The docstring credited the parity to replacing RAVEN's Monte-Carlo background, which RAVEN no longer uses, said RAVEN sorts by p-value and reports both tails, and cited an IMPROVEMENTS entry that does not exist. It now states the population standard deviation, the one-sided p-value and the z-score sort, and the constant factor a RAVEN version with the sample standard deviation differs by. --- src/raven_toolbox/analysis/reporter.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/raven_toolbox/analysis/reporter.py b/src/raven_toolbox/analysis/reporter.py index 5150968..4d9e911 100644 --- a/src/raven_toolbox/analysis/reporter.py +++ b/src/raven_toolbox/analysis/reporter.py @@ -92,15 +92,17 @@ def reporter_metabolites( is given, two extra results are returned for the up- (fc ≥ 0) and down- (fc < 0) regulated gene subsets, in addition to ``"all"``. - Parity with RAVEN's ``reporterMetabolites``: the ``z_score`` and underlying - background correction match exactly (exact closed-form instead of RAVEN's - Monte-Carlo, see IMPROVEMENTS RM1). The reported ``p_value`` is the - *one-sided* (``"up"``) enrichment ``1 - Φ(z)`` and the result is sorted by - ``z_score`` descending. RAVEN sorts by p-value and reports both tails - (``allPValues``, ``allUpPValues``, ``allDownPValues``); the up/down splits - here come from the ``gene_fold_changes`` subset partition instead, so the - same information is available via the three returned ``ReporterResult`` - rows. + The background correction uses the population standard deviation + (``ddof=0``) of the scored Z-scores, which the closed form requires. The + reported ``p_value`` is the one-sided enrichment ``1 - Φ(z)`` and each + table is sorted by ``z_score`` descending. + + Parity with RAVEN's ``reporterMetabolites``: the same closed-form + correction, one-sided p-value and sort order, with ``"up"``/``"down"`` + matching its ``"only up"``/``"only down"`` entries. A RAVEN version that + divides by the sample standard deviation (``std`` with N−1) gives + z-scores smaller by a factor ``√(N/(N−1))`` for N scored genes; the + ranking is the same. """ model_genes = {g.id for g in model.genes} scored = { From 0ca92c43d1f50cb97cab8cccf738f69d5110da96 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:47:22 +0200 Subject: [PATCH 06/17] docs: bring the sampling module and random_sampling docstrings up to date The module docstring presented the module as random-objective sampling only, said RAVEN has no seed and hard-codes two objectives, and cited an IMPROVEMENTS entry that does not exist. It now describes the method dispatch and the actual differences from RAVEN's random-objective method. The convergence study link points to CompMetEng/raven-docs. --- src/raven_toolbox/analysis/sampling.py | 68 +++++++++++++------------- 1 file changed, 35 insertions(+), 33 deletions(-) diff --git a/src/raven_toolbox/analysis/sampling.py b/src/raven_toolbox/analysis/sampling.py index b7cc230..2d9b224 100644 --- a/src/raven_toolbox/analysis/sampling.py +++ b/src/raven_toolbox/analysis/sampling.py @@ -1,31 +1,33 @@ -"""Random-objective flux sampling — RAVEN's ``randomSampling`` (port + improvements). - -Samples the flux solution space by the **random-objective** method of Bordel et al. -(2010, PLoS Comput Biol, doi:10.1371/journal.pcbi.1000859), as ported from RAVEN's -``randomSampling``: each sample maximises a small random linear combination of -reactions, so every sample is an *extreme point* (vertex) of the flux polytope. - -This is a different statistical object from cobrapy's ``cobra.sampling`` (OptGP / -ACHR), which draw a (near-)uniform Markov-chain sample of the polytope *interior*. -Use cobra's samplers when you need the uniform flux distribution; use this when you -want a fast, robust spread of diverse optimal states — the workflow RAVEN uses to -compare conditions, and one that stays well-behaved on large or tightly-constrained -models where MCMC mixing is poor. cobrapy has no equivalent, so this is a genuine -addition, not a wrapper. - -Improvements over RAVEN (see IMPROVEMENTS SAMP1): - -* **`good_reactions` via one FVA pass**, not a hand-rolled per-reaction ``parfor`` - loop. A reaction is usable as a random objective if it can carry flux and is not - stuck in a stoichiometrically-infeasible loop (its range blows past the arbitrary - large bound). ``cobra``'s FVA computes exactly that, faster and in far less code, - and can optionally be made ``loopless``. -* **Reproducible** via ``seed`` (RAVEN has no seed control). -* **`n_objectives` is a parameter** (RAVEN hard-codes 2, though its docstring claims - 3). -* **Tidy output**: a ``samples`` DataFrame shaped samples × reactions (matching - ``cobra.sampling``), plus the reusable ``good_reactions`` list — instead of a - reactions × samples matrix and a parallel index vector. +"""Flux sampling — RAVEN's ``randomSampling``. + +:func:`random_sampling` is the entry point for every sampling method and, like +RAVEN's ``randomSampling``, dispatches on ``method``: + +* ``"achr"`` (default) and ``"chrr"`` draw a (near-)uniform Markov-chain sample + of the polytope interior; they are implemented in + :mod:`raven_toolbox.analysis.flux_sampling`. +* ``"random_objective"`` is the random-objective method of Bordel et al. (2010, + PLoS Comput Biol, doi:10.1371/journal.pcbi.1000859), implemented here: each + sample maximises a small random linear combination of reactions, so every + sample is an *extreme point* (vertex) of the flux polytope. cobrapy has no + equivalent. Use it for a spread of diverse optimal states, not for the + uniform flux distribution. + +Differences from RAVEN's random-objective method: + +* **Reactions usable as random objectives** (``good_reactions``) come from one + FVA pass, loopless by default (:func:`find_good_reactions`), which drops + reactions that carry flux only through a loop. RAVEN runs + ``getAllowedBounds`` and drops reactions that reach 0.999 × the model's + largest bound, and on GECKO ecModels also the ``usage_prot_`` reactions. +* **``replace_max_bound``** defaults to False; RAVEN's ``replaceBoundsWithInf`` + defaults to true. +* **Output**: a ``samples`` DataFrame shaped samples × reactions (matching + ``cobra.sampling``), plus the reusable ``good_reactions`` list, instead of a + reactions × samples matrix and a vector of reaction indexes. + +``seed`` and ``n_objectives`` (default 2) behave as RAVEN's ``seed`` and +``nObjectives``. """ from __future__ import annotations @@ -113,7 +115,7 @@ def random_sampling( n_samples to compensate, try ``method='chrr'`` (raven-toolbox's rounding-based alternative, not wired into ACHR), or use an ESS/R-hat diagnostic to assess sample quality (see the `sampling convergence calibration study - `_ + `_ on raven-docs). warmup: ``chrr`` — burn-in steps discarded before the first recorded sample. @@ -134,10 +136,9 @@ def random_sampling( ``random_objective`` only — worker processes for drawing samples in parallel (each sample is an independent LP solve). Defaults to ``cobra.Configuration().processes``; set to 1 to sample serially. - Reproducible for a given ``seed`` regardless of ``n_proc`` — each + Reproducible for a given ``seed`` regardless of ``n_proc``: each sample draws from its own independent RNG stream spawned from - ``seed``, not from one continuously-advancing generator, so results - no longer match pre-parallelization output for the same seed value. + ``seed``, not from one continuously-advancing generator. Returns ------- @@ -336,7 +337,8 @@ def _random_objective( n_samples Number of flux vectors to return. n_objectives - Reactions combined into each random objective (RAVEN's fixed 2). + Reactions combined into each random objective (default 2, as RAVEN's + ``nObjectives``). good_reactions Reaction ids eligible as objectives. If ``None`` they are computed once with :func:`find_good_reactions` and returned for reuse. From 64c7f34ebfb743e09a5484c2cdff4e053c032cd6 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:48:28 +0200 Subject: [PATCH 07/17] docs: say which bounds remove_no_flux_reactions tests The docstrings said it removes reactions that cannot carry flux in any steady state and never one that merely lacks an open boundary. Only model.exchanges are opened; sinks, demands and internal bounds keep their values. --- src/raven_toolbox/manipulation/simplify.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/raven_toolbox/manipulation/simplify.py b/src/raven_toolbox/manipulation/simplify.py index f363531..2457398 100644 --- a/src/raven_toolbox/manipulation/simplify.py +++ b/src/raven_toolbox/manipulation/simplify.py @@ -8,9 +8,10 @@ (one that, given reaction directions, can only be produced or only consumed, or participates in a single reaction); RAVEN ``deleteInaccessible``. Iterates to a fixpoint and prunes orphaned metabolites. -* ``remove_no_flux_reactions`` — reactions that cannot carry flux in *any* steady state, - found by FVA (both min and max flux zero); RAVEN ``deleteMinMax``. Catches blocked - reactions that are not topological dead-ends. +* ``remove_no_flux_reactions`` — reactions that cannot carry flux under the model's + bounds (exchange reactions opened by default), found by FVA (both min and max flux + zero); RAVEN ``deleteMinMax``. Catches blocked reactions that are not topological + dead-ends. * ``remove_duplicate_reactions`` — all-but-one of each set of reactions with identical stoichiometry/bounds; RAVEN ``deleteDuplicates`` (detection-only: ``find_duplicate_reactions``). * ``constrain_reversible_reactions`` — does not remove reactions; tightens bounds, making @@ -131,11 +132,14 @@ def remove_no_flux_reactions( """Remove reactions that cannot carry any flux (RAVEN ``simplifyModel`` ``deleteMinMax``). Runs FVA and drops every reaction whose minimum and maximum flux are both zero. - With ``open_exchanges`` (default) boundary reactions are opened first, so only - *structurally* blocked reactions are removed, never ones that merely lack an open - boundary. A no-flux reaction cannot carry flux under any tighter constraint either, - so removing it before task discovery / the merge is safe. Modifies in place; returns - the removed reaction ids. + With ``open_exchanges`` (default) the reactions in ``model.exchanges`` are first + opened to at least ``(-1000, 1000)``; sink and demand reactions and all internal + bounds keep their current values, so a reaction blocked only by one of those is + still removed. With ``open_exchanges=False`` the current bounds are used, as RAVEN's + ``deleteMinMax`` does, which can remove more. A reaction that is blocked under these + bounds stays blocked under any tighter ones, so removing it before task discovery / + the merge is safe as long as no later step relaxes a sink, demand or internal + bound. Modifies in place; returns the removed reaction ids. """ blocked = find_blocked_reactions(model, open_exchanges=open_exchanges) if blocked: From 507259810ef3f81ecba4c66f87faaade5b826a19 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:48:28 +0200 Subject: [PATCH 08/17] docs: qualify the RAVEN parity of the check_tasks close_boundaries default RAVEN's checkTasks leaves the model's exchanges open only on a model without boundary metabolites; prepINITModel calls closeModel first, which corresponds to close_boundaries=True. --- src/raven_toolbox/tasks/check.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/raven_toolbox/tasks/check.py b/src/raven_toolbox/tasks/check.py index bcf293e..51807b4 100644 --- a/src/raven_toolbox/tasks/check.py +++ b/src/raven_toolbox/tasks/check.py @@ -10,10 +10,11 @@ and an output allows / requires net production (``Sv ≤ UB``, and ``≥ LB`` if ``LB > 0``). By default, existing boundary reactions are left exactly as they are: a task's inputs/outputs *add to* whatever the model's own open exchanges already allow, -rather than replacing them (RAVEN's ``checkTasks`` semantics, and what the published -task lists this module reads were curated against). Pass ``close_boundaries=True`` for -the stricter reading, where a task's declared inputs/outputs are the complete boundary -of the system for that check. +rather than replacing them. That is what RAVEN's ``checkTasks`` does on a model without +boundary metabolites, such as one read from YAML. Pass ``close_boundaries=True`` for the +stricter reading, where a task's declared inputs/outputs are the complete boundary of the +system for that check; it corresponds to RAVEN's ``checkTasks(closeModel(model), ...)``, +which is how ``prepINITModel`` runs it. """ from __future__ import annotations @@ -212,9 +213,10 @@ def check_tasks( ``tasks`` is a parsed list of :class:`Task` or a path to a task-list file. By default (``close_boundaries=False``), existing exchange/sink/demand reactions are left open, so a task's inputs/outputs add to whatever the model's own boundary - already allows — matching RAVEN's ``checkTasks``, which published task lists were - curated against. Pass ``close_boundaries=True`` for the stricter reading, where - inputs/outputs are defined purely by the tasks. + already allows, as RAVEN's ``checkTasks`` does on a model without boundary + metabolites. Pass ``close_boundaries=True`` for the stricter reading, where + inputs/outputs are defined purely by the tasks, as with RAVEN's + ``checkTasks(closeModel(model), ...)``. """ tasks = _as_tasks(tasks) base, name_to_id, comp_to_ids = _prepare_base(model, close_boundaries) @@ -387,8 +389,10 @@ def find_task_essential_reactions( kept regardless of expression score and made irreversible in their forced direction. When a reaction is essential in several tasks with conflicting directions, the majority wins (ties → forward), matching RAVEN's ``pos < neg``. ``close_boundaries`` defaults to - ``False``, matching RAVEN: a task's inputs/outputs add to whatever the model's own - open exchanges already allow, rather than replacing them. + ``False``: a task's inputs/outputs add to whatever the model's own open exchanges + already allow, rather than replacing them, as in :func:`check_tasks`. RAVEN's + ``prepINITModel`` calls ``closeModel`` first, which corresponds to + ``close_boundaries=True``. The solver is configured per task to match RAVEN's ``solveLP`` (``FeasibilityTol = 1e-9``, single-threaded with a fixed seed) so the degenerate min-flux vertex — and thus From ae79c403bb57199304fcdc2d76de018354da884f Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:48:48 +0200 Subject: [PATCH 09/17] docs: state that blocked reactions are blocked under the tested bounds connect_blocked_reactions and the remove_no_flux_reactions summary described blocked reactions without saying they are found under the model's current (or exchange-opened) bounds. --- src/raven_toolbox/gapfilling/fill.py | 2 +- src/raven_toolbox/manipulation/simplify.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/raven_toolbox/gapfilling/fill.py b/src/raven_toolbox/gapfilling/fill.py index 1c22f36..bdb849d 100644 --- a/src/raven_toolbox/gapfilling/fill.py +++ b/src/raven_toolbox/gapfilling/fill.py @@ -138,7 +138,7 @@ def connect_blocked_reactions( """Add template reactions so blocked draft reactions can carry flux. Finds reactions that - cannot carry flux in ``model``, then adds the minimum-penalty set of template + cannot carry flux in ``model`` under its current bounds, then adds the minimum-penalty set of template reactions that lets the (irreversible) ones carry flux, and returns the filled model. Like RAVEN, only irreversible blocked reactions are forced — reversible ones can carry flux trivially in the split formulation, so forcing them is diff --git a/src/raven_toolbox/manipulation/simplify.py b/src/raven_toolbox/manipulation/simplify.py index 2457398..3f0a35a 100644 --- a/src/raven_toolbox/manipulation/simplify.py +++ b/src/raven_toolbox/manipulation/simplify.py @@ -129,7 +129,7 @@ def remove_zero_interval_reactions(model: cobra.Model) -> list[str]: def remove_no_flux_reactions( model: cobra.Model, *, open_exchanges: bool = True ) -> list[str]: - """Remove reactions that cannot carry any flux (RAVEN ``simplifyModel`` ``deleteMinMax``). + """Remove reactions that cannot carry flux under the tested bounds (RAVEN ``deleteMinMax``). Runs FVA and drops every reaction whose minimum and maximum flux are both zero. With ``open_exchanges`` (default) the reactions in ``model.exchanges`` are first From 0c82f655a22b9a459a7b2eb33fe50e8d088bc4c1 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:49:11 +0200 Subject: [PATCH 10/17] docs: say that find_leak_metabolite needs a closed model and when gap_report raises find_leak_metabolite uses the model's boundary reactions as they are, so the caller has to close them, and an infeasible search does not show the model is leak-free. gap_report raises Infeasible when the model, or its closed copy for the mass-balancing section, is infeasible. --- src/raven_toolbox/gapfilling/leak.py | 6 ++++++ src/raven_toolbox/gapfilling/report.py | 8 ++++++++ 2 files changed, 14 insertions(+) diff --git a/src/raven_toolbox/gapfilling/leak.py b/src/raven_toolbox/gapfilling/leak.py index 45c0db8..7dbd984 100644 --- a/src/raven_toolbox/gapfilling/leak.py +++ b/src/raven_toolbox/gapfilling/leak.py @@ -135,6 +135,12 @@ def find_leak_metabolite( ---------- model: Model to test. Never modified; the search runs on an internal copy. + Its boundary reactions are used as they are: close every exchange, + sink and demand reaction first (``for r in model.boundary: r.bounds = + (0, 0)`` on a copy), or an open uptake or excretion counts as a leak. + Positive lower bounds (such as on a maintenance reaction) still + apply and can make the search infeasible; the result then has status + ``"infeasible"``, which does not show that the model is leak-free. direction: ``"produce"`` to look for freely excreted metabolites, or ``"consume"`` to look for freely consumed metabolites. diff --git a/src/raven_toolbox/gapfilling/report.py b/src/raven_toolbox/gapfilling/report.py index 65b2b5b..9365f96 100644 --- a/src/raven_toolbox/gapfilling/report.py +++ b/src/raven_toolbox/gapfilling/report.py @@ -277,6 +277,14 @@ def gap_report( Returns ------- GapReportResult + + Raises + ------ + cobra.exceptions.Infeasible + If ``model`` is infeasible under its own bounds, or once its boundary + reactions are closed for ``can_produce_without_input`` and + ``can_consume_without_output`` (for example when a maintenance + reaction has a positive lower bound). """ no_flux_reactions = [r.id for r in model.reactions if r.id in _blocked_ids(model)] no_flux_reactions_relaxed = [r.id for r in model.reactions if r.id in _blocked_ids(_relaxed_copy(model))] From 70f58209e6fd16d797e95d217704f801485e198c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:49:50 +0200 Subject: [PATCH 11/17] docs: fix the README documentation claims and links, and the CONTRIBUTING docstring convention The README linked the MkDocs URL of the RAVEN 3 comparison page, which gives 404, and said every example on raven-docs is executed on each commit. CONTRIBUTING said all docstrings are NumPy style because the site's generator expects it; most are prose and the generator does not parse sections. raven-docs links point to CompMetEng. --- CONTRIBUTING.md | 7 ++++--- README.md | 8 ++++---- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index f7338e1..760086e 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -1,6 +1,6 @@ # Contributing to raven-toolbox -User documentation lives on [raven-docs](https://github.com/edkerk/raven-docs), +User documentation lives on [raven-docs](https://github.com/CompMetEng/raven-docs), which covers both this package and MATLAB RAVEN. This file is for working on the package itself, and for the maintainer tasks that have nowhere else to live. @@ -32,7 +32,8 @@ documentation belongs in the docstrings: raven-docs generates its API reference from them at build time, so a docstring edit reaches the site without anything being copied across. -Docstrings are NumPy style, which is what the site's generator expects. +Write new and changed public docstrings in NumPy style (Parameters, Returns and, where +useful, Examples sections). Many existing docstrings are still prose. --- @@ -45,7 +46,7 @@ consumes the published output. End users do none of this. They download the published, version-pinned artefacts automatically; that side is documented in -[Downloaded data and binaries](https://github.com/edkerk/raven-docs/blob/main/docs/installation/data-and-binaries.md). +[Downloaded data and binaries](https://github.com/CompMetEng/raven-docs/blob/main/docs/installation/data-and-binaries.md). Publishing the results, and the versioning model behind the releases, is a shared convention rather than a raven-toolbox concern, and is documented in diff --git a/README.md b/README.md index 8b47b98..e3b02a4 100644 --- a/README.md +++ b/README.md @@ -27,7 +27,7 @@ functionality that is unique to RAVEN: The status of every RAVEN function (ported, mapped to a cobra equivalent, or explicitly not ported) is documented function-by-function in -**[RAVEN 3 and raven-toolbox](https://raven-docs.readthedocs.io/en/latest/raven3-vs-raven-toolbox/)**. +**[RAVEN 3 and raven-toolbox](https://raven-docs.readthedocs.io/en/latest/raven3-vs-raven-toolbox.html)**. ## Design principle @@ -44,7 +44,7 @@ COBRA ecosystem. RAVEN 3.0 beta ([release notes](https://github.com/CompMetEng/raven-toolbox/releases/tag/v3.0.0b1), [changelog](CHANGELOG.md)). It has been validated against MATLAB RAVEN on **Human-GEM** (5 Hart2015 cell-line models, Jaccard 0.975–0.980 — see -[the Human-GEM validation study](https://github.com/edkerk/raven-docs/blob/main/docs/parameter-tuning/studies/humangem-validation.md) +[the Human-GEM validation study](https://github.com/CompMetEng/raven-docs/blob/main/docs/parameter-tuning/studies/humangem-validation.md) on raven-docs). Two deliberate scope decisions, not pending work: @@ -81,7 +81,7 @@ pip install -e ".[dev]" ``` raven-toolbox requires Python ≥ 3.11. Genome-scale ftINIT MILPs currently require **Gurobi** -([details on solver portability](https://github.com/edkerk/raven-docs/blob/main/docs/parameter-tuning/studies/init-solver-benchmark.md) +([details on solver portability](https://github.com/CompMetEng/raven-docs/blob/main/docs/parameter-tuning/studies/init-solver-benchmark.md) on raven-docs); toy and unit-test work runs on the open-source GLPK. ### External command-line tools (BLAST, DIAMOND, HMMER, MAFFT, CD-HIT) @@ -122,7 +122,7 @@ command, which is not downloaded automatically and runs on Linux and macOS: [**raven-docs**](https://raven-docs.readthedocs.io/) is the documentation, and it covers MATLAB RAVEN and raven-toolbox side by side: the same task shown in both -languages, with every example executed on each commit. +languages, with the examples run in both languages by the documentation's CI. This repository carries no prose documentation of its own. The API reference on the site is generated from the docstrings here at build time, so a docstring edit From 62d527432f6436af0c7ec7358c0967701fca2d0b Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:49:54 +0200 Subject: [PATCH 12/17] docs: point the raven-docs study links in ftinit and homology to CompMetEng --- src/raven_toolbox/init/ftinit.py | 18 +++++++++--------- .../reconstruction/homology/homology.py | 2 +- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index c4080cf..c7052b2 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -32,11 +32,11 @@ Needs a MILP solver (cobra's configured optlang solver; only Gurobi is fully viable at genome scale — see the `INIT solver benchmark -`_ +`_ on raven-docs). Magic numbers (``force_on``/``force_on_ess`` = 0.1, ``big_m`` = 100) are exposed and scale-dependent; calibration tables are in the `INIT parameter calibration study -`_ +`_ on raven-docs. ``big_m`` caps a *scored* reaction's flux in its on/off (direction) constraint — using a fixed 100 rather than the reaction's ±1000 bound keeps the LP relaxation tight (what makes the genome-scale @@ -150,7 +150,7 @@ def run_ftinit( — solve time here is dominated by model construction, so a tight gap is nearly free (measured on genome-scale Human-GEM, see the `INIT parameter calibration study - `_ + `_ on raven-docs). For the genome-scale staged schedule (:func:`raven_toolbox.init.ftinit`), an essential-forced step can run away without a cap (one severely-degraded case ran >75 min unbounded); @@ -214,7 +214,7 @@ def run_ftinit( knowingly trades determinism for speed. Measured to have no effect on the *swing* between different seeds (raising it is not a substitute for ``resolve_ties``) — see the `ftINIT reproducibility study - `_ + `_ on raven-docs. """ scores = dict(rxn_scores or {}) @@ -547,7 +547,7 @@ def _resolve_ties(opt, prob, obj_expr, indicators, primary, time_limit) -> bool: reduces the run-to-run spread relative to no tie-break at all, but it is no longer a guarantee. A warning is raised when this happens; see the `ftINIT reproducibility study - `_ + `_ on raven-docs. Returns ``True`` if phase 2 produced a usable solution (the caller then reads the @@ -805,12 +805,12 @@ def ftinit( incorrectly caps some reversible/negative reactions above magnitude 100; raven-toolbox never had this because its own ``big_m=100`` default relies on the model being rescaled first, not on a per-reaction cap — see the `INIT parameter calibration study - `_). + `_). ``mip_gap``/``time_limit`` are forwarded to each :func:`run_ftinit` solve. On genome-scale models they are essential for tractability — see the `INIT parameter calibration study - `_ + `_ on raven-docs for the calibration table. ``resolve_ties`` and ``prove_abs_gap`` (both opt-in, default off → exact RAVEN @@ -838,7 +838,7 @@ def ftinit( ``seed`` nor ``threads`` makes the extraction *stable* under small input changes, which is a separate property — see the `ftINIT reproducibility study - `_ + `_ on raven-docs. Measured there on Human-GEM/DLD1: ``resolve_ties`` makes repeated builds at a fixed seed identical and halves the seed-to-seed spread in predicted essential genes (16 → 8 flips); its own tie-break phases can exhaust ``time_limit`` at genome @@ -858,7 +858,7 @@ def ftinit( essential-gene drift on another, traced to a network-topology regime swap the reference-matching objective cannot see coming, and was removed. See the `reference_reactions postmortem - `_ + `_ on raven-docs for the full account. """ steps = steps if steps is not None else get_init_steps(series) diff --git a/src/raven_toolbox/reconstruction/homology/homology.py b/src/raven_toolbox/reconstruction/homology/homology.py index f516560..5d1094f 100644 --- a/src/raven_toolbox/reconstruction/homology/homology.py +++ b/src/raven_toolbox/reconstruction/homology/homology.py @@ -274,7 +274,7 @@ def get_model_from_homology( ``ValueError``; hits involving other organisms are ignored. Defaults for the three filters come from a `calibration against KEGG orthology - `_ + `_ on raven-docs, scored with precision weighted above recall (a wrongly transferred reaction is harder to undo than a missing one): From 36e95039fd75bb75a2d3c01c0f107206956cc21c Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 08:50:04 +0200 Subject: [PATCH 13/17] docs: correct the Kumar et al. citation and describe the reversal candidates in fill_gaps_kumar_milp The method is from BMC Bioinformatics 8:212. The module docstring now says that every irreversible draft reaction, exchanges and maintenance drains included, is a reversal candidate. --- src/raven_toolbox/gapfilling/kumar_milp.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/raven_toolbox/gapfilling/kumar_milp.py b/src/raven_toolbox/gapfilling/kumar_milp.py index 765daaf..8653d46 100644 --- a/src/raven_toolbox/gapfilling/kumar_milp.py +++ b/src/raven_toolbox/gapfilling/kumar_milp.py @@ -1,13 +1,16 @@ """Objective-based gap-filling via a global growth-floor MILP. Implements the single-level MILP formulation from Kumar et al. 2007 -(Bioinformatics 23:1626–1635) to find the minimum-cost set of model +(BMC Bioinformatics 8:212, doi:10.1186/1471-2105-8-212) to find the minimum-cost set of model modifications that make the objective (biomass) >= *min_growth*. Two repair mechanisms are supported: - **Database reactions** (``y_db``): add reactions from a universal template. - **Directionality reversal** (``y_rev``): flip an existing draft reaction so - it can also carry flux in the reverse direction. + it can also carry flux in the reverse direction. Every irreversible draft + reaction is a candidate, exchange reactions and maintenance drains such as + ATP hydrolysis included, and with the default weights a reversal is cheaper + than an addition; check ``reversed_reactions`` before accepting a fill. The universal template should include any exchange or transport reactions that are desired as repair candidates; they are treated identically to other From 0c3090227121d18cf26ab69790d42d79027eeeec Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 09:23:53 +0200 Subject: [PATCH 14/17] fix: stop gap_report on infeasible models with a pointer to forced fluxes gap_report raised a bare cobra Infeasible partway through when the model was infeasible under its own bounds or once its boundary was closed. It now raises a ValueError up front naming the reactions with a forced flux and asking for ATP maintenance and similar bounds to be set to zero. --- src/raven_toolbox/gapfilling/report.py | 22 +++++++++++++++++++--- tests/test_gapfilling_report.py | 15 +++++++++++++++ 2 files changed, 34 insertions(+), 3 deletions(-) diff --git a/src/raven_toolbox/gapfilling/report.py b/src/raven_toolbox/gapfilling/report.py index 9365f96..6677764 100644 --- a/src/raven_toolbox/gapfilling/report.py +++ b/src/raven_toolbox/gapfilling/report.py @@ -29,6 +29,7 @@ """ from __future__ import annotations +import math from collections.abc import Iterable from dataclasses import dataclass from typing import Literal @@ -258,6 +259,18 @@ def _check_production( return not_produced, needed_for_production, min_to_connect +def _require_feasible(model: cobra.Model, condition: str) -> None: + """Raise a ValueError naming the forced-flux reactions if ``model`` is infeasible.""" + if not math.isnan(model.slim_optimize(error_value=float("nan"))): + return + forced = [r.id for r in model.reactions if r.lower_bound > 0 or r.upper_bound < 0] + hint = f" Reactions with a forced flux: {', '.join(forced)}." if forced else "" + raise ValueError( + f"The model is infeasible {condition}, so gap_report cannot analyse it.{hint} " + "Set the bounds of ATP maintenance and any other forced flux to zero first." + ) + + def gap_report( model: cobra.Model, *, @@ -280,12 +293,14 @@ def gap_report( Raises ------ - cobra.exceptions.Infeasible + ValueError If ``model`` is infeasible under its own bounds, or once its boundary reactions are closed for ``can_produce_without_input`` and - ``can_consume_without_output`` (for example when a maintenance - reaction has a positive lower bound). + ``can_consume_without_output``. The usual cause is a reaction with a + forced flux, such as ATP maintenance with a positive lower bound; set + its bounds to zero before calling ``gap_report``. """ + _require_feasible(model, "under its own bounds") no_flux_reactions = [r.id for r in model.reactions if r.id in _blocked_ids(model)] no_flux_reactions_relaxed = [r.id for r in model.reactions if r.id in _blocked_ids(_relaxed_copy(model))] @@ -298,6 +313,7 @@ def gap_report( closed = model.copy() for rxn in closed.boundary: rxn.bounds = (0.0, 0.0) + _require_feasible(closed, "once its boundary reactions are closed") can_produce_without_input = _can_exchange(closed, "produce") can_consume_without_output = _can_exchange(closed, "consume") diff --git a/tests/test_gapfilling_report.py b/tests/test_gapfilling_report.py index a845155..c727c70 100644 --- a/tests/test_gapfilling_report.py +++ b/tests/test_gapfilling_report.py @@ -148,3 +148,18 @@ def test_template_models_connects_a_blocked_reaction(): assert result.gap_fill is not None assert "R3" in result.gap_fill.newly_connected assert "EX_x" in result.gap_fill.added_reactions + + +def test_infeasible_model_raises_with_maintenance_hint(gap_report_model): + atpm = cobra.Reaction("ATPM", lower_bound=5.0, upper_bound=1000.0) + atpm.add_metabolites({gap_report_model.metabolites.get_by_id("z"): -1}) + gap_report_model.add_reactions([atpm]) + with pytest.raises(ValueError, match="infeasible under its own bounds.*ATPM.*ATP maintenance"): + gap_report(gap_report_model) + + +def test_forced_flux_infeasible_once_closed_raises(gap_report_model): + # Feasible with the uptake open, infeasible once every boundary is closed. + gap_report_model.reactions.get_by_id("R1").lower_bound = 1.0 + with pytest.raises(ValueError, match="once its boundary reactions are closed.*R1"): + gap_report(gap_report_model) From 94cf9b7e3c78746883398901c75a92a6b35880ca Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 13:25:26 +0200 Subject: [PATCH 15/17] fill_gaps_kumar_milp adds template reactions only Draft reactions keep their bounds; the MILP minimises the number of added template reactions. reversed_reactions and the weights argument are removed. --- src/raven_toolbox/gapfilling/__init__.py | 8 +- src/raven_toolbox/gapfilling/kumar_milp.py | 96 ++++------------------ tests/test_gapfilling_kumar_milp.py | 71 +++++++++++----- 3 files changed, 72 insertions(+), 103 deletions(-) diff --git a/src/raven_toolbox/gapfilling/__init__.py b/src/raven_toolbox/gapfilling/__init__.py index 913cf3c..ac49ae7 100644 --- a/src/raven_toolbox/gapfilling/__init__.py +++ b/src/raven_toolbox/gapfilling/__init__.py @@ -8,11 +8,9 @@ required; scales better to large models (fastGapFill / swiftGapFill). - :func:`analyse_topology` (no solver) — BFS metabolite-producibility pre-screening; identifies unreachable metabolites and prunes the candidate reaction pool. -- :func:`fill_gaps_kumar_milp` (objective MILP) — growth-floor MILP with directionality - reversal repair (Kumar et al. 2007). - -For the objective-feasibility flavour without directionality repair, also see -``cobra.flux_analysis.gapfill``. +- :func:`fill_gaps_kumar_milp` (objective MILP) — fewest template reactions that let the + objective reach a growth floor (Kumar et al. 2007); ``cobra.flux_analysis.gapfill`` + solves a similar problem. Separately, :func:`find_leak_metabolite` tests whether the model itself contains a stoichiometric leak (a metabolite it can produce or consume for free) rather than diff --git a/src/raven_toolbox/gapfilling/kumar_milp.py b/src/raven_toolbox/gapfilling/kumar_milp.py index 8653d46..3e4477e 100644 --- a/src/raven_toolbox/gapfilling/kumar_milp.py +++ b/src/raven_toolbox/gapfilling/kumar_milp.py @@ -1,16 +1,10 @@ """Objective-based gap-filling via a global growth-floor MILP. -Implements the single-level MILP formulation from Kumar et al. 2007 -(BMC Bioinformatics 8:212, doi:10.1186/1471-2105-8-212) to find the minimum-cost set of model -modifications that make the objective (biomass) >= *min_growth*. - -Two repair mechanisms are supported: -- **Database reactions** (``y_db``): add reactions from a universal template. -- **Directionality reversal** (``y_rev``): flip an existing draft reaction so - it can also carry flux in the reverse direction. Every irreversible draft - reaction is a candidate, exchange reactions and maintenance drains such as - ATP hydrolysis included, and with the default weights a reversal is cheaper - than an addition; check ``reversed_reactions`` before accepting a fill. +Finds the smallest set of template reactions that makes the objective +(biomass) >= *min_growth*, following the database repair of Kumar et al. 2007 +(BMC Bioinformatics 8:212, doi:10.1186/1471-2105-8-212). Reactions already in +the draft keep their bounds: a reaction with the wrong direction is a curation +issue and is not changed here. The universal template should include any exchange or transport reactions that are desired as repair candidates; they are treated identically to other @@ -18,12 +12,9 @@ This is computationally more expensive than LP-based gap-filling (:mod:`~raven_toolbox.gapfilling.fast_lp`) because it involves a MILP solve, -but it supports the directionality-reversal repair mechanism and optimises -for objective feasibility rather than mere connectivity. - -Note: ``cobra.flux_analysis.gapfill`` also solves an objective-based gap-fill -(by default without directionality reversal). Use that function if reversal -is not needed, as it has a more mature solver interface. +but it optimises for objective feasibility rather than mere connectivity. +It is the counterpart of RAVEN's ``gapFillMILP``; ``cobra.flux_analysis.gapfill`` +solves a similar problem. """ from __future__ import annotations @@ -50,12 +41,8 @@ class KumarGapFillResult: ---------- added_reactions: Template reaction IDs selected by the MILP. - reversed_reactions: - Draft reaction IDs whose directionality is reversed (lower_bound - changed to allow negative flux). model: - Draft model with *added_reactions* incorporated and *reversed_reactions* - bounds updated. + Draft model with *added_reactions* incorporated. exit_status: Solver status string (e.g. ``"optimal"``, ``"infeasible"``), or ``"not_restored"`` when the solver reported a solution but the @@ -64,7 +51,6 @@ class KumarGapFillResult: """ added_reactions: list[str] = field(default_factory=list) - reversed_reactions: list[str] = field(default_factory=list) model: cobra.Model = field(default_factory=cobra.Model) exit_status: str = "not_solved" @@ -91,15 +77,14 @@ def fill_gaps_kumar_milp( templates: cobra.Model | Iterable[cobra.Model], *, min_growth: float | None = None, - weights: tuple[float, float] = (1.0, 2.0), big_m: float = 1000.0, verbose: bool = True, ) -> KumarGapFillResult: """Gap-fill to make the model objective (biomass) >= min_growth. - Solves a global MILP that selects the minimum-weight combination of - directionality reversals and template reactions to make the model's - objective reaction feasible at *min_growth*. + Solves a global MILP that selects the fewest template reactions that make + the model's objective reaction feasible at *min_growth*. The draft + reactions keep their bounds. Parameters ---------- @@ -111,9 +96,6 @@ def fill_gaps_kumar_milp( min_growth: Minimum required objective value. If ``None``, runs FBA on the merged model and uses 10% of the maximum as the floor. - weights: - ``(w_rev, w_db)`` — cost per reversed draft reaction and cost per - added template reaction. Lower-weight repairs are preferred. big_m: Big-M constant for coupling constraints. Must exceed any expected flux magnitude in the model. @@ -125,7 +107,6 @@ def fill_gaps_kumar_milp( KumarGapFillResult """ templates = _as_models(templates) - w_rev, w_db = weights # ---- Merge draft + templates ---- working, template_ids = _merge_templates(model, templates) @@ -147,20 +128,11 @@ def fill_gaps_kumar_milp( if min_growth is None: # FBA on merged model (all template reactions available). - # Also try with all draft reversals relaxed, since the model may only - # achieve growth after directionality repair. max_val = working.slim_optimize() - if not math.isfinite(max_val) or max_val <= 0: - # Retry with all irreversible draft reactions temporarily reversed - with working as w: - for rxn in w.reactions: - if rxn.id in {r.id for r in model.reactions} and rxn.lower_bound >= 0: - rxn.lower_bound = -big_m - max_val = w.slim_optimize() if not math.isfinite(max_val) or max_val <= 0: if verbose: print( - f"fill_gaps_kumar_milp: FBA on merged+relaxed model gives {max_val}. " + f"fill_gaps_kumar_milp: FBA on the merged model gives {max_val}. " "Cannot establish a growth floor." ) return KumarGapFillResult(model=model.copy(), exit_status="infeasible") @@ -177,30 +149,6 @@ def fill_gaps_kumar_milp( extra_cons: list = [] obj_terms: list = [] - # Reversal candidates: draft reactions that are currently irreversible (lb >= 0) - rev_indicators: dict[str, Any] = {} - draft_ids = {r.id for r in model.reactions} - - for rxn in list(working.reactions): - if rxn.id not in draft_ids or rxn.lower_bound < 0: - continue - # Binary y_rev: when 1, allows the reaction to carry negative flux - y = prob.Variable(f"_gf_rev_{rxn.id}", type="binary") - indicators.append(y) - rev_indicators[rxn.id] = y - obj_terms.append(w_rev * y) - - # Coupling: v_j + bigM * y_rev >= 0 - # (preserves v_j >= 0 when y_rev = 0; allows v_j >= -bigM when y_rev = 1) - # Relax the reaction's lower bound to allow potential reversal - rxn.lower_bound = -big_m - c = prob.Constraint( - rxn.flux_expression + big_m * y, - lb=0, - name=f"_gf_revc_{rxn.id}", - ) - extra_cons.append(c) - # Database reaction candidates db_indicators: dict[str, Any] = {} finite_bounds = [ @@ -216,7 +164,7 @@ def fill_gaps_kumar_milp( y = prob.Variable(f"_gf_db_{rid}", type="binary") indicators.append(y) db_indicators[rid] = y - obj_terms.append(w_db * y) + obj_terms.append(y) ub = rxn.upper_bound if math.isfinite(rxn.upper_bound) else effective_big_m lb = rxn.lower_bound if math.isfinite(rxn.lower_bound) else -effective_big_m @@ -245,7 +193,7 @@ def fill_gaps_kumar_milp( elif coeff < 0: obj_rxn.upper_bound = min(obj_rxn.upper_bound, -min_growth) - # MILP objective: minimise weighted sum of binary variables + # MILP objective: minimise the number of added template reactions if obj_terms: working.objective = prob.Objective(sum(obj_terms), direction="min") else: @@ -257,7 +205,7 @@ def fill_gaps_kumar_milp( if verbose: print( f"fill_gaps_kumar_milp: solving MILP " - f"({len(rev_indicators)} reversal, {len(db_indicators)} database candidates)..." + f"({len(db_indicators)} database candidates)..." ) # A fractional indicator within the solver's default integrality tolerance @@ -277,27 +225,18 @@ def fill_gaps_kumar_milp( # ---- Extract solution ---- threshold = 0.5 - reversed_rxn_ids = [ - rid for rid, y in rev_indicators.items() if (y.primal or 0) > threshold - ] added_rxn_ids = [ rid for rid, y in db_indicators.items() if (y.primal or 0) > threshold ] if verbose: print( - f"fill_gaps_kumar_milp: reversed {len(reversed_rxn_ids)} reaction(s), " - f"added {len(added_rxn_ids)} template reaction(s)." + f"fill_gaps_kumar_milp: added {len(added_rxn_ids)} template reaction(s)." ) # ---- Build output model ---- filled = model.copy() - # Apply directionality reversals - for rid in reversed_rxn_ids: - if rid in filled.reactions: - filled.reactions.get_by_id(rid).lower_bound = -big_m - # Add selected template reactions remaining = set(added_rxn_ids) for t in templates: @@ -318,7 +257,6 @@ def fill_gaps_kumar_milp( return KumarGapFillResult( added_reactions=sorted(added_rxn_ids), - reversed_reactions=sorted(reversed_rxn_ids), model=filled, exit_status=status, ) diff --git a/tests/test_gapfilling_kumar_milp.py b/tests/test_gapfilling_kumar_milp.py index 23a5533..7e8c001 100644 --- a/tests/test_gapfilling_kumar_milp.py +++ b/tests/test_gapfilling_kumar_milp.py @@ -52,8 +52,11 @@ def directionality_gap(): draft.add_reactions([exb, r1, bio]) draft.objective = "biomass" - # Empty template — reversal is the only repair - template = cobra.Model("template_empty") + # The template offers nothing that makes A + template = cobra.Model("template") + other = cobra.Reaction("r_other", lower_bound=0, upper_bound=1000) + other.add_metabolites({_met("B_c"): -1}) + template.add_reactions([other]) return draft, template @@ -69,12 +72,6 @@ def test_added_reactions_list(growth_gap): assert isinstance(result.added_reactions, list) -def test_reversed_reactions_list(growth_gap): - draft, template = growth_gap - result = fill_gaps_kumar_milp(draft, template, verbose=False) - assert isinstance(result.reversed_reactions, list) - - def test_model_returned(growth_gap): draft, template = growth_gap result = fill_gaps_kumar_milp(draft, template, verbose=False) @@ -100,13 +97,40 @@ def test_r_missing_added(growth_gap): assert "r_missing" in result.added_reactions -def test_directionality_reversal(directionality_gap): - """r1 should be reversed when that is the only repair option.""" +def test_draft_directions_are_kept(directionality_gap): + """A draft reaction in the wrong direction is not reversed: there is no fill.""" draft, template = directionality_gap - result = fill_gaps_kumar_milp(draft, template, verbose=False) - if result.exit_status != "optimal": - pytest.skip(f"Solver status: {result.exit_status}") - assert "r1" in result.reversed_reactions + result = fill_gaps_kumar_milp(draft, template, min_growth=1.0, verbose=False) + assert result.exit_status == "infeasible" + assert result.added_reactions == [] + assert result.model.reactions.r1.bounds == (0, 1000) + + +def test_forced_flux_is_kept(): + """A draft reaction with a forced flux (maintenance) keeps its lower bound.""" + A, B = _met("A_c"), _met("B_c") + draft = cobra.Model("draft") + exa = cobra.Reaction("EX_A", lower_bound=0, upper_bound=10) + exa.add_metabolites({A: 1}) + r1 = cobra.Reaction("r1", lower_bound=0, upper_bound=1000) + r1.add_metabolites({A: -1, B: 1}) + grow = cobra.Reaction("grow", lower_bound=0, upper_bound=1000) + grow.add_metabolites({B: -1}) + maint = cobra.Reaction("maint", lower_bound=5, upper_bound=1000) + maint.add_metabolites({A: -1}) + draft.add_reactions([exa, r1, grow, maint]) + draft.objective = "grow" + + template = cobra.Model("template") + u1 = cobra.Reaction("u1", lower_bound=0, upper_bound=10) + u1.add_metabolites({_met("A_c"): 1}) + template.add_reactions([u1]) + + result = fill_gaps_kumar_milp(draft, template, min_growth=8.0, verbose=False) + assert result.exit_status == "optimal" + assert result.added_reactions == ["u1"] + assert result.model.reactions.maint.lower_bound == 5 + assert result.model.slim_optimize() >= 8 - 1e-6 def test_no_objective_raises(growth_gap): @@ -152,14 +176,23 @@ def test_small_floor_is_not_met_by_a_fractional_indicator(growth_gap, solver): assert result.model.slim_optimize() >= 0.005 +def test_floor_check_on_filled_model(growth_gap): + """The floor check fails on the draft and passes once r_missing is added.""" + from raven_toolbox.gapfilling import kumar_milp + + draft, template = growth_gap + obj = {draft.reactions.biomass: 1.0} + assert not kumar_milp._meets_floor(draft, obj, 0.005) + draft.add_reactions([template.reactions.r_missing.copy()]) + assert kumar_milp._meets_floor(draft, obj, 0.005) + + def test_selection_that_does_not_restore_growth_is_reported(growth_gap, monkeypatch): - """A rounded MILP selection is checked on the filled model.""" + """A selection that fails the floor check is returned as 'not_restored'.""" from raven_toolbox.gapfilling import kumar_milp draft, template = growth_gap - draft.solver = "glpk" - # A loose integrality tolerance lets the indicator stay fractional. - monkeypatch.setattr(kumar_milp, "_INTEGRALITY_TOL", 1e-4) + monkeypatch.setattr(kumar_milp, "_meets_floor", lambda *args: False) result = fill_gaps_kumar_milp(draft, template, min_growth=0.005, verbose=False) - assert result.added_reactions == [] + assert result.added_reactions == ["r_missing"] assert result.exit_status == "not_restored" From 1603f022ceb47d8efea79fdad33b7424250470a3 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 13:34:30 +0200 Subject: [PATCH 16/17] fill_gaps_kumar_milp: report an infeasible MILP on Gurobi instead of crashing --- src/raven_toolbox/gapfilling/kumar_milp.py | 5 +++-- tests/test_gapfilling_kumar_milp.py | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/raven_toolbox/gapfilling/kumar_milp.py b/src/raven_toolbox/gapfilling/kumar_milp.py index 3e4477e..905661c 100644 --- a/src/raven_toolbox/gapfilling/kumar_milp.py +++ b/src/raven_toolbox/gapfilling/kumar_milp.py @@ -215,8 +215,9 @@ def fill_gaps_kumar_milp( with contextlib.suppress(Exception): working.solver.configuration.tolerances.integrality = _INTEGRALITY_TOL - working.optimize() - status = working.solver.status + # Solve through optlang: cobra's optimize() reads primal values even when + # there are none, which Gurobi refuses for an infeasible MILP. + status = working.solver.optimize() if status != "optimal": if verbose: diff --git a/tests/test_gapfilling_kumar_milp.py b/tests/test_gapfilling_kumar_milp.py index 7e8c001..4221cff 100644 --- a/tests/test_gapfilling_kumar_milp.py +++ b/tests/test_gapfilling_kumar_milp.py @@ -101,7 +101,7 @@ def test_draft_directions_are_kept(directionality_gap): """A draft reaction in the wrong direction is not reversed: there is no fill.""" draft, template = directionality_gap result = fill_gaps_kumar_milp(draft, template, min_growth=1.0, verbose=False) - assert result.exit_status == "infeasible" + assert result.exit_status != "optimal" assert result.added_reactions == [] assert result.model.reactions.r1.bounds == (0, 1000) From cd65beece509d2b7e443ec8dfc176cf0a1fefafa Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Sat, 17 Oct 2026 21:27:10 +0200 Subject: [PATCH 17/17] Rename fill_gaps_kumar_milp to fill_gaps_milp Without directionality reversal the function is no longer Kumar's GapFill; the name now matches RAVEN's gapFillMILP. KumarGapFillResult becomes MILPGapFillResult and the module gapfilling.kumar_milp becomes gapfilling.milp. --- src/raven_toolbox/gapfilling/__init__.py | 8 ++-- .../gapfilling/{kumar_milp.py => milp.py} | 36 +++++++++--------- ..._kumar_milp.py => test_gapfilling_milp.py} | 38 +++++++++---------- 3 files changed, 41 insertions(+), 41 deletions(-) rename src/raven_toolbox/gapfilling/{kumar_milp.py => milp.py} (89%) rename tests/{test_gapfilling_kumar_milp.py => test_gapfilling_milp.py} (82%) diff --git a/src/raven_toolbox/gapfilling/__init__.py b/src/raven_toolbox/gapfilling/__init__.py index ac49ae7..0230114 100644 --- a/src/raven_toolbox/gapfilling/__init__.py +++ b/src/raven_toolbox/gapfilling/__init__.py @@ -8,7 +8,7 @@ required; scales better to large models (fastGapFill / swiftGapFill). - :func:`analyse_topology` (no solver) — BFS metabolite-producibility pre-screening; identifies unreachable metabolites and prunes the candidate reaction pool. -- :func:`fill_gaps_kumar_milp` (objective MILP) — fewest template reactions that let the +- :func:`fill_gaps_milp` (objective MILP) — fewest template reactions that let the objective reach a growth floor (Kumar et al. 2007); ``cobra.flux_analysis.gapfill`` solves a similar problem. @@ -19,8 +19,8 @@ """ from raven_toolbox.gapfilling.fast_lp import FastLPResult, fill_gaps_fast_lp from raven_toolbox.gapfilling.fill import GapFillResult, connect_blocked_reactions -from raven_toolbox.gapfilling.kumar_milp import KumarGapFillResult, fill_gaps_kumar_milp from raven_toolbox.gapfilling.leak import LeakMetaboliteResult, find_leak_metabolite +from raven_toolbox.gapfilling.milp import MILPGapFillResult, fill_gaps_milp from raven_toolbox.gapfilling.report import ( GapReportResult, MinToConnect, @@ -36,8 +36,8 @@ "fill_gaps_fast_lp", "TopologicalAnalysisResult", "analyse_topology", - "KumarGapFillResult", - "fill_gaps_kumar_milp", + "MILPGapFillResult", + "fill_gaps_milp", "LeakMetaboliteResult", "find_leak_metabolite", "GapReportResult", diff --git a/src/raven_toolbox/gapfilling/kumar_milp.py b/src/raven_toolbox/gapfilling/milp.py similarity index 89% rename from src/raven_toolbox/gapfilling/kumar_milp.py rename to src/raven_toolbox/gapfilling/milp.py index 905661c..aa11f4e 100644 --- a/src/raven_toolbox/gapfilling/kumar_milp.py +++ b/src/raven_toolbox/gapfilling/milp.py @@ -34,8 +34,8 @@ @dataclass -class KumarGapFillResult: - """Outcome of a Kumar 2007 MILP gap-fill. +class MILPGapFillResult: + """Outcome of a growth-floor MILP gap-fill. Parameters ---------- @@ -72,14 +72,14 @@ def _merge_templates( return working, template_ids -def fill_gaps_kumar_milp( +def fill_gaps_milp( model: cobra.Model, templates: cobra.Model | Iterable[cobra.Model], *, min_growth: float | None = None, big_m: float = 1000.0, verbose: bool = True, -) -> KumarGapFillResult: +) -> MILPGapFillResult: """Gap-fill to make the model objective (biomass) >= min_growth. Solves a global MILP that selects the fewest template reactions that make @@ -104,7 +104,7 @@ def fill_gaps_kumar_milp( Returns ------- - KumarGapFillResult + MILPGapFillResult """ templates = _as_models(templates) @@ -113,7 +113,7 @@ def fill_gaps_kumar_milp( if verbose: print( - f"fill_gaps_kumar_milp: merged model has {len(working.reactions)} reactions " + f"fill_gaps_milp: merged model has {len(working.reactions)} reactions " f"({len(model.reactions)} draft, {len(template_ids)} template)." ) @@ -122,7 +122,7 @@ def fill_gaps_kumar_milp( obj_coeffs = linear_reaction_coefficients(working) if not obj_coeffs: raise ValueError( - "fill_gaps_kumar_milp: model has no objective function. " + "fill_gaps_milp: model has no objective function. " "Set a non-zero objective before calling this function." ) @@ -132,14 +132,14 @@ def fill_gaps_kumar_milp( if not math.isfinite(max_val) or max_val <= 0: if verbose: print( - f"fill_gaps_kumar_milp: FBA on the merged model gives {max_val}. " + f"fill_gaps_milp: FBA on the merged model gives {max_val}. " "Cannot establish a growth floor." ) - return KumarGapFillResult(model=model.copy(), exit_status="infeasible") + return MILPGapFillResult(model=model.copy(), exit_status="infeasible") min_growth = 0.1 * max_val if verbose: print( - f"fill_gaps_kumar_milp: setting min_growth = {min_growth:.4g} " + f"fill_gaps_milp: setting min_growth = {min_growth:.4g} " f"(10% of FBA max {max_val:.4g})." ) @@ -198,13 +198,13 @@ def fill_gaps_kumar_milp( working.objective = prob.Objective(sum(obj_terms), direction="min") else: if verbose: - print("fill_gaps_kumar_milp: no repair candidates found.") - return KumarGapFillResult(model=model.copy(), exit_status="infeasible") + print("fill_gaps_milp: no repair candidates found.") + return MILPGapFillResult(model=model.copy(), exit_status="infeasible") # ---- Solve ---- if verbose: print( - f"fill_gaps_kumar_milp: solving MILP " + f"fill_gaps_milp: solving MILP " f"({len(db_indicators)} database candidates)..." ) @@ -221,8 +221,8 @@ def fill_gaps_kumar_milp( if status != "optimal": if verbose: - print(f"fill_gaps_kumar_milp: solver returned '{status}'.") - return KumarGapFillResult(model=model.copy(), exit_status=status) + print(f"fill_gaps_milp: solver returned '{status}'.") + return MILPGapFillResult(model=model.copy(), exit_status=status) # ---- Extract solution ---- threshold = 0.5 @@ -232,7 +232,7 @@ def fill_gaps_kumar_milp( if verbose: print( - f"fill_gaps_kumar_milp: added {len(added_rxn_ids)} template reaction(s)." + f"fill_gaps_milp: added {len(added_rxn_ids)} template reaction(s)." ) # ---- Build output model ---- @@ -251,12 +251,12 @@ def fill_gaps_kumar_milp( if not _meets_floor(filled, obj_coeffs, min_growth): if verbose: print( - "fill_gaps_kumar_milp: the selected repairs do not reach " + "fill_gaps_milp: the selected repairs do not reach " f"min_growth {min_growth:.4g} in the filled model." ) status = "not_restored" - return KumarGapFillResult( + return MILPGapFillResult( added_reactions=sorted(added_rxn_ids), model=filled, exit_status=status, diff --git a/tests/test_gapfilling_kumar_milp.py b/tests/test_gapfilling_milp.py similarity index 82% rename from tests/test_gapfilling_kumar_milp.py rename to tests/test_gapfilling_milp.py index 4221cff..cf26b66 100644 --- a/tests/test_gapfilling_kumar_milp.py +++ b/tests/test_gapfilling_milp.py @@ -1,8 +1,8 @@ -"""Tests for growth-floor MILP gap-filling (gapfilling/kumar_milp.py).""" +"""Tests for growth-floor MILP gap-filling (gapfilling/milp.py).""" import cobra import pytest -from raven_toolbox.gapfilling import KumarGapFillResult, fill_gaps_kumar_milp +from raven_toolbox.gapfilling import MILPGapFillResult, fill_gaps_milp def _met(mid): @@ -62,26 +62,26 @@ def directionality_gap(): def test_result_type(growth_gap): draft, template = growth_gap - result = fill_gaps_kumar_milp(draft, template, verbose=False) - assert isinstance(result, KumarGapFillResult) + result = fill_gaps_milp(draft, template, verbose=False) + assert isinstance(result, MILPGapFillResult) def test_added_reactions_list(growth_gap): draft, template = growth_gap - result = fill_gaps_kumar_milp(draft, template, verbose=False) + result = fill_gaps_milp(draft, template, verbose=False) assert isinstance(result.added_reactions, list) def test_model_returned(growth_gap): draft, template = growth_gap - result = fill_gaps_kumar_milp(draft, template, verbose=False) + result = fill_gaps_milp(draft, template, verbose=False) assert isinstance(result.model, cobra.Model) def test_gap_filled_enables_growth(growth_gap): """After adding r_missing, the model should be able to produce biomass.""" draft, template = growth_gap - result = fill_gaps_kumar_milp(draft, template, verbose=False) + result = fill_gaps_milp(draft, template, verbose=False) if result.exit_status != "optimal": pytest.skip(f"Solver status: {result.exit_status}") val = result.model.slim_optimize() @@ -91,7 +91,7 @@ def test_gap_filled_enables_growth(growth_gap): def test_r_missing_added(growth_gap): """r_missing should be in the returned model.""" draft, template = growth_gap - result = fill_gaps_kumar_milp(draft, template, verbose=False) + result = fill_gaps_milp(draft, template, verbose=False) if result.exit_status != "optimal": pytest.skip(f"Solver status: {result.exit_status}") assert "r_missing" in result.added_reactions @@ -100,7 +100,7 @@ def test_r_missing_added(growth_gap): def test_draft_directions_are_kept(directionality_gap): """A draft reaction in the wrong direction is not reversed: there is no fill.""" draft, template = directionality_gap - result = fill_gaps_kumar_milp(draft, template, min_growth=1.0, verbose=False) + result = fill_gaps_milp(draft, template, min_growth=1.0, verbose=False) assert result.exit_status != "optimal" assert result.added_reactions == [] assert result.model.reactions.r1.bounds == (0, 1000) @@ -126,7 +126,7 @@ def test_forced_flux_is_kept(): u1.add_metabolites({_met("A_c"): 1}) template.add_reactions([u1]) - result = fill_gaps_kumar_milp(draft, template, min_growth=8.0, verbose=False) + result = fill_gaps_milp(draft, template, min_growth=8.0, verbose=False) assert result.exit_status == "optimal" assert result.added_reactions == ["u1"] assert result.model.reactions.maint.lower_bound == 5 @@ -138,7 +138,7 @@ def test_no_objective_raises(growth_gap): draft, template = growth_gap draft.objective = {} # clear objective with pytest.raises(ValueError, match="no objective"): - fill_gaps_kumar_milp(draft, template, verbose=False) + fill_gaps_milp(draft, template, verbose=False) def test_exit_status_on_infeasible(): @@ -151,7 +151,7 @@ def test_exit_status_on_infeasible(): m.add_reactions([bio]) m.objective = "bio" template = cobra.Model("empty_template") - result = fill_gaps_kumar_milp(m, template, verbose=False) + result = fill_gaps_milp(m, template, verbose=False) assert result.exit_status != "optimal" @@ -170,7 +170,7 @@ def test_small_floor_is_not_met_by_a_fractional_indicator(growth_gap, solver): """With a floor below 1e-5 x big-M, a near-zero indicator must not pass as 'not added'.""" draft, template = growth_gap draft.solver = solver - result = fill_gaps_kumar_milp(draft, template, min_growth=0.005, verbose=False) + result = fill_gaps_milp(draft, template, min_growth=0.005, verbose=False) assert result.exit_status == "optimal" assert result.added_reactions == ["r_missing"] assert result.model.slim_optimize() >= 0.005 @@ -178,21 +178,21 @@ def test_small_floor_is_not_met_by_a_fractional_indicator(growth_gap, solver): def test_floor_check_on_filled_model(growth_gap): """The floor check fails on the draft and passes once r_missing is added.""" - from raven_toolbox.gapfilling import kumar_milp + from raven_toolbox.gapfilling import milp draft, template = growth_gap obj = {draft.reactions.biomass: 1.0} - assert not kumar_milp._meets_floor(draft, obj, 0.005) + assert not milp._meets_floor(draft, obj, 0.005) draft.add_reactions([template.reactions.r_missing.copy()]) - assert kumar_milp._meets_floor(draft, obj, 0.005) + assert milp._meets_floor(draft, obj, 0.005) def test_selection_that_does_not_restore_growth_is_reported(growth_gap, monkeypatch): """A selection that fails the floor check is returned as 'not_restored'.""" - from raven_toolbox.gapfilling import kumar_milp + from raven_toolbox.gapfilling import milp draft, template = growth_gap - monkeypatch.setattr(kumar_milp, "_meets_floor", lambda *args: False) - result = fill_gaps_kumar_milp(draft, template, min_growth=0.005, verbose=False) + monkeypatch.setattr(milp, "_meets_floor", lambda *args: False) + result = fill_gaps_milp(draft, template, min_growth=0.005, verbose=False) assert result.added_reactions == ["r_missing"] assert result.exit_status == "not_restored"