From 37a67876ac75ec4c7035038e04ea4999d86b0454 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Mon, 20 Jul 2026 10:31:29 +0200 Subject: [PATCH 1/2] feat(ftinit): opt-in deterministic extraction (strict gap + canonical optimum) The ftINIT MILP is degenerate: many reaction subsets reach the same score optimum and the solver returns an arbitrary one. This is reproducible for a fixed solver build + seed but fragile to the Gurobi version or platform, so a downstream metric such as gene essentiality can drift on solver noise rather than on the model. Measured on Human-GEM/DLD1: changing only the solver seed moves ~125 of 7766 kept reactions and flips 5 essential-gene calls. Add two opt-in flags (default off, exact RAVEN behaviour preserved): - strict_gap: one near-proven-optimal solve per step instead of RAVEN's loose absolute-gap escalation, whose final near-zero-objective run accepts an arbitrary within-gap incumbent. - canonical: a lexicographic phase 2 that, holding the score objective at its optimum, minimises the id-ordered count of kept reactions, selecting the unique sparsest optimum instead of an arbitrary tie-break. Threaded through run_ftinit / _solve_step / ftinit. Tests cover degenerate-tie resolution, no regression on the unique-optimum oracle, and staged stability. --- src/raven_toolbox/init/ftinit.py | 166 ++++++++++++++++++++++++++--- src/raven_toolbox/init/taskfill.py | 66 +++++++++++- tests/test_init_ftinit.py | 42 ++++++++ tests/test_init_pipeline.py | 23 ++++ tests/test_init_taskfill.py | 23 ++++ 5 files changed, 298 insertions(+), 22 deletions(-) diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index 5d45b5e..62acd52 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -62,6 +62,7 @@ _FORCE_ON = 0.1 # min flux for a reaction to count as "on" (RAVEN forceOnLim) _BIG_M = 100.0 # indicator/direction big-M cap on a *scored* reaction's flux (RAVEN's 100) +_STRICT_ABS_GAP = 0.05 # absolute gap for the opt-in strict mode (below the 0.1 score granularity) def _dbg(msg: str) -> None: @@ -101,6 +102,8 @@ def run_ftinit( mip_gap: float | None = None, mip_gap_abs: float | None = None, time_limit: float | None = None, + strict_abs_gap: float | None = None, + canonical: bool = False, ) -> FtInitResult: """Run the single-step ftINIT MILP and return the extracted model. @@ -112,6 +115,12 @@ def run_ftinit( per-step "simple metabolite" removal, e.g. H2O/H+). See the module docstring for the formulation. This is the single-step variant; the staged schedule (:func:`raven_toolbox.init.ftinit`) calls it per step. + + ``canonical`` (opt-in, default off to preserve RAVEN parity) resolves the MILP's + degeneracy deterministically: after the score optimum is found, a lexicographic + phase 2 holds the objective and minimises the id-ordered count of "on" reactions, so + the kept set is the unique sparsest optimum, independent of solver seed/version, + instead of an arbitrary tie-break. See :func:`_canonicalize`. """ scores = dict(rxn_scores or {}) essential = set(essential_rxns or []) @@ -217,9 +226,8 @@ def add_constraint(expr, **kw): add_constraint(add(termlist), lb=0.0, ub=None if allow_excretion else 0.0) opt.add(variables + constraints) - opt.objective = prob.Objective( - add([mul([Real(score), ind]) for ind, score in indicators.values()]), direction="max" - ) + obj_expr = add([mul([Real(score), ind]) for ind, score in indicators.values()]) + opt.objective = prob.Objective(obj_expr, direction="max") try: # Gurobi-specific; harmless if the backend differs. Match RAVEN's optimizeProb # defaults exactly, because the ftINIT MILP is highly degenerate and the chosen # incumbent (hence which reactions are kept) depends on these: @@ -238,13 +246,23 @@ def add_constraint(expr, **kw): opt.problem.Params.Seed = 1234 except Exception: # noqa: BLE001 pass - if mip_gap is not None: + if strict_abs_gap is not None: + # Strict mode: prove the optimum to a fixed *absolute* gap below the reaction-score + # granularity (scores are nudged to |score| ≥ 0.1), so the kept-set objective is + # proven optimal at any objective scale. A relative gap is meaningless on the + # near-zero-objective final step; this also bypasses the loose escalation below. + try: # Gurobi-specific; harmless if the backend differs + opt.problem.Params.MIPGap = 0.0 + opt.problem.Params.MIPGapAbs = strict_abs_gap + except Exception: # noqa: BLE001 + pass + elif mip_gap is not None: try: # Gurobi-specific; harmless if the backend differs opt.problem.Params.MIPGap = mip_gap except Exception: # noqa: BLE001 pass - if mip_gap_abs is not None: + if strict_abs_gap is None and mip_gap_abs is not None: # RAVEN's multi-run gap strategy (ftINIT.m). The final staged step has a # near-zero objective (it mostly removes small negative-score reactions), so a # fixed *relative* gap becomes an almost-zero absolute gap the solver cannot @@ -276,27 +294,115 @@ def add_constraint(expr, **kw): _achieved = None _dbg(f"[ftinit] final solve: obj={opt.objective.value} status={opt.status} " f"achieved_gap={_achieved}") - # Accept a near-optimal incumbent (when a MIP gap / time limit is set), as RAVEN does. - if opt.status not in ("optimal", "feasible", "suboptimal", "time_limit"): - raise OptimizationError(f"ftINIT MILP did not solve (status: {opt.status}).") + # Accept a near-optimal incumbent (when a MIP gap / time limit is set), as RAVEN does, + # but only if the solver actually holds one to read. + if opt.status not in ("optimal", "feasible", "suboptimal", "time_limit") \ + or not _has_solution(opt): + raise OptimizationError( + f"ftINIT MILP produced no usable solution (status: {opt.status}); " + "increase time_limit or disable strict_gap." + ) + + # Report the *primary* (score) objective; a canonical phase 2 replaces the objective + # in place, so capture it before that. + primary_obj = float(opt.objective.value) if opt.objective.value is not None else 0.0 # RAVEN: a reaction is "on" iff its indicator ≥ 0.5 (positive indicators are # continuous and can land fractionally when a reaction can carry only tiny flux). - on = {rid for rid, (ind, _) in indicators.items() if (ind.primal or 0.0) >= 0.5} + def _read_solution(): + on = {rid for rid, (ind, _) in indicators.items() if (ind.primal or 0.0) >= 0.5} + fluxes = {rid: sum(sign * (var.primal or 0.0) for var, sign in terms) + for rid, terms in flux_terms.items()} + return on, fluxes + + on, fluxes = _read_solution() # the primary optimum + # canonical is best-effort: keep its result only if phase 2 actually converged. + if canonical and indicators and _canonicalize(opt, prob, obj_expr, indicators, + primary_obj, time_limit): + on, fluxes = _read_solution() + kept = free_or_essential | on deleted = [r.id for r in model.reactions if r.id not in kept] - fluxes: dict[str, float] = { - rid: sum(sign * (var.primal or 0.0) for var, sign in terms) - for rid, terms in flux_terms.items() - } out = model.copy() out.remove_reactions(deleted, remove_orphans=True) return FtInitResult(out, sorted(kept), sorted(deleted), fluxes, - float(opt.objective.value), on_reactions=on, + primary_obj, on_reactions=on, achieved_gap=_achieved) +def _has_solution(opt) -> bool: + """Whether the solver currently holds a readable primal solution. + + A MILP can finish with an accepted status (notably ``time_limit``) yet no incumbent, + so reading ``.primal`` would raise. On Gurobi this is the solution count; on other + backends the status is taken as authoritative. + """ + try: + return opt.problem.SolCount > 0 + except Exception: # noqa: BLE001 - non-Gurobi backend + return True + + +def _canonicalize(opt, prob, obj_expr, indicators, primary, time_limit) -> bool: + """Pin ftINIT's degenerate optimum to a single canonical solution (in place on ``opt``). + + The ftINIT MILP is highly degenerate: many reaction subsets reach the same score + optimum and the solver returns an arbitrary one — reproducible for a fixed solver + build + seed, but fragile to the Gurobi version or platform (a changed tie-break + moves ~1-2% of the kept reactions, which flips a handful of gene-essentiality calls). + + This runs a lexicographic phase 2, holding the score objective at its optimum with a + floor constraint: first minimise the count of kept removable reactions (the sparsest + optimum), then, among the sparsest, minimise their summed id rank (prefer lower ids). + The result is a stable, near-unique optimum independent of seed/solver version. ``opt`` + is left holding the phase-2 solution, which the caller reads for the "on" set and + fluxes; the reported objective stays the phase-1 ``primary`` value. + + Only the negative-score reactions carry a true 0/1 "keep" binary (the positive + indicators are continuous and pinned near 1 by the score objective), so the two phases + run over those: their count and id-sum are integers, provable with a cheap absolute gap + below 1 rather than the near-full proof a tiny relative gap would need. This + canonicalises the removable-reaction choices — the genuine seed-fragile degeneracy; a + residual flux-distribution degeneracy (which only feeds the next step's small + ``ess_force`` clamp) is left to the solver. + + Returns ``True`` if phase 2 produced a usable solution (the caller then reads the + canonical "on" set and fluxes), ``False`` if it did not converge to an incumbent (the + caller keeps the phase-1 optimum). Canonicalisation is best-effort: it never fails the + extraction. + """ + binaries = {rid: ind for rid, (ind, score) in indicators.items() if score < 0} + if not binaries: # only positive / free reactions: nothing removable to canonicalise + return _has_solution(opt) + tol = max(abs(primary) * 1e-7, 1e-7) + opt.add(prob.Constraint(obj_expr, lb=primary - tol, name="_canon_obj_floor")) + count = add([mul([Real(1.0), ind]) for ind in binaries.values()]) + + def _phase(objective) -> bool: + opt.objective = objective + try: # integer objective: an absolute gap < 1 proves the optimum cheaply. + opt.problem.Params.MIPGap = 0.0 + opt.problem.Params.MIPGapAbs = 0.4 + except Exception: # noqa: BLE001 - GLPK solves exactly; harmless + pass + if time_limit is not None: + opt.configuration.timeout = int(time_limit) + opt.optimize() + return (opt.status in ("optimal", "feasible", "suboptimal", "time_limit") + and _has_solution(opt)) + + # Phase 2a — parsimony: the fewest kept removable reactions. + if not _phase(prob.Objective(count, direction="min")): + return False + kmin = opt.objective.value or 0.0 + # Phase 2b — among the sparsest, prefer lower reaction ids (deterministic tie-break). + opt.add(prob.Constraint(count, ub=kmin + 0.5, name="_canon_count_cap")) + ranks = {rid: i for i, rid in enumerate(sorted(binaries))} + idsum = add([mul([Real(float(1 + ranks[rid])), ind]) for rid, ind in binaries.items()]) + return _phase(prob.Objective(idsum, direction="min")) + + def _nudge_scores(rxn_scores: Mapping[str, float]) -> dict[str, float]: """Push tiny reaction scores off zero (RAVEN ``ftINIT.m:160-161``). @@ -317,7 +423,7 @@ def _nudge_scores(rxn_scores: Mapping[str, float]) -> dict[str, float]: def _solve_step( min_model, scores, step, *, essential, directions, ess_force, force_on, big_m, - mip_gap, mip_gap_abs, time_limit, + mip_gap, mip_gap_abs, time_limit, strict_gap=False, canonical=False, ) -> FtInitResult: """Solve one ftINIT step, following RAVEN's multi-run gap-escalation schedule. @@ -328,15 +434,23 @@ def _solve_step( caller ``time_limit`` caps each run's own limit. With no schedule (``step.milp_runs`` empty, e.g. the ``'full'`` series) this is a single solve at the caller's gap. """ - def _run(mg, mga, tl): + def _run(mg, mga, tl, strict=None): return run_ftinit( min_model, scores, essential_rxns=essential, essential_directions=directions, essential_force=ess_force, allow_excretion=step.allow_met_secr, rem_pos_rev=step.pos_rev_off, ignore_mets=step.mets_to_ignore, force_on=force_on, force_on_ess=force_on, big_m=big_m, mip_gap=mg, mip_gap_abs=mga, time_limit=tl, + strict_abs_gap=strict, canonical=canonical, ) + # Strict mode: one solve per step proven to a fixed absolute gap (below the 0.1 score + # granularity), bypassing RAVEN's loose relative escalation — whose final + # near-zero-objective run otherwise accepts an arbitrary within-gap incumbent. Trades + # runtime for a stable, well-defined optimum (pairs naturally with ``canonical``). + if strict_gap: + return _run(None, None, time_limit, strict=_STRICT_ABS_GAP) + if not step.milp_runs: return _run(mip_gap, mip_gap_abs, time_limit) @@ -372,6 +486,8 @@ def ftinit( mip_gap: float | None = None, mip_gap_abs: float | None = 10.0, time_limit: float | None = None, + strict_gap: bool = False, + canonical: bool = False, ) -> cobra.Model: """Run the full ftINIT pipeline on prepData and return the context-specific model. @@ -402,6 +518,20 @@ def ftinit( ``mip_gap``/``time_limit`` are forwarded to each :func:`run_ftinit` solve. On genome-scale models they are essential for tractability — see ``docs/init_param_calibration.md`` for the calibration table. + + ``strict_gap`` and ``canonical`` (both opt-in, default off → exact RAVEN behaviour) + make the extracted model reproducible across solver versions/platforms instead of + merely across identical runs, which matters when a downstream metric (e.g. gene + essentiality) must not drift on solver noise: + + * ``strict_gap`` replaces the loose relative-gap escalation with a single solve per + step proven to a fixed *absolute* gap (below the 0.1 reaction-score granularity), + so the kept-set objective is optimal at any scale — not an arbitrary within-gap + incumbent. Slower, but removes the largest source of tie-break drift. + * ``canonical`` adds a lexicographic phase 2 that selects the unique sparsest (then + lowest-id) optimum, so the degenerate choice is pinned rather than left to the + solver — applied both to each extraction step and to the task gap-fill. Best used + together with ``strict_gap`` (a well-defined primary optimum to canonicalise). """ if metabolomics: raise NotImplementedError( @@ -436,6 +566,7 @@ def ftinit( min_model, scores, step, essential=essential, directions=directions, ess_force=ess_force, force_on=force_on, big_m=big_m, mip_gap=mip_gap, mip_gap_abs=mip_gap_abs, time_limit=time_limit, + strict_gap=strict_gap, canonical=canonical, ) for rid in res.on_reactions: turned_on[rid] = res.fluxes[rid] @@ -460,7 +591,8 @@ def ftinit( if fill_gaps and prep.tasks: # add reactions back so every task is feasible # The gap-fill MILP is its own problem (RAVEN ftINITFillGaps); it uses RAVEN's # fixed per-task 300 s limit and seed, not the main extraction's time_limit. - out = fill_tasks(out, prep.ref_model, prep.tasks, rxn_scores=rxn_scores).model + out = fill_tasks(out, prep.ref_model, prep.tasks, rxn_scores=rxn_scores, + canonical=canonical).model if gene_scores is not None: # prune negative-scoring genes from the GPRs out, _ = remove_low_score_genes(out, gene_scores) return out diff --git a/src/raven_toolbox/init/taskfill.py b/src/raven_toolbox/init/taskfill.py index ba9bab2..f40b4ee 100644 --- a/src/raven_toolbox/init/taskfill.py +++ b/src/raven_toolbox/init/taskfill.py @@ -125,9 +125,53 @@ def _set_fill_solver(model: cobra.Model, time_limit: float | None, seed: int) -> model.solver.configuration.timeout = int(time_limit) +def _canonicalize_fill(work, prob, candidates, cost_expr, time_limit) -> list[str] | None: + """Pin the degenerate min-cost gap-fill to a single canonical set (in place on ``work``). + + Like the extraction MILP, the fill has many equal-cost solutions and the solver returns + an arbitrary one (seed/version dependent). Hold the cost at its optimum, then lexicographically + minimise the number of added reactions and then their summed id rank — both integer + objectives, provable with a cheap absolute gap below 1. Returns the chosen reaction ids, + or ``None`` if a phase did not converge (the caller then keeps the arbitrary min-cost fill). + """ + yvars = {cid: work.variables[f"_fill_{cid}"] for cid in candidates} + primary_cost = work.objective.value or 0.0 + tol = max(abs(primary_cost) * 1e-7, 1e-7) + work.add_cons_vars([prob.Constraint(cost_expr, ub=primary_cost + tol, name="_fill_cost_floor")]) + count = add([mul([Real(1.0), y]) for y in yvars.values()]) + + def _phase(objective) -> bool: + work.objective = objective + try: # integer objective → an absolute gap < 1 proves the optimum cheaply. + work.solver.problem.Params.MIPGap = 0.0 + work.solver.problem.Params.MIPGapAbs = 0.4 + except Exception: # noqa: BLE001 - harmless on other backends + pass + work.slim_optimize() + if work.solver.status not in ("optimal", "feasible", "suboptimal", "time_limit"): + return False + try: + return work.solver.problem.SolCount > 0 + except Exception: # noqa: BLE001 - non-Gurobi backend: status is authoritative + return True + + # fewest added reactions (parsimony) ... + if not _phase(prob.Objective(count, direction="min")): + return None + kmin = work.objective.value or 0.0 + # ... then, among the sparsest, the unique lowest-id set. + work.add_cons_vars([prob.Constraint(count, ub=kmin + 0.5, name="_fill_count_cap")]) + ranks = {cid: i for i, cid in enumerate(sorted(candidates))} + idsum = add([mul([Real(float(1 + ranks[cid])), yvars[cid]]) for cid in candidates]) + if not _phase(prob.Objective(idsum, direction="min")): + return None + return [cid for cid in candidates if (yvars[cid].primal or 0.0) > 0.5] + + def _gap_fill_task( reference_model: cobra.Model, present_ids: set[str], task: Task, costs: dict[str, float], *, time_limit: float | None, seed: int, + canonical: bool = False, ) -> list[str]: """Min-cost reference reactions that make ``task`` feasible (RAVEN ``ftINITFillGaps``). @@ -169,7 +213,8 @@ def _gap_fill_task( work.add_cons_vars(extras) # add() over a flat list, not Python sum() — the latter is O(n²) in sympy and with # thousands of candidates dominates gap-fill runtime (see ftINIT/tINIT, same fix). - work.objective = prob.Objective(add(objective_terms), direction="min") + cost_expr = add(objective_terms) + work.objective = prob.Objective(cost_expr, direction="min") _set_fill_solver(work, time_limit, seed) work.slim_optimize() # Accept a near-optimal incumbent (time_limit); only a truly infeasible fill (no @@ -177,8 +222,16 @@ def _gap_fill_task( if work.solver.status not in ("optimal", "feasible", "suboptimal", "time_limit") or \ work.variables[f"_fill_{candidates[0]}"].primal is None: raise OptimizationError(f"gap-filling found no way to make task {task.id!r} feasible.") - return [cid for cid in candidates - if (work.variables[f"_fill_{cid}"].primal or 0.0) > 0.5] + + chosen = [cid for cid in candidates + if (work.variables[f"_fill_{cid}"].primal or 0.0) > 0.5] + # canonical (best-effort): pin the degenerate min-cost fill to the fewest, lowest-id + # reactions so the added set does not depend on the solver seed/version. + if canonical: + canon = _canonicalize_fill(work, prob, candidates, cost_expr, time_limit) + if canon is not None: + chosen = canon + return chosen def fill_tasks( @@ -189,6 +242,7 @@ def fill_tasks( rxn_scores: Mapping[str, float] | None = None, time_limit: float | None = _FILL_TIME_LIMIT, seed: int = _FILL_SEED, + canonical: bool = False, ) -> TaskFillResult: """Add minimum-cost reference reactions so every task is feasible in ``model``. @@ -199,7 +253,9 @@ def fill_tasks( the model, excluding exchange/boundary reactions); ``rxn_scores`` (original reaction id → score) sets each candidate's cost as ``−min(score, −0.1)`` (missing → cost 1). ``should_fail`` tasks are ignored. Each gap-fill MILP is single-threaded with a fixed - ``seed`` and bounded by ``time_limit`` (RAVEN's 300 s). + ``seed`` and bounded by ``time_limit`` (RAVEN's 300 s). ``canonical`` (opt-in) pins the + degenerate min-cost fill to the fewest, lowest-id reactions so the added set does not + depend on the solver seed/version — see :func:`_canonicalize_fill`. Boundary reactions are closed while testing/solving each task, so task inputs and outputs come solely from the task's ranged metabolite bounds (RAVEN gap-fills the exchange-free @@ -225,7 +281,7 @@ def fill_tasks( for r in reference_model.reactions if r.id not in present and not r.boundary} try: chosen = _gap_fill_task(reference_model, present, task, costs, - time_limit=time_limit, seed=seed) + time_limit=time_limit, seed=seed, canonical=canonical) except OptimizationError: failed.append(task.id) continue diff --git a/tests/test_init_ftinit.py b/tests/test_init_ftinit.py index 2a58c43..c434db0 100644 --- a/tests/test_init_ftinit.py +++ b/tests/test_init_ftinit.py @@ -137,3 +137,45 @@ def test_forced_flux_lower_bound_is_respected(): res = run_ftinit(model, scores) assert res.fluxes["R6"] >= 2.0 - 1e-6 assert "R6" not in res.deleted_reactions + + +# --------------------------------------------------------------------------- # +# canonical (deterministic uniqueness) — Tier 2 item 4. +# --------------------------------------------------------------------------- # +def _degenerate_model(): + """Two interchangeable negative-score reactions (R1, R2) both feed an essential E. + + Exactly one is needed to make E carry flux, so the score optimum (-1) is degenerate + between keeping R1 or R2. Which one a plain solve keeps is an arbitrary tie-break. + """ + m = cobra.Model("degen") + a, mm, p = (cobra.Metabolite(x, name=x, compartment="s") for x in ("a", "m", "p")) + m.add_metabolites([a, mm, p]) + R1 = cobra.Reaction("R1", lower_bound=0, upper_bound=1000); R1.add_metabolites({a: -1, mm: 1}) + R2 = cobra.Reaction("R2", lower_bound=0, upper_bound=1000); R2.add_metabolites({a: -1, mm: 1}) + E = cobra.Reaction("E", lower_bound=0, upper_bound=1000); E.add_metabolites({mm: -1, p: 1}) + EXa = cobra.Reaction("EX_a", lower_bound=-1000, upper_bound=1000); EXa.add_metabolites({a: -1}) + EXp = cobra.Reaction("EX_p", lower_bound=-1000, upper_bound=1000); EXp.add_metabolites({p: -1}) + m.add_reactions([R1, R2, E, EXa, EXp]) + m.objective = "E" + return m + + +def test_canonical_breaks_degenerate_tie_by_id(): + """canonical selects the unique sparsest, lowest-id optimum among equal alternatives.""" + m = _degenerate_model() + res = run_ftinit(m, {"R1": -1.0, "R2": -1.0}, essential_rxns=["E"], canonical=True) + # exactly one of the degenerate pair is kept (the tie is resolved, not doubled) ... + assert len({"R1", "R2"} & set(res.kept_reactions)) == 1 + # ... and it is deterministically the lower-id one. + assert "R1" in res.kept_reactions and "R2" not in res.kept_reactions + # the reported objective is the primary (score) optimum, not the phase-2 secondary. + assert res.objective == pytest.approx(-1.0, abs=1e-6) + + +def test_canonical_safe_on_unique_optimum(): + """On a non-degenerate model canonical returns the same optimum (no regression).""" + model = make_test_model() + res = run_ftinit(model, _scores(model), canonical=True) + assert set(res.kept_reactions) == _LOOP + assert res.objective == pytest.approx(8.0, abs=1e-6) diff --git a/tests/test_init_pipeline.py b/tests/test_init_pipeline.py index 910d727..e5dee74 100644 --- a/tests/test_init_pipeline.py +++ b/tests/test_init_pipeline.py @@ -159,3 +159,26 @@ def test_essential_merged_away_is_skipped(): prep = prep_init_model(m, [task], ext_comp="s") # must not raise assert "REV" not in prep.essential_rxns # merged into a collapsed group + + +# --------------------------------------------------------------------------- # +# strict_gap / canonical (deterministic extraction) — Tier 2 items 4 & 5. +# These are opt-in; the default path stays exact-RAVEN. On the toy oracle they must +# reproduce T0001 (no regression) and be run-to-run identical. +# --------------------------------------------------------------------------- # +def test_ftinit_strict_gap_matches_oracle(): + """strict_gap: one near-proven-optimal solve per step, still gives T0001.""" + model = make_test_model() + prep = prep_init_model(model, ext_comp="s") + out = ftinit(prep, _scores(model), strict_gap=True) + assert {r.id for r in out.reactions} == set(TEST_MODEL_FTINIT_NO_TASKS) + + +def test_ftinit_canonical_matches_oracle_and_is_stable(): + """canonical (+ strict_gap): preserves T0001 and is identical across repeated runs.""" + model = make_test_model() + prep = prep_init_model(model, ext_comp="s") + out1 = ftinit(prep, _scores(model), strict_gap=True, canonical=True) + out2 = ftinit(prep, _scores(model), strict_gap=True, canonical=True) + assert {r.id for r in out1.reactions} == set(TEST_MODEL_FTINIT_NO_TASKS) + assert {r.id for r in out1.reactions} == {r.id for r in out2.reactions} diff --git a/tests/test_init_taskfill.py b/tests/test_init_taskfill.py index add5693..2683908 100644 --- a/tests/test_init_taskfill.py +++ b/tests/test_init_taskfill.py @@ -157,3 +157,26 @@ def test_additions_carry_forward_to_later_tasks(): res = fill_tasks(gapped, ref, [first, second]) assert res.added_reactions == ["R7"] # added once; the second task saw R7 already there assert not res.failed_tasks + + +def test_canonical_gap_fill_breaks_tie_by_id(): + """Two equal-cost candidate fills: canonical adds the lower-id one deterministically.""" + import cobra + + from raven_toolbox.tasks import Task + + ref = cobra.Model("ref") + A, M, P = (cobra.Metabolite(x, name=x, compartment="s") for x in ("A", "M", "P")) + ref.add_metabolites([A, M, P]) + RA = cobra.Reaction("RA", lower_bound=0, upper_bound=1000); RA.add_metabolites({A: -1, M: 1}) + RB = cobra.Reaction("RB", lower_bound=0, upper_bound=1000); RB.add_metabolites({A: -1, M: 1}) + RP = cobra.Reaction("RP", lower_bound=0, upper_bound=1000); RP.add_metabolites({M: -1, P: 1}) + ref.add_reactions([RA, RB, RP]) + task = Task(id="mkP", inputs=[("A[s]", 0.0, 1000.0)], outputs=[("P[s]", 1.0, 1000.0)]) + + gapped = ref.copy() + gapped.remove_reactions(["RA", "RB"], remove_orphans=False) # both routes to M removed + # exactly one equal-cost route is needed; canonical must pick the lower id (RA). + res = fill_tasks(gapped, ref, [task], canonical=True) + assert res.added_reactions == ["RA"] + assert not res.failed_tasks From 6a33a0bc2ce24bb6de19c0a888a4916b06f515c6 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Wed, 22 Jul 2026 00:22:02 +0200 Subject: [PATCH 2/2] docs(ftinit): clarify strict_gap/canonical reduce MILP fragility, not biological accuracy --- src/raven_toolbox/init/ftinit.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/raven_toolbox/init/ftinit.py b/src/raven_toolbox/init/ftinit.py index 62acd52..59b08e0 100644 --- a/src/raven_toolbox/init/ftinit.py +++ b/src/raven_toolbox/init/ftinit.py @@ -520,9 +520,11 @@ def ftinit( ``docs/init_param_calibration.md`` for the calibration table. ``strict_gap`` and ``canonical`` (both opt-in, default off → exact RAVEN behaviour) - make the extracted model reproducible across solver versions/platforms instead of - merely across identical runs, which matters when a downstream metric (e.g. gene - essentiality) must not drift on solver noise: + reduce the arbitrariness of the degenerate MILP's tie-break, yielding a more + parsimonious and more *reproducible* extracted model. They do **not** make the model + biologically more accurate: the alternative optima they choose among are equally + consistent with the expression data, so this only pins *which* optimum is returned. It + reduces the fragility of the MILP, not its correctness: * ``strict_gap`` replaces the loose relative-gap escalation with a single solve per step proven to a fixed *absolute* gap (below the 0.1 reaction-score granularity), @@ -532,6 +534,16 @@ def ftinit( lowest-id) optimum, so the degenerate choice is pinned rather than left to the solver — applied both to each extraction step and to the task gap-fill. Best used together with ``strict_gap`` (a well-defined primary optimum to canonicalise). + + Caveats worth carrying into any workflow that uses these: they reduce run-to-run and + platform fragility but do *not* guarantee reproducibility across Gurobi versions + (proving the genome-scale optimum is intractable, so ``strict_gap`` may fall back to an + incumbent), and a downstream metric such as gene essentiality can even shift or worsen, + because the sparser ``canonical`` model is more sensitive to the residual (mostly + transport) degeneracy. For reproducible gene essentiality specifically, pin the solver + stack (raven-toolbox commit + ``gurobipy`` version) rather than relying on these flags. + Rule of thumb: for many models the baseline is fine; reach for these when one or a few + stable, parsimonious model artifacts are wanted. """ if metabolomics: raise NotImplementedError(