From f9334b3df735eae66f26004256f543141abe0d8c Mon Sep 17 00:00:00 2001
From: Eduard Kerkhoven <eduardk@chalmers.se>
Date: Fri, 17 Jul 2026 13:55:53 +0200
Subject: [PATCH] fix: derive the loop-detection threshold from the model
 bounds
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

randomSampling excluded loop reactions by testing maxF > 999 | minF < -999,
a hardcoded ±1000 while the bound replacement eight lines below derives its
threshold from the model. On any model whose bounds are on another scale
(ecModels, 100-scale bounds) nothing matched, so no reaction was excluded and
loops were sampled while the function reported success.

The threshold now comes from the model's own bounds, and reactions that reach
an infinite bound are excluded explicitly. Models bounded at ±1000 keep the
previous behaviour exactly.

Adds a test with a loop that runs to a 100 cap, which the old threshold missed
entirely.
---
 analysis/randomSampling.m          | 14 +++++++++--
 testing/function_tests/tSampling.m | 38 ++++++++++++++++++++++++++++++
 2 files changed, 50 insertions(+), 2 deletions(-)

diff --git a/analysis/randomSampling.m b/analysis/randomSampling.m
index 92b2d0a8..7042370f 100755
--- a/analysis/randomSampling.m
+++ b/analysis/randomSampling.m
@@ -190,8 +190,18 @@
 if isempty(goodRxns)
     [minF, maxF] = getAllowedBounds(model, 'runParallel', runParallel);
     goodRxns = true(nRxns, 1);
-    % Reactions that reach ±1000 are involved in loops
-    goodRxns(maxF > 999 | minF < -999) = false;
+    % Reactions that reach the model's own highest bound are involved in
+    % loops. The threshold has to be derived from the model, as the bound
+    % replacement below does: hardcoding ±1000 matches nothing on a model
+    % whose bounds are on another scale (ecModels, 100-scale bounds), so no
+    % reaction would be excluded and loops would be sampled while the
+    % function reports success.
+    loopUb = max(model.ub);
+    loopLb = min(model.lb);
+    goodRxns(maxF == Inf | maxF > loopUb*0.999) = false;
+    if loopLb < 0
+        goodRxns(minF == -Inf | minF < loopLb*0.999) = false;
+    end
     % Reactions that cannot carry any non-zero flux
     goodRxns(~(maxF > 0 | minF < 0)) = false;
     % In ecModels do not sample from usage_prot reactions
diff --git a/testing/function_tests/tSampling.m b/testing/function_tests/tSampling.m
index 27cf06c3..2e984f20 100644
--- a/testing/function_tests/tSampling.m
+++ b/testing/function_tests/tSampling.m
@@ -136,5 +136,43 @@ function randomSamplingUnknownMethodErrors(testCase)
                 'RAVEN:badInput');
         end
 
+        function randomSamplingLoopDetectionUsesModelBounds(testCase)
+            % Loop detection must take its threshold from the model. R2/R3
+            % form an a -> b -> a loop that runs up to this model's 100 cap,
+            % while the linear path R1 -> R4 is held to 10. A threshold that
+            % assumes ±1000 bounds matches neither, leaving the loop
+            % reactions to be sampled as if they were loop-free.
+            m = tSampling.loopModel();
+            evalc(['[~, goodRxns] = randomSampling(m, 2, ''method'', ' ...
+                '''randomObjective'', ''seed'', 1);']);
+            testCase.verifyEqual(sort(goodRxns(:))', [1 4]);
+        end
+
+    end
+
+    methods (Static, Access = private)
+
+        function m = loopModel()
+            m = struct();
+            m.id       = 'loop';
+            m.rxns     = {'R1';'R2';'R3';'R4'};
+            m.rxnNames = m.rxns;
+            m.mets     = {'a';'b'};
+            m.metNames = {'a';'b'};
+            m.metComps = [1;1];
+            m.comps    = {'c'};
+            m.compNames= {'cytosol'};
+            %             R1  R2  R3  R4
+            m.S = sparse([ 1  -1   1   0;    % a
+                           0   1  -1  -1]);  % b
+            m.lb  = [0;0;0;0];
+            m.ub  = [10;100;100;100];
+            m.rev = [0;0;0;0];
+            m.c   = [0;0;0;1];
+            m.b   = zeros(2,1);
+            m.genes = {}; m.grRules = {'';'';'';''}; m.rxnGeneMat = sparse(4,0);
+            m.metFormulas = {'C';'C'};
+        end
+
     end
 end