fix: align raven-toolbox with RAVEN (parity review)

.github/workflows/parity-nightly.yml

@@ -57,7 +57,7 @@ jobs:
           cache: pip
           cache-dependency-path: pyproject.toml
 
-      - run: pip install --upgrade pip
+      - run: python -m pip install --upgrade pip
       - run: pip install -e ".[dev]" gurobipy
 
       # Fail here, loudly, rather than let the tests skip and the job go green

.github/workflows/release.yml

@@ -31,7 +31,7 @@ jobs:
           python-version: "3.12"
           cache: pip
           cache-dependency-path: pyproject.toml
-      - run: pip install --upgrade pip
+      - run: python -m pip install --upgrade pip
       - run: pip install build
       - run: python -m build
 

CHANGELOG.md

@@ -4,6 +4,18 @@ Milestones in the raven-toolbox port. For function-level status see
 [docs/raven_migration.md](https://github.com/CompMetEng/raven-toolbox/blob/develop/docs/reference/migration.md); for open work see
 [docs/todo.md](https://github.com/CompMetEng/raven-toolbox/blob/develop/docs/reference/todo.md).
 
+## Unreleased
+
+* `init.ftinit` and `init.run_ftinit` give each scored reaction its own big-M, its flux
+  capacity `max(min(ub, big_m), min(-lb, big_m), force_on)` with an infinite bound counted as
+  `big_m`, as RAVEN's `ftINITInternalAlg` does. `big_m` is the cap on that capacity; its
+  default is 1000 (was 100).
+* `tasks.find_task_essential_reactions` (and so `init.prep_init_model`) re-solves a knock-out
+  from a cold start when the warm-started solve does not report a feasible status, and
+  counts the reaction as essential only if that solve is not feasible either, as RAVEN's
+  `solveLP` does. Gurobi could report a feasible knock-out as infeasible when started from
+  the previous knock-out's basis, which added reactions to the essential set.
+
 ## 3.0.0b1 — 2026-09-07
 
 First beta toward 3.0, matching the upcoming MATLAB RAVEN 3.0 beta.

scripts/analyze_init_params.py

@@ -20,7 +20,7 @@
 * ``ftinit_milp``  — single staged-MILP step (step 0 of series '1+1') on the merged model.
                      Cheap (~30-200 s each); the core sweep for ``mip_gap``/``big_m``/``force_on``.
 * ``prep_scale``   — rescaleModelForINIT on/off and its ``max_stoich_diff``, fed into the
-                     same step-0 MILP. Shows why scaling is needed for a fixed big-M.
+                     same step-0 MILP. Shows why scaling is needed when the big-M is capped.
 * ``ftinit_full``  — the whole ``ftinit`` pipeline (both steps + gap-fill). Sweeps
                      ``mip_gap``/``big_m``. Expensive (~200 s+/config).
 
@@ -185,7 +185,7 @@ def cfg(label, **kw):
             save()
         return store[key]
 
-    # mip_gap sweep (big_m=100, force_on=0.1)
+    # mip_gap sweep (force_on=0.1; big_m=100 caps the per-reaction big-M)
     res = [cfg(f"gap={g}", mip_gap=g, big_m=100.0, force_on=0.1, time_limit=900) for g in MIP_GAPS]
     _finalize(res)
     doc += _print_table("ftINIT step-0: mip_gap (big_m=100, force_on=0.1)", res)
@@ -194,9 +194,10 @@ def cfg(label, **kw):
     res = [cfg(f"big_m={int(b)}", mip_gap=0.001, big_m=b, force_on=0.1, time_limit=900) for b in BIG_MS]
     _finalize(res)
     doc += _print_table("ftINIT step-0: big_m (gap=0.001, force_on=0.1)", res,
-                        "big_m caps a scored reaction's flux; large values weaken the LP relaxation.")
+                        "big_m caps the per-reaction big-M (a scored reaction's flux capacity); "
+                        "large values weaken the LP relaxation.")
 
-    # force_on sweep (gap=0.001, big_m=100) — changes the model (connectivity threshold)
+    # force_on sweep (gap=0.001, big_m=100 as the big-M cap); changes the model (connectivity threshold)
     res = [cfg(f"force_on={fo}", mip_gap=0.001, big_m=100.0, force_on=fo, time_limit=900) for fo in FORCE_ONS]
     _finalize(res)
     doc += _print_table("ftINIT step-0: force_on (gap=0.001, big_m=100)", res,
@@ -236,7 +237,7 @@ def cfg(label, msd):
     res.append(cfg("scale=off", None))
     _finalize(res)
     doc += _print_table("prep scaling: rescaleModelForINIT max_stoich_diff (+scale off), big_m=100", res,
-                        "With big_m=100 fixed, scale=off / poor conditioning is expected to be "
+                        "With the big-M capped at big_m=100, scale=off / poor conditioning is expected to be "
                         "infeasible or far slower — that is the reason scaling is on by default.")
     return doc
 
@@ -250,7 +251,7 @@ def cfg(label, **kw):
             print(f"[ftinit_full] {label} ...", flush=True)
             t = time.time()
             try:
-                out = ftinit(prep, rxn_scores, gene_scores=gene_scores, series="1+1", **kw)
+                out = ftinit(prep, rxn_scores=rxn_scores, gene_scores=gene_scores, series="1+1", **kw)
                 r = Result(label=label, seconds=time.time() - t, status="ok",
                            objective=0.0, n_kept=len(out.reactions),
                            reactions=sorted(x.id for x in out.reactions))

scripts/analyze_init_robustness.py

@@ -188,7 +188,7 @@ def cached(key, fn):
     def model_for(e, **kw):
         g = gene_scores_from_expression(e, 1.0)
         r = score_reactions_from_genes(ref, g, no_gene_score=kw.get("no_gene_score", -2.0))
-        return ftinit(prep, r, gene_scores=g, series="1+1",
+        return ftinit(prep, rxn_scores=r, gene_scores=g, series="1+1",
                       force_on=kw.get("force_on", 0.1), mip_gap=MIP_GAP, time_limit=TIME_LIMIT)
 
     phases = set(args.phase.split(","))

scripts/analyze_init_solvers.py

@@ -101,7 +101,7 @@ def run(solver: str) -> Result:
             local_prep.min_model.solver = solver
             g = gene_scores_from_expression(expr, 1.0)
             r = score_reactions_from_genes(ref, g)
-            model = ftinit(local_prep, r, gene_scores=g, series="1+1",
+            model = ftinit(local_prep, rxn_scores=r, gene_scores=g, series="1+1",
                            mip_gap=args.mip_gap, time_limit=args.time_limit)
             rset = sorted(x.id for x in model.reactions)
             res = Result(solver, time.time() - t, "ok", len(rset), rset)

src/raven_toolbox/analysis/fseof.py

@@ -1,21 +1,19 @@
-"""Flux Scanning based on Enforced Objective Flux — FSEOF (port + redesign).
+"""Flux Scanning based on Enforced Objective Flux (FSEOF).
 
 FSEOF (Choi et al., Appl Environ Microbiol 2010) finds metabolic-engineering targets
 for over-producing a metabolite: enforce an increasing flux toward the target product
-while optimising growth, and watch how each reaction's flux responds. This is a port
-of RAVEN's ``FSEOF`` with a substantially richer, more robust output (RAVEN's
-weaknesses are noted in IMPROVEMENTS, FS1–FS4):
-
-* **Robust trend, not strict monotonicity.** Each reaction's flux is regressed against
-  the enforced product flux across the scan; the **slope** is the response and the
-  **correlation** (|r|) is a quality score. A reaction is a target if it tracks the
-  product cleanly (|r| ≥ ``correlation_threshold``), so one noisy step from LP
-  alternative optima does not discard it (and pFBA per step keeps the scan stable).
-* **Direction classification RAVEN lacks.** Targets are labelled ``amplify`` (|flux|
-  rises with the product → over-express), ``knockdown`` (|flux| falls), or ``knockout``
-  (|flux| → ~0 → delete). RAVEN only ever reports the amplification targets.
-* **Gene-level view** via :attr:`FSEOFResult.gene_targets`, and the full flux scan is
-  retained in :attr:`FSEOFResult.scan` — all as DataFrames, not a printed TSV.
+while optimising growth, and watch how each reaction's flux responds. This follows
+RAVEN's ``FSEOF`` and returns tables instead of a printed list:
+
+* **Trend of the absolute flux.** Each reaction's absolute flux is regressed against the
+  enforced product flux across the scan; the **slope** is the response and the
+  **correlation** (|r|) is a quality score. A reaction is a target when its track is not
+  constant and |r| >= ``correlation_threshold``. Every step is solved with pFBA.
+* **Direction classification.** Targets are labelled ``amplify`` (the absolute flux
+  rises with the product), ``knockdown`` (it falls) or ``knockout`` (it falls to at most
+  1e-8), like the Amplify, Knockdown and Knockout labels of RAVEN's ``FSEOF``.
+* **Gene-level view** via :attr:`FSEOFResult.gene_targets`, and the full flux scan in
+  :attr:`FSEOFResult.scan`, both as DataFrames.
 """
 from __future__ import annotations
 
@@ -97,7 +95,17 @@ def fseof(
         product flux for it to be reported as a target; raise it to keep only
         the most cleanly-trending reactions.
     flux_eps:
-        Tolerance below which a flux or slope is treated as zero.
+        Tolerance below which the theoretical maximum of the target flux is treated as zero
+        (the scan then raises). Reaction fluxes are classified as in RAVEN's ``FSEOF``: a
+        reaction is a target when some absolute flux exceeds 1e-10, the standard deviation of
+        its absolute flux is above 0 and the correlation reaches ``correlation_threshold``;
+        it is ``knockout`` when the slope is negative and the last absolute flux is at most
+        1e-8, ``knockdown`` for a negative slope otherwise.
+
+        Differences from ``FSEOF.m``: the scan stops at the first infeasible enforced level,
+        so the last absolute flux used for the knockout label is that of the last feasible
+        level; ``targets`` is sorted by ``|slope| * |r|`` and not grouped by label; and
+        ``flux_eps`` only gates the check of the theoretical maximum.
     min_target:
         Floor of the enforced-flux range, evenly spaced up to the scan's ceiling
         (``max_fraction`` of the theoretical maximum) in ``n_steps`` steps.
@@ -127,40 +135,39 @@ def fseof(
             try:
                 columns[level] = pfba(model).fluxes
             except OptimizationError:
-                break  # enforced flux became infeasible — stop scanning
+                break  # enforced flux became infeasible: stop scanning
             enforced.append(level)
     if len(enforced) < 2:
         raise RuntimeError("FSEOF needs at least two feasible enforced-flux levels.")
 
     scan = pd.DataFrame(columns)
-    targets = _classify(model, scan, np.asarray(enforced), correlation_threshold, flux_eps)
+    targets = _classify(model, scan, np.asarray(enforced), correlation_threshold)
     return FSEOFResult(scan=scan, enforced=enforced, targets=targets)
 
 
-def _classify(model, scan, enforced, corr_threshold, flux_eps) -> pd.DataFrame:
+_TRACK_EPS = 1e-10  # a track must exceed this somewhere to be fitted, as in FSEOF.m
+_KNOCKOUT_EPS = 1e-8  # last absolute flux at or below this is a knockout, as in FSEOF.m
+
+
+def _classify(model, scan, enforced, corr_threshold) -> pd.DataFrame:
     rows = []
     for rxn in model.reactions:
         flux = scan.loc[rxn.id, enforced.tolist() if hasattr(enforced, "tolist") else enforced]
         flux = flux.to_numpy(dtype=float)
         initial, final = flux[0], flux[-1]
-        if flux.std() < flux_eps:  # flat -> no response
+        # Slope and correlation are taken on absolute flux, as in RAVEN's FSEOF, so a
+        # reversible reaction running backwards counts by its magnitude.
+        magnitude = np.abs(flux)
+        if not (magnitude > _TRACK_EPS).any() or magnitude.std() <= 0:
             continue
-        fit = linregress(enforced, flux)
+        fit = linregress(enforced, magnitude)
         slope, corr = float(fit.slope), float(fit.rvalue)
-        if abs(corr) < corr_threshold or abs(slope) < flux_eps:
+        if abs(corr) < corr_threshold or slope == 0:
             continue
-        # Classify on the slope of |flux| vs the enforced product flux — the
-        # criterion the docstring states (|flux| rises = amplify, etc.), not an
-        # endpoint-only ``abs(final) vs abs(initial)`` comparison, which can
-        # mislabel a track whose first/last values straddle a peak/trough even
-        # though its overall trend runs the other way. Keep ``knockout`` for
-        # tracks the regression drives essentially to zero.
-        abs_fit = linregress(enforced, np.abs(flux))
-        abs_slope = float(abs_fit.slope)
-        if abs(final) < flux_eps and abs_slope < 0:
-            ttype = "knockout"
-        elif abs_slope > 0:
+        if slope > 0:
             ttype = "amplify"
+        elif abs(final) <= _KNOCKOUT_EPS:
+            ttype = "knockout"
         else:
             ttype = "knockdown"
         rows.append({

src/raven_toolbox/analysis/reporter.py

@@ -47,8 +47,17 @@ def _gene_z(pvalues: dict[str, float]) -> dict[str, float]:
 
 def _reporter_one(model: cobra.Model, gene_z: dict[str, float], test: str) -> ReporterResult:
     z_values = np.fromiter(gene_z.values(), dtype=float)
-    mu = float(z_values.mean()) if z_values.size else 0.0
-    sigma = float(z_values.std(ddof=0)) if z_values.size else 0.0
+    # Moments are taken on z-scores centred on c, the first z-score plus the mean of the
+    # z-scores relative to it (two passes). A spread of a few ulps between otherwise identical
+    # z-scores is then not lost in the rounding of the mean, also when the first z-score is an outlier.
+    z_first = float(z_values[0]) if z_values.size else 0.0
+    centre = z_first + float((z_values - z_first).mean()) if z_values.size else 0.0
+    centred = z_values - centre
+    mu_c = float(centred.mean()) if z_values.size else 0.0
+    sigma = float(centred.std(ddof=0)) if z_values.size else 0.0
+    # The spread is tested for being 0 up to 100 * eps(max|z|), as in RAVEN.
+    spread = float(z_values.max() - z_values.min()) if z_values.size else 0.0
+    has_spread = spread > 100 * float(np.spacing(np.abs(z_values).max())) if z_values.size else False
 
     rows = []
     for met in model.metabolites:
@@ -57,9 +66,9 @@ def _reporter_one(model: cobra.Model, gene_z: dict[str, float], test: str) -> Re
             continue
         zs = np.array([gene_z[g] for g in neighbours])
         n = zs.size
-        raw = zs.sum() / math.sqrt(n)
+        raw_c = (zs - centre).sum() / math.sqrt(n)
         # Exact background correction for sampling-with-replacement (see module doc).
-        corrected = (raw - math.sqrt(n) * mu) / sigma if sigma > 0 else 0.0
+        corrected = (raw_c - math.sqrt(n) * mu_c) / sigma if has_spread else 0.0
         if not math.isfinite(corrected):  # never on clamped z-scores; guards the p-value
             corrected = 0.0
         rows.append(
@@ -70,7 +79,7 @@ def _reporter_one(model: cobra.Model, gene_z: dict[str, float], test: str) -> Re
                 "p_value": float(1.0 - norm.cdf(corrected)),
                 "n_genes": n,
                 "mean_z": float(zs.mean()),
-                "std_z": float(zs.std(ddof=1)) if n > 1 else float("nan"),
+                "std_z": float(zs.std(ddof=1)) if n > 1 else 0.0,
             }
         )
     table = pd.DataFrame(rows, columns=["metabolite", "name", "z_score", "p_value", "n_genes", "mean_z", "std_z"])
@@ -86,11 +95,14 @@ def reporter_metabolites(
 ) -> list[ReporterResult]:
     """Compute Reporter Metabolites from per-gene differential-expression p-values.
 
-    ``gene_pvalues`` maps gene id → p-value (genes not in the model, or with a NaN or
-    out-of-``[0, 1]`` p-value, are dropped — a stray invalid p-value would otherwise
-    turn the whole result NaN). If ``gene_fold_changes`` (gene id → log fold change)
-    is given, two extra results are returned for the up- (fc ≥ 0) and down- (fc < 0)
-    regulated gene subsets, in addition to ``"all"``.
+    ``gene_pvalues`` maps gene id → p-value. Genes not in the model and genes with a
+    NaN (or ``None``) p-value are dropped. Any other p-value outside ``[0, 1]``,
+    including ±inf and for genes that are not in the model, raises ``ValueError``
+    (RAVEN raises ``RAVEN:badInput``). If ``gene_fold_changes`` (gene id → log fold change)
+    has at least one non-zero value, two extra results are returned for the up- (fc ≥ 0)
+    and down- (fc < 0) regulated gene subsets, in addition to ``"all"``; with all fold
+    changes zero (or none given) only ``"all"`` is returned, as in RAVEN. A metabolite
+    with a single scored neighbouring gene has ``std_z`` 0.
 
     The background correction uses the population standard deviation
     (``ddof=0``) of the scored Z-scores, which the closed form requires. The
@@ -104,16 +116,25 @@ def reporter_metabolites(
     z-scores smaller by a factor ``√(N/(N−1))`` for N scored genes; the
     ranking is the same.
     """
+    given = {g: float(p) for g, p in gene_pvalues.items() if p is not None and not math.isnan(p)}
+    invalid = [(g, p) for g, p in given.items() if not 0.0 <= p <= 1.0]
+    if invalid:
+        gene, value = invalid[0]
+        msg = (
+            f"The P-values must lie between 0 and 1, but {len(invalid)} of them do not. "
+            f"The first is {value:g} (gene {gene})"
+        )
+        raise ValueError(msg)
     model_genes = {g.id for g in model.genes}
-    scored = {
-        g: float(p)
-        for g, p in gene_pvalues.items()
-        if g in model_genes and p is not None and not math.isnan(p) and 0.0 <= p <= 1.0
-    }
+    scored = {g: p for g, p in given.items() if g in model_genes}
     gene_z = _gene_z(scored)
     results = [_reporter_one(model, gene_z, "all")]
 
-    if gene_fold_changes is not None:
+    # RAVEN's any(geneFoldChanges) runs after genes outside the model or without a
+    # p-value are dropped, and ignores NaN.
+    if gene_fold_changes is not None and any(
+        v != 0 and not math.isnan(v) for g, v in gene_fold_changes.items() if g in scored
+    ):
         up = {g: z for g, z in gene_z.items() if gene_fold_changes.get(g, 0.0) >= 0}
         down = {g: z for g, z in gene_z.items() if gene_fold_changes.get(g, 0.0) < 0}
         results.append(_reporter_one(model, up, "up"))

src/raven_toolbox/biomass/gam.py

@@ -22,11 +22,13 @@ def set_gam(
 
     Scales every metabolite participating in the biomass pseudoreaction
     whose ``name`` is in ``cofactor_met_names`` (e.g. ATP, ADP, H2O,
-    H+, phosphate) to ``±value``, preserving the sign of its current
+    H+, phosphate) to exactly ``±value``, preserving the sign of its current
     coefficient.
 
     If both ``ngam_rxn`` and ``ngam_value`` are given, the NGAM
-    reaction's bounds are fixed at ``(ngam_value, ngam_value)``.
+    reaction's bounds are fixed at ``(ngam_value, ngam_value)``. An empty
+    ``ngam_rxn`` counts as not given, and ``ngam_value`` without an ``ngam_rxn`` is an
+    error. A given ``ngam_rxn`` must be in the model, also when no ``ngam_value`` is given.
 
     Parameters
     ----------
@@ -50,21 +52,30 @@ def set_gam(
     -------
     cobra.Model
         The (mutated) model, for chaining.
+
+    Raises
+    ------
+    KeyError
+        If ``biomass_rxn``, or a given ``ngam_rxn``, is not in the model. Both reactions
+        are looked up before anything is changed.
+    ValueError
+        If ``ngam_value`` is given without an ``ngam_rxn``.
     """
+    if ngam_value is not None and not ngam_rxn:
+        raise ValueError("ngam_value needs an ngam_rxn: the NGAM reaction to set is not given")
     rxn = model.reactions.get_by_id(biomass_rxn)
+    ngam = model.reactions.get_by_id(ngam_rxn) if ngam_rxn else None
     cofactor_set = set(cofactor_met_names)
 
-    deltas: dict[cobra.Metabolite, float] = {}
+    targets: dict[cobra.Metabolite, float] = {}
     for met, coef in rxn.metabolites.items():
         if met.name in cofactor_set:
             target = (1 if coef > 0 else -1) * float(value)
-            # combine=True adds delta → new total = target
-            deltas[met] = target - coef
-    if deltas:
-        rxn.add_metabolites(deltas, combine=True)
+            targets[met] = target
+    if targets:
+        rxn.add_metabolites(targets, combine=False)  # replaces the coefficient, so it is exactly target
 
-    if ngam_rxn is not None and ngam_value is not None:
-        ngam = model.reactions.get_by_id(ngam_rxn)
+    if ngam is not None and ngam_value is not None:
         ngam.bounds = (float(ngam_value), float(ngam_value))
 
     return model