fix: derive the loop-detection threshold from the model bounds

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

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