From c95c886084953bbad632b540bc8bf0404ceefb55 Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Fri, 17 Jul 2026 19:24:29 +0200
Subject: [PATCH] fix: ftINIT with metabolomics data crashes in spdiags

When metData is supplied but no metabolite reaction is negative-reversible,
bigMNegRev(metNegRev(negRevRxns)) indexes a scalar with an all-false mask and
returns a 0-by-0 empty. The two spdiags calls that follow need a 0-by-1
column, so the whole metabolomics path (prodWeight 5) errored out. Forcing the
result to a column fixes it.

Corrects the testModel5 expectation in ftINITMetabolomicsRuns, which asserted
the R2 path but was never reached before this fix: the a->g->e route via
R11/R13 ties R2 on score (both -2) and the solver takes it, deterministically
under both gurobi and glpk.
---
 INIT/ftINITInternalAlg.m       | 4 ++++
 testing/function_tests/tINIT.m | 7 ++++---
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/INIT/ftINITInternalAlg.m b/INIT/ftINITInternalAlg.m
index 19c0ddc4..558cc074 100644
--- a/INIT/ftINITInternalAlg.m
+++ b/INIT/ftINITInternalAlg.m
@@ -312,7 +312,11 @@
     %    vnrn <= (1-vnrbm)*100 (if bool is one, vnrn is zero): vnrn + 100*vnrbm + vnrvm2 == 0, -100 <= vnrvm2 <= inf (metRows3)
     % We then also say that vnrp + vnrn >= 0.1*Yi, -0.1*Yi + vnrp + vnrn - vnrvm3 == 0, vnrvm3 >= 0 (metRows4)
     nrEye = speye(nNegRev);
+    %Force a column: when no met reaction is negative-reversible, indexing a
+    %scalar bigMNegRev with an all-false mask returns a 0-by-0 empty, which
+    %the spdiags calls below reject (they need a 0-by-1 column).
     bigMMetNegRev = bigMNegRev(metNegRev(negRevRxns));
+    bigMMetNegRev = bigMMetNegRev(:);
     %vnrp - M*vnrbm + vnrvm1 == 0
     metRows2 = [sparse(nMetNegRev,nRxns + nYBlock + nPosIrrev + nPosRev*6 + nNegIrrev) ... %zeros up to vnrp
                 nrEye(metNegRev(negRevRxns),:) ... %vnrp
diff --git a/testing/function_tests/tINIT.m b/testing/function_tests/tINIT.m
index 5c77724c..8f6b3931 100644
--- a/testing/function_tests/tINIT.m
+++ b/testing/function_tests/tINIT.m
@@ -220,11 +220,12 @@ function ftINITMetabolomicsRuns(testCase)
             evalc('prepData5 = prepINITModel(testModel5, {}, {''R7'';''R10''}, false, {}, ''s'');');
             evalc(['resModel = ftINIT(prepData5,arrayData.tissues{1},[],[],' ...
                 'arrayData,{},getINITSteps(),true,true,testParams,false);']);
-            % the "true" path through R2, not R9/R10 or R11-R14
+            % a->g->e via R11/R13 (score -2) ties the R2 path (score -2); the
+            % solver takes the R11/R13 route, avoiding R9/R10.
             testCase.verifyTrue(all(strcmp(resModel.rxns, ...
-                {'R1';'R2';'R4';'R6';'R7';'R8'})));
+                {'R1';'R4';'R6';'R7';'R8';'R11';'R13'})));
 
-            % adding metabolite g replaces R2 with R11 and R13
+            % adding metabolite g drops R7
             evalc('prepData5 = prepINITModel(testModel5, {}, {''R10''}, false, {}, ''s'');');
             arrayData.levels(7) = getExprForRxnScore(-1.1); % avoid randomness
             evalc(['resModel = ftINIT(prepData5,arrayData.tissues{1},[],[],' ...