feat(ftinit): opt-in deterministic extraction (strict gap + canonical optimum)

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,32 @@ 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)
+    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),
+        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).
+
+    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(
@@ -436,6 +578,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 +603,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