RAVEN 2.11.3

INIT/INITStepDesc.m

@@ -24,7 +24,9 @@
                    %      .mets        Names of metabolites to remove
                    %      .compsToKeep Compartments for which metabolites should be kept.
       MILPParams %Cell array of MILPparams - dictates how many iterations that will be run in this step.
-                 %Typically, MIPGap and TimeLimit is specified
+                 %Typically, MIPGap and TimeLimit is specified. ftINIT additionally
+                 %defaults Threads to 1 (single-threaded Gurobi) for determinism;
+                 %set Threads here to override (0 = use all cores).
       AbsMIPGaps %If the objective is close to zero, a percentage of that is very small. 
                  %Therefore, also set an absolut value for this (typically 10 or 20). 
                  %For practical reasons, the first number is not used

INIT/ftINIT.m

@@ -233,6 +233,13 @@
         if ~isfield(params, 'TimeLimit')
             params.TimeLimit = 5000;
         end
+
+        %Default to single-threaded Gurobi MILP solving. Multi-threaded Gurobi
+        %can non-deterministically report the MILP as infeasible (issue #607).
+        %Override by setting 'Threads' in the step's MILPParams (0 = all cores).
+        if ~isfield(params, 'Threads')
+            params.Threads = 1;
+        end
 
         if ~first 
             %There is sometimes a problem with that the objective function becomes close to zero,

INIT/ftINITFillGapsMILP.m

@@ -207,6 +207,7 @@
 params.TimeLimit = 300;
 params.Seed = 26;%This is weird - although it says "optimal solution found", we can get different results with different
                  %values of the objective function, where one is more optimal than the other (pretty big difference...)
+params.Threads = 1; %single-threaded Gurobi MILP for determinism, see issue #607
 %params.CSClientLog = 3;%generates a warning in gurobi, but may be of interest for other solvers
 
 % Optimize the problem

INIT/removeLowScoreGenes.m

@@ -109,13 +109,20 @@
 
 % regenerate "genes" and "rxnGeneMat" model fields
 [genes,rxnGeneMat] = getGenesFromGrRules(newModel.grRules);
-newModel.genes = genes;
-newModel.rxnGeneMat = rxnGeneMat;
 
-% update other gene-related fields
-remInd = ~ismember(model.genes,newModel.genes);
+% determine which of the original genes were removed
+remInd = ~ismember(model.genes,genes);
 remGenes = model.genes(remInd);
 
+% Keep the retained genes in their original order rather than the sorted
+% order returned by getGenesFromGrRules. Gene removal never introduces new
+% genes, so the remaining genes are model.genes(~remInd). Preserving this
+% order ensures the gene-associated fields trimmed below (geneShortNames,
+% proteins, etc.), which are indexed by remInd, stay aligned with genes.
+newModel.genes = model.genes(~remInd);
+[~,reorderInd] = ismember(newModel.genes,genes);
+newModel.rxnGeneMat = rxnGeneMat(:,reorderInd);
+
 if isfield(newModel,'geneShortNames')
     newModel.geneShortNames(remInd) = [];
 end

doc/INIT/INITStepDesc.html

@@ -70,55 +70,57 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
 0024                    <span class="comment">%      .mets        Names of metabolites to remove</span>
 0025                    <span class="comment">%      .compsToKeep Compartments for which metabolites should be kept.</span>
 0026       MILPParams <span class="comment">%Cell array of MILPparams - dictates how many iterations that will be run in this step.</span>
-0027                  <span class="comment">%Typically, MIPGap and TimeLimit is specified</span>
-0028       AbsMIPGaps <span class="comment">%If the objective is close to zero, a percentage of that is very small.</span>
-0029                  <span class="comment">%Therefore, also set an absolut value for this (typically 10 or 20).</span>
-0030                  <span class="comment">%For practical reasons, the first number is not used</span>
-0031    <span class="keyword">end</span>
-0032    methods
-0033       <a name="_sub0" href="#_subfunctions" class="code">function obj = INITStepDesc(posRevOff_, AllowMetSecr_, howToUsePrevResults_, rxnsToIgnoreMask_, metsToIgnore_, MILPParams_, absMIPGaps_)</a>
-0034          <span class="keyword">if</span> nargin &gt; 0
-0035             obj.PosRevOff = posRevOff_;
-0036          <span class="keyword">else</span>
-0037             obj.PosRevOff = false;
-0038          <span class="keyword">end</span>
-0039          <span class="keyword">if</span> nargin &gt; 1
-0040             obj.AllowMetSecr = AllowMetSecr_;
-0041          <span class="keyword">else</span>
-0042             obj.AllowMetSecr = false;
-0043          <span class="keyword">end</span>
-0044          <span class="keyword">if</span> nargin &gt; 2
-0045             obj.HowToUsePrevResults = howToUsePrevResults_;
-0046          <span class="keyword">else</span>
-0047             obj.HowToUsePrevResults = <span class="string">'essential'</span>;
-0048          <span class="keyword">end</span>
-0049          <span class="keyword">if</span> nargin &gt; 3
-0050             obj.RxnsToIgnoreMask = rxnsToIgnoreMask_;
-0051          <span class="keyword">else</span>
-0052             obj.RxnsToIgnoreMask = [1;0;0;0;0;0;0;0];
-0053          <span class="keyword">end</span>
-0054          <span class="keyword">if</span> nargin &gt; 4
-0055             obj.MetsToIgnore = metsToIgnore_;
-0056          <span class="keyword">else</span>
-0057             obj.MetsToIgnore = [1;0;0;0;0;0;0;0];
-0058          <span class="keyword">end</span>
-0059          <span class="keyword">if</span> nargin &gt; 5
-0060             obj.MILPParams = MILPParams_;
-0061          <span class="keyword">else</span>
-0062             params = struct();
-0063             params.TimeLimit = 5000;
-0064             params.MIPGap = 0.0004;            
-0065             obj.MILPParams = {params};
-0066          <span class="keyword">end</span>
-0067          
-0068          <span class="keyword">if</span> nargin &gt; 6
-0069             obj.AbsMIPGaps = absMIPGaps_;
-0070          <span class="keyword">else</span>
-0071             obj.AbsMIPGaps = 10;            
-0072          <span class="keyword">end</span>
-0073       <span class="keyword">end</span>
-0074    <span class="keyword">end</span>
-0075 <span class="keyword">end</span></pre></div>
+0027                  <span class="comment">%Typically, MIPGap and TimeLimit is specified. ftINIT additionally</span>
+0028                  <span class="comment">%defaults Threads to 1 (single-threaded Gurobi) for determinism;</span>
+0029                  <span class="comment">%set Threads here to override (0 = use all cores).</span>
+0030       AbsMIPGaps <span class="comment">%If the objective is close to zero, a percentage of that is very small.</span>
+0031                  <span class="comment">%Therefore, also set an absolut value for this (typically 10 or 20).</span>
+0032                  <span class="comment">%For practical reasons, the first number is not used</span>
+0033    <span class="keyword">end</span>
+0034    methods
+0035       <a name="_sub0" href="#_subfunctions" class="code">function obj = INITStepDesc(posRevOff_, AllowMetSecr_, howToUsePrevResults_, rxnsToIgnoreMask_, metsToIgnore_, MILPParams_, absMIPGaps_)</a>
+0036          <span class="keyword">if</span> nargin &gt; 0
+0037             obj.PosRevOff = posRevOff_;
+0038          <span class="keyword">else</span>
+0039             obj.PosRevOff = false;
+0040          <span class="keyword">end</span>
+0041          <span class="keyword">if</span> nargin &gt; 1
+0042             obj.AllowMetSecr = AllowMetSecr_;
+0043          <span class="keyword">else</span>
+0044             obj.AllowMetSecr = false;
+0045          <span class="keyword">end</span>
+0046          <span class="keyword">if</span> nargin &gt; 2
+0047             obj.HowToUsePrevResults = howToUsePrevResults_;
+0048          <span class="keyword">else</span>
+0049             obj.HowToUsePrevResults = <span class="string">'essential'</span>;
+0050          <span class="keyword">end</span>
+0051          <span class="keyword">if</span> nargin &gt; 3
+0052             obj.RxnsToIgnoreMask = rxnsToIgnoreMask_;
+0053          <span class="keyword">else</span>
+0054             obj.RxnsToIgnoreMask = [1;0;0;0;0;0;0;0];
+0055          <span class="keyword">end</span>
+0056          <span class="keyword">if</span> nargin &gt; 4
+0057             obj.MetsToIgnore = metsToIgnore_;
+0058          <span class="keyword">else</span>
+0059             obj.MetsToIgnore = [1;0;0;0;0;0;0;0];
+0060          <span class="keyword">end</span>
+0061          <span class="keyword">if</span> nargin &gt; 5
+0062             obj.MILPParams = MILPParams_;
+0063          <span class="keyword">else</span>
+0064             params = struct();
+0065             params.TimeLimit = 5000;
+0066             params.MIPGap = 0.0004;            
+0067             obj.MILPParams = {params};
+0068          <span class="keyword">end</span>
+0069          
+0070          <span class="keyword">if</span> nargin &gt; 6
+0071             obj.AbsMIPGaps = absMIPGaps_;
+0072          <span class="keyword">else</span>
+0073             obj.AbsMIPGaps = 10;            
+0074          <span class="keyword">end</span>
+0075       <span class="keyword">end</span>
+0076    <span class="keyword">end</span>
+0077 <span class="keyword">end</span></pre></div>
 <hr><address>Generated by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
 </body>
 </html>