edkerk · Oct 17, 2026 · Oct 17, 2026 · Oct 17, 2026 · Oct 17, 2026 · Oct 17, 2026
diff --git 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
@@ -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

diff --git 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 = {

diff --git 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
-        <https://github.com/edkerk/raven-docs/blob/main/docs/parameter-tuning/studies/sampling-convergence-calibration.md>`_
+        <https://github.com/CompMetEng/raven-docs/blob/main/docs/parameter-tuning/studies/sampling-convergence-calibration.md>`_
         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.

diff --git 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/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[<header>] = <cell>``).
+
+``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
@@ -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,19 +394,40 @@ 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
 
 
 # --- 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/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_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
@@ -21,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,
@@ -38,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/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