Deterministic, score-aligned reaction placement (fixes 52.8% agreement regression)

CHANGELOG.md

@@ -6,6 +6,19 @@ Milestones in the raven-toolbox port. For function-level status see
 
 ## Unreleased
 
+* **Deterministic, score-aligned reaction placement in `assign_compartments`.** The placement master
+  maximises a gene-localisation objective that never mentions the per-reaction placement variable, so
+  each reaction's compartment was a free co-optimum: the pinned solver returned it reproducibly but
+  *arbitrarily*, and reaction agreement with curated yeast-GEM was only **52.8 %** (an earlier un-pinned
+  build happened to land ~72 %). A lexicographic second pass now fixes the gene layout to the primary
+  optimum, then places each reaction in the compartment its own enzymes are predicted to occupy — the
+  summed DeepLoc score of the reaction's genes, with a small `default_compartment` prior so genes-free
+  and score-tied reactions fall there deterministically. Yeast reaction agreement rises to **72.5 %**
+  (1408/1943) and now rests on the localisation evidence rather than a solver tie-break; **gene agreement
+  is unchanged** (88.7 %, 716/807 — the gene layout is untouched); coherent placement fragments fewer
+  metabolites, so it adds **fewer transports** (1001 → 967); growth and blocked-fraction are unchanged;
+  the second solve is warm-started and adds ~4 s. Reproducible run-to-run (identical placement across
+  independent runs).
 * **`diff_models` compares grRules as logic, not text.** The GPR check now DNF-expands each rule (via the
   existing `manipulation.gpr_to_dnf`), sorts the genes within each isozyme clause and sorts the clauses, so
   operand order no longer registers as a difference: `a and b` == `b and a` and `a or b` == `b or a`. The

docs/studies/localization_redesign.md

@@ -182,6 +182,13 @@ Reading:
 - **The honest agreement-vs-transport frontier is now measurable.** 72.5 % agreement costs ~1261
   reconnecting transports (983 FVA-usable), against curated yeast-GEM's own 1467 inter-compartment
   reactions — high agreement honestly requires many transports (there is no free lunch).
+- **Reaction placement is deterministic and score-aligned, not an arbitrary co-optimum.** The master's
+  objective is over gene-localisation variables only, so each reaction's compartment is a free
+  co-optimum — a pinned solver returns it reproducibly, but *arbitrarily* (an un-tie-broken solve scored
+  only 52.8 %). A lexicographic second pass fixes the gene layout to the primary optimum, then places
+  each reaction in the compartment its own enzymes are predicted to occupy (summed DeepLoc score of the
+  reaction's genes, small default-compartment prior for genes-free and tied reactions). The gene layout
+  is untouched, so the 72.5 % rests on the localisation evidence itself, not on a solver's tie-break.
 
 Every number is authoritative because the growth column *is* the certificate: there is no separate MILP
 floor to disagree with it. Raw JSON in `.research_tmp/certify_yeast_gem.json`.

docs/studies/multiorganism_validation.md

@@ -19,6 +19,13 @@ only a compartment-label map and a metabolite base-key per organism.
 | *A. thaliana* | plant | AraCore v2.1 | 585 | 6 | **yes** | 35 s | **yes** | **81.8 %** | **93.7 %** | 215 | 0.0 % |
 | *C. reinhardtii* | algal | iCre1355 | 2394 | 11 | **yes** | 93 s | **yes** | 53.4 % | 70.5 % | 1108 | 19.0 % |
 
+> **Reaction-agreement column predates the deterministic placement tie-break** (Unreleased; see
+> [yeast_validation.md](yeast_validation.md)). Before it, the per-reaction compartment was a free
+> co-optimum the solver resolved arbitrarily; it is now placed by each reaction's own enzyme scores.
+> These reaction-agreement figures will be re-measured on the new placement (yeast, re-run, went 72.0 %
+> → 72.5 %). The **gene**-agreement column, growth, and blocked fraction are unaffected — the tie-break
+> leaves the gene layout untouched.
+
 Agreement is scored against each model's **own** compartment annotation (a *circular* truth — the model
 being re-placed is also the reference). This measures self-consistency and functionality, not accuracy;
 the independent-ground-truth work (P2: HPA for human, SUBA5 for Arabidopsis, the Chlamydomonas

docs/studies/yeast_validation.md

@@ -21,19 +21,29 @@ parent organelle, matching DeepLoc's resolution).
 
 | metric | `assign_compartments` |
 |---|--:|
-| reaction-level agreement | **72.0 %** (1400/1943) |
+| reaction-level agreement | **72.5 %** (1408/1943) |
 | gene-level agreement | **88.7 %** (716/807) |
-| transports added | 1001 (curated 1467) |
+| transports added | 967 (curated 1467) |
 | blocked re-placed reactions | 29.7 % |
 | materialised growth | 0.1426 /h (floor 0.040) |
 | certified (real FBA) | **yes** |
 
 Gene agreement (88.7 %) sits well above DeepLoc's ~64 % raw ceiling because functionality corrects
-placements; the model stays functional because the method provisions the ~1000 transports the placement
+placements; the model stays functional because the method provisions the ~970 transports the placement
 actually needs, and the blocked 29.7 % is yeast-GEM's *intrinsic* dead-end floor, not a placement
 artefact. Every number is authoritative because the growth column is a real FBA on the materialised
 model — the certificate, not a separate solver floor.
 
+Reaction placement is **deterministic and score-aligned.** The placement master maximises a
+gene-localisation objective, which never mentions the per-reaction placement variable — so each
+reaction's compartment was a free co-optimum: reproducible once the solver is pinned, but an *arbitrary*
+vertex (which on its own scored only 52.8 %). A lexicographic second pass fixes the gene layout to the
+primary optimum, then places each reaction in the compartment its own enzymes are predicted to occupy
+(the summed DeepLoc score of the reaction's genes, with a small default-compartment prior for genes-free
+and score-tied reactions). The gene layout — hence gene agreement — is untouched; the 72.5 % reaction
+agreement now rests on the localisation evidence itself, not on a solver's arbitrary tie-break, and the
+coherent placement needs fewer transports (967 vs 1001).
+
 ## 2. Head-to-head with CarveFungi (Comparison 2)
 
 CarveFungi's own `minmax_reduction` carve (CPLEX, unmodified) on its real universal-DB candidate set +

src/raven_toolbox/localization/certify.py

@@ -158,8 +158,14 @@ def _score(score_df, g: str, c: str) -> float:
     return 0.0 if v is None or (isinstance(v, float) and math.isnan(v)) else float(v)
 
 
+# Tie-break weight for the placement master's second pass: small enough never to override a real
+# per-reaction score difference (DeepLoc scores are ~[0, 1] summed over a reaction's genes), just enough
+# to send a genes-free or score-tied reaction to the default compartment deterministically.
+_DEFAULT_COMPARTMENT_PRIOR = 1e-3
+
+
 def _solve_placement_master(
-    model, movable, genes_in_scope, gene_rxns, score_df, compartments, *,
+    model, movable, genes_in_scope, gene_rxns, score_df, compartments, default_compartment, *,
     multi_compartment_penalty, forced, colocation_groups, time_limit,
 ):
     """Flux-free score-maximising placement MILP (mono-localisation).
@@ -169,6 +175,11 @@ def _solve_placement_master(
     compartment (``x[a, c] = x[b, c]``), letting the score objective pick which one. No flux and no
     growth floor at all in the master — the placement can never harvest a compartment's score through
     leaked flux.
+
+    Solved in two lexicographic passes: the primary objective places genes by score (and penalises
+    spread), which leaves the per-reaction placement ``x`` a free co-optimum; the second pass then fixes
+    the gene layout and places each reaction in the compartment its own enzymes score highest for, so the
+    reaction placement is deterministic and evidence-aligned rather than an arbitrary co-optimal vertex.
     """
     model.solver  # noqa: B018 — initialise the solver so model.problem is usable
     prob = model.problem
@@ -212,21 +223,52 @@ def _solve_placement_master(
 
     opt.add(list(x.values()) + list(y.values()) + cons)
 
-    obj = []
+    # Primary objective: place each gene in the compartment(s) its DeepLoc score favours, penalising
+    # spread across compartments. This fixes the gene layout but leaves the *reaction* placement x a free
+    # co-optimum (the objective never mentions x), which a deterministic solver resolves to an arbitrary
+    # vertex -- yeast-GEM reaction agreement 52.8%. A lexicographically-lower second pass then chooses x.
+    primary = []
     for g in genes_sorted:
         for c in compartments:
             s = _score(score_df, g, c)
             if s:
-                obj.append(mul([Real(s), y[g, c]]))
+                primary.append(mul([Real(s), y[g, c]]))
     if multi_compartment_penalty:
         for v in y.values():
-            obj.append(mul([Real(-multi_compartment_penalty), v]))
-    opt.objective = prob.Objective(add(obj) if obj else Real(0.0), direction="max")
+            primary.append(mul([Real(-multi_compartment_penalty), v]))
+    opt.objective = prob.Objective(add(primary) if primary else Real(0.0), direction="max")
     if time_limit is not None:
         opt.configuration.timeout = int(time_limit)
     _pin_deterministic(prob, opt)
     opt.optimize()
+    if opt.status not in ("optimal", "feasible", "suboptimal", "time_limit"):
+        return opt.status, {}, {}
 
+    # Lexicographic second pass. Fix the gene layout to the primary optimum -- fixing the *solution* (each
+    # y binary), not the objective value, so there is no near-optimal tolerance to tune -- then place each
+    # reaction in the compartment its own enzymes are predicted to occupy: reward x[r, c] by the summed
+    # gene score of r's genes for c, with a small prior for `default_compartment` so genes-free and
+    # score-tied reactions fall there deterministically rather than to an arbitrary co-optimal vertex.
+    # The gene layout (agreement, multi-compartment consolidation) is untouched; only reaction placement,
+    # which was free, is now meaningful.
+    try:  # read every primal before touching any bound -- the first bound change discards the solution
+        y_star = {k: 1.0 if (v.primal or 0.0) >= 0.5 else 0.0 for k, v in y.items()}
+    except (AttributeError, ValueError):
+        return "no_incumbent", {}, {}
+    for k, v in y.items():
+        v.lb = v.ub = y_star[k]
+    secondary = []
+    for r in movable:
+        r_genes = sorted({gg.id for gg in r.genes} & genes_in_scope)
+        for c in compartments:
+            w = sum(_score(score_df, g, c) for g in r_genes)
+            if c == default_compartment:
+                w += _DEFAULT_COMPARTMENT_PRIOR
+            if w:
+                secondary.append(mul([Real(w), x[r.id, c]]))
+    opt.objective = prob.Objective(add(secondary) if secondary else Real(0.0), direction="max")
+    _pin_deterministic(prob, opt)
+    opt.optimize()
     if opt.status not in ("optimal", "feasible", "suboptimal", "time_limit"):
         return opt.status, {}, {}
     try:
@@ -621,7 +663,7 @@ def assign_compartments(
         seen.add(signature)
 
         status, placements, gene_comps = _solve_placement_master(
-            model, movable, genes_in_scope, gene_rxns, score_df, compartments,
+            model, movable, genes_in_scope, gene_rxns, score_df, compartments, default_compartment,
             multi_compartment_penalty=multi_compartment_penalty,
             forced=forced, colocation_groups=groups, time_limit=time_limit)
         if not placements: