Skip to content

Commit

Permalink
fix: parallelization (#541)
Browse files Browse the repository at this point in the history 
* fix: code tweaks to allow for parallelization

* fix: exportModel valid libSBML output

* revert: haveFlux faster if not parallelized

* fix: solverTests

* fix: tinitTests check if gurobi is functional

* fix: fillGaps tests check if Gurobi is functional

* fix: rename SBML binaries to ensure being used

* fix: checkInstallation new libSBML filenames

* fix: exportForGit SBML binaries filenames
  • Loading branch information
@edkerk
edkerk committed May 26, 2024
1 parent f92f6fa commit d95e03e
Show file tree
Hide file tree
Showing 37 changed files with 2,907 additions and 2,893 deletions.
21 changes: 11 additions & 10 deletions core/getAllowedBounds.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,24 +49,25 @@
tmpModel = model;
tmpModel.c = c;

% Get maximal flux
tmpModel.c(rxns(i)) = 1;
[solMax,hsSol]=solveLP(tmpModel);
if ~isempty(solMax.f)
maxFluxes(i) = solMax.x(rxns(i));
else
maxFluxes(i) = NaN;
end

% Get minimal flux
tmpModel.c(rxns(i)) = -1;
solMin = solveLP(tmpModel);
solMin = solveLP(tmpModel,[],[],hsSol);
if ~isempty(solMin.f)
minFluxes(i) = solMin.x(rxns(i));
else
minFluxes(i) = NaN;
end

% Get maximal flux
tmpModel.c(rxns(i)) = 1;
solMax=solveLP(tmpModel);
exitFlags(i,:) = [solMin.stat solMax.stat];
if ~isempty(solMax.f)
maxFluxes(i) = solMax.x(rxns(i));
else
maxFluxes(i) = NaN;
end

end
% Reset original Parallel setting
ps.Pool.AutoCreate = oldPoolAutoCreateSetting;
Expand Down
2 changes: 1 addition & 1 deletion core/getGenesFromGrRules.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,7 +54,7 @@

% construct new gene list
nonEmpty = ~cellfun(@isempty,rxnGenes);
genes = unique([rxnGenes{nonEmpty}]');
genes = unique(transpose([rxnGenes{nonEmpty}]));
genes(cellfun(@isempty,genes)) = [];

if ~isempty(originalGenes)
Expand Down
20 changes: 11 additions & 9 deletions core/haveFlux.m
View file
Original file line number Diff line number Diff line change
@@ -1,23 +1,25 @@
function I=haveFlux(model,cutOff,rxns)
% haveFlux
% Checks which reactions can carry a (positive or negative) flux.
% Is used as a faster version of getAllowedBounds if it is only interesting
% Checks which reactions can carry a (positive or negative) flux. Is used
% as a faster version of getAllowedBounds if it is only interesting
% whether the reactions can carry a flux or not
%
% Input:
% model a model structure
% cutOff the flux value that a reaction has to carry to be
% identified as positive (optional, default 10^-8)
% rxns either a cell array of IDs, a logical vector with the
% same number of elements as metabolites in the model,
% of a vector of indexes (optional, default model.rxns)
% same number of elements as metabolites in the model, or a
% vector of indexes (optional, default model.rxns)
%
% I logical array with true if the corresponding
% reaction can carry a flux
% Output:
% I logical array with true if the corresponding reaction can
% carry a flux
%
% NOTE: If a model has +/- Inf bounds then those are replaced with an
% arbitary large value of +/- 10000 prior to solving
% If a model has +/- Inf bounds then those are replaced with an arbitary
% large value of +/- 10000 prior to solving
%
% Usage: I=haveFlux(model,cutOff, rxns)
% Usage: I = haveFlux(model, cutOff, rxns)

if nargin<2
cutOff=10^-6;
Expand Down
39 changes: 20 additions & 19 deletions doc/core/getAllowedBounds.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,28 +118,29 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0049 tmpModel = model;
0050 tmpModel.c = c;
0051
0052 <span class="comment">% Get minimal flux</span>
0053 tmpModel.c(rxns(i)) = -1;
0054 solMin = solveLP(tmpModel);
0055 <span class="keyword">if</span> ~isempty(solMin.f)
0056 minFluxes(i) = solMin.x(rxns(i));
0052 <span class="comment">% Get maximal flux</span>
0053 tmpModel.c(rxns(i)) = 1;
0054 [solMax,hsSol]=solveLP(tmpModel);
0055 <span class="keyword">if</span> ~isempty(solMax.f)
0056 maxFluxes(i) = solMax.x(rxns(i));
0057 <span class="keyword">else</span>
0058 minFluxes(i) = NaN;
0058 maxFluxes(i) = NaN;
0059 <span class="keyword">end</span>
0060
0061 <span class="comment">% Get maximal flux</span>
0062 tmpModel.c(rxns(i)) = 1;
0063 solMax=solveLP(tmpModel);
0064 exitFlags(i,:) = [solMin.stat solMax.stat];
0065 <span class="keyword">if</span> ~isempty(solMax.f)
0066 maxFluxes(i) = solMax.x(rxns(i));
0067 <span class="keyword">else</span>
0068 maxFluxes(i) = NaN;
0069 <span class="keyword">end</span>
0070 <span class="keyword">end</span>
0071 <span class="comment">% Reset original Parallel setting</span>
0072 ps.Pool.AutoCreate = oldPoolAutoCreateSetting;
0073 <span class="keyword">end</span></pre></div>
0061 <span class="comment">% Get minimal flux</span>
0062 tmpModel.c(rxns(i)) = -1;
0063 solMin = solveLP(tmpModel,[],[],hsSol);
0064 <span class="keyword">if</span> ~isempty(solMin.f)
0065 minFluxes(i) = solMin.x(rxns(i));
0066 <span class="keyword">else</span>
0067 minFluxes(i) = NaN;
0068 <span class="keyword">end</span>
0069 exitFlags(i,:) = [solMin.stat solMax.stat];
0070
0071 <span class="keyword">end</span>
0072 <span class="comment">% Reset original Parallel setting</span>
0073 ps.Pool.AutoCreate = oldPoolAutoCreateSetting;
0074 <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>
2 changes: 1 addition & 1 deletion doc/core/getGenesFromGrRules.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0054
0055 <span class="comment">% construct new gene list</span>
0056 nonEmpty = ~cellfun(@isempty,rxnGenes);
0057 genes = unique([rxnGenes{nonEmpty}]');
0057 genes = unique(transpose([rxnGenes{nonEmpty}]));
0058 genes(cellfun(@isempty,genes)) = [];
0059
0060 <span class="keyword">if</span> ~isempty(originalGenes)
Expand Down
172 changes: 88 additions & 84 deletions doc/core/haveFlux.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,24 +28,26 @@ <h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre class="comment"> haveFlux
Checks which reactions can carry a (positive or negative) flux.
Is used as a faster version of getAllowedBounds if it is only interesting
Checks which reactions can carry a (positive or negative) flux. Is used
as a faster version of getAllowedBounds if it is only interesting
whether the reactions can carry a flux or not

Input:
model a model structure
cutOff the flux value that a reaction has to carry to be
identified as positive (optional, default 10^-8)
rxns either a cell array of IDs, a logical vector with the
same number of elements as metabolites in the model,
of a vector of indexes (optional, default model.rxns)
same number of elements as metabolites in the model, or a
vector of indexes (optional, default model.rxns)

I logical array with true if the corresponding
reaction can carry a flux
Output:
I logical array with true if the corresponding reaction can
carry a flux

NOTE: If a model has +/- Inf bounds then those are replaced with an
arbitary large value of +/- 10000 prior to solving
If a model has +/- Inf bounds then those are replaced with an arbitary
large value of +/- 10000 prior to solving

Usage: I=haveFlux(model,cutOff, rxns)</pre></div>
Usage: I = haveFlux(model, cutOff, rxns)</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
Expand All @@ -62,83 +64,85 @@ <h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^
<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function I=haveFlux(model,cutOff,rxns)</a>
0002 <span class="comment">% haveFlux</span>
0003 <span class="comment">% Checks which reactions can carry a (positive or negative) flux.</span>
0004 <span class="comment">% Is used as a faster version of getAllowedBounds if it is only interesting</span>
0003 <span class="comment">% Checks which reactions can carry a (positive or negative) flux. Is used</span>
0004 <span class="comment">% as a faster version of getAllowedBounds if it is only interesting</span>
0005 <span class="comment">% whether the reactions can carry a flux or not</span>
0006 <span class="comment">%</span>
0007 <span class="comment">% model a model structure</span>
0008 <span class="comment">% cutOff the flux value that a reaction has to carry to be</span>
0009 <span class="comment">% identified as positive (optional, default 10^-8)</span>
0010 <span class="comment">% rxns either a cell array of IDs, a logical vector with the</span>
0011 <span class="comment">% same number of elements as metabolites in the model,</span>
0012 <span class="comment">% of a vector of indexes (optional, default model.rxns)</span>
0013 <span class="comment">%</span>
0014 <span class="comment">% I logical array with true if the corresponding</span>
0015 <span class="comment">% reaction can carry a flux</span>
0016 <span class="comment">%</span>
0017 <span class="comment">% NOTE: If a model has +/- Inf bounds then those are replaced with an</span>
0018 <span class="comment">% arbitary large value of +/- 10000 prior to solving</span>
0019 <span class="comment">%</span>
0020 <span class="comment">% Usage: I=haveFlux(model,cutOff, rxns)</span>
0021
0022 <span class="keyword">if</span> nargin&lt;2
0023 cutOff=10^-6;
0024 <span class="keyword">end</span>
0025 <span class="keyword">if</span> isempty(cutOff)
0026 cutOff=10^-6;
0027 <span class="keyword">end</span>
0028 <span class="keyword">if</span> nargin&lt;3
0029 rxns=model.rxns;
0030 <span class="keyword">elseif</span> ~islogical(rxns) &amp;&amp; ~isnumeric(rxns)
0031 rxns=<a href="convertCharArray.html" class="code" title="function inputConverted = convertCharArray(funcInput)">convertCharArray</a>(rxns);
0032 <span class="keyword">end</span>
0033
0034 <span class="comment">%This is since we're maximizing for the sum of fluxes, which isn't possible</span>
0035 <span class="comment">%when there are infinite bounds</span>
0036 model.lb(model.lb==-inf)=-10000;
0037 model.ub(model.ub==inf)=10000;
0038
0039 <span class="comment">%Get the reaction IDs. A bit of an awkward way, but fine.</span>
0040 indexes=<a href="getIndexes.html" class="code" title="function indexes=getIndexes(model, objects, type, returnLogical)">getIndexes</a>(model,rxns,<span class="string">'rxns'</span>);
0041 rxns=model.rxns(indexes);
0042
0043 <span class="comment">%First remove all dead-end reactions</span>
0044 smallModel=<a href="simplifyModel.html" class="code" title="function [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,deleteUnconstrained, deleteDuplicates, deleteZeroInterval, deleteInaccessible, deleteMinMax, groupLinear, constrainReversible, reservedRxns, suppressWarnings)">simplifyModel</a>(model,false,false,true,true);
0045
0046 <span class="comment">%Get the indexes of the reactions in the connected model</span>
0047 indexes=<a href="getIndexes.html" class="code" title="function indexes=getIndexes(model, objects, type, returnLogical)">getIndexes</a>(smallModel,intersect(smallModel.rxns,rxns),<span class="string">'rxns'</span>);
0048 J=false(numel(indexes),1);
0049 mixIndexes=indexes(randperm(numel(indexes)));
0050
0051 <span class="comment">%Maximize for all fluxes first in order to get fewer rxns to test</span>
0052 smallModel.c=ones(numel(smallModel.c),1);
0053 sol=solveLP(smallModel);
0054 <span class="keyword">if</span> ~isempty(sol.x)
0055 J(abs(sol.x(mixIndexes))&gt;cutOff)=true;
0056 <span class="keyword">end</span>
0057
0058 <span class="comment">%Loop through and maximize then minimize each rxn if it does not already</span>
0059 <span class="comment">%have a flux</span>
0060 Z=zeros(numel(smallModel.c),1);
0061 hsSolOut=[];
0062 <span class="keyword">for</span> i=[1 -1]
0063 <span class="keyword">for</span> j=1:numel(J)
0064 <span class="keyword">if</span> J(j)==false
0065 <span class="comment">%Only minimize for reversible reactions</span>
0066 <span class="keyword">if</span> i==1 || smallModel.rev(mixIndexes(j))~=0
0067 smallModel.c=Z;
0068 smallModel.c(mixIndexes(j))=i;
0069 [sol, hsSolOut]=solveLP(smallModel,0,[],hsSolOut);
0070 <span class="keyword">if</span> any(sol.x)
0071 J(abs(sol.x(mixIndexes))&gt;cutOff)=true;
0072 <span class="keyword">end</span>
0073 <span class="keyword">end</span>
0074 <span class="keyword">end</span>
0075 <span class="keyword">end</span>
0076 <span class="keyword">end</span>
0077
0078 I=ismember(rxns,smallModel.rxns(mixIndexes(J)));
0079 <span class="keyword">end</span></pre></div>
0007 <span class="comment">% Input:</span>
0008 <span class="comment">% model a model structure</span>
0009 <span class="comment">% cutOff the flux value that a reaction has to carry to be</span>
0010 <span class="comment">% identified as positive (optional, default 10^-8)</span>
0011 <span class="comment">% rxns either a cell array of IDs, a logical vector with the</span>
0012 <span class="comment">% same number of elements as metabolites in the model, or a</span>
0013 <span class="comment">% vector of indexes (optional, default model.rxns)</span>
0014 <span class="comment">%</span>
0015 <span class="comment">% Output:</span>
0016 <span class="comment">% I logical array with true if the corresponding reaction can</span>
0017 <span class="comment">% carry a flux</span>
0018 <span class="comment">%</span>
0019 <span class="comment">% If a model has +/- Inf bounds then those are replaced with an arbitary</span>
0020 <span class="comment">% large value of +/- 10000 prior to solving</span>
0021 <span class="comment">%</span>
0022 <span class="comment">% Usage: I = haveFlux(model, cutOff, rxns)</span>
0023
0024 <span class="keyword">if</span> nargin&lt;2
0025 cutOff=10^-6;
0026 <span class="keyword">end</span>
0027 <span class="keyword">if</span> isempty(cutOff)
0028 cutOff=10^-6;
0029 <span class="keyword">end</span>
0030 <span class="keyword">if</span> nargin&lt;3
0031 rxns=model.rxns;
0032 <span class="keyword">elseif</span> ~islogical(rxns) &amp;&amp; ~isnumeric(rxns)
0033 rxns=<a href="convertCharArray.html" class="code" title="function inputConverted = convertCharArray(funcInput)">convertCharArray</a>(rxns);
0034 <span class="keyword">end</span>
0035
0036 <span class="comment">%This is since we're maximizing for the sum of fluxes, which isn't possible</span>
0037 <span class="comment">%when there are infinite bounds</span>
0038 model.lb(model.lb==-inf)=-10000;
0039 model.ub(model.ub==inf)=10000;
0040
0041 <span class="comment">%Get the reaction IDs. A bit of an awkward way, but fine.</span>
0042 indexes=<a href="getIndexes.html" class="code" title="function indexes=getIndexes(model, objects, type, returnLogical)">getIndexes</a>(model,rxns,<span class="string">'rxns'</span>);
0043 rxns=model.rxns(indexes);
0044
0045 <span class="comment">%First remove all dead-end reactions</span>
0046 smallModel=<a href="simplifyModel.html" class="code" title="function [reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,deleteUnconstrained, deleteDuplicates, deleteZeroInterval, deleteInaccessible, deleteMinMax, groupLinear, constrainReversible, reservedRxns, suppressWarnings)">simplifyModel</a>(model,false,false,true,true);
0047
0048 <span class="comment">%Get the indexes of the reactions in the connected model</span>
0049 indexes=<a href="getIndexes.html" class="code" title="function indexes=getIndexes(model, objects, type, returnLogical)">getIndexes</a>(smallModel,intersect(smallModel.rxns,rxns),<span class="string">'rxns'</span>);
0050 J=false(numel(indexes),1);
0051 mixIndexes=indexes(randperm(numel(indexes)));
0052
0053 <span class="comment">%Maximize for all fluxes first in order to get fewer rxns to test</span>
0054 smallModel.c=ones(numel(smallModel.c),1);
0055 sol=solveLP(smallModel);
0056 <span class="keyword">if</span> ~isempty(sol.x)
0057 J(abs(sol.x(mixIndexes))&gt;cutOff)=true;
0058 <span class="keyword">end</span>
0059
0060 <span class="comment">%Loop through and maximize then minimize each rxn if it does not already</span>
0061 <span class="comment">%have a flux</span>
0062 Z=zeros(numel(smallModel.c),1);
0063 hsSolOut=[];
0064 <span class="keyword">for</span> i=[1 -1]
0065 <span class="keyword">for</span> j=1:numel(J)
0066 <span class="keyword">if</span> J(j)==false
0067 <span class="comment">%Only minimize for reversible reactions</span>
0068 <span class="keyword">if</span> i==1 || smallModel.rev(mixIndexes(j))~=0
0069 smallModel.c=Z;
0070 smallModel.c(mixIndexes(j))=i;
0071 [sol, hsSolOut]=solveLP(smallModel,0,[],hsSolOut);
0072 <span class="keyword">if</span> any(sol.x)
0073 J(abs(sol.x(mixIndexes))&gt;cutOff)=true;
0074 <span class="keyword">end</span>
0075 <span class="keyword">end</span>
0076 <span class="keyword">end</span>
0077 <span class="keyword">end</span>
0078 <span class="keyword">end</span>
0079
0080 I=ismember(rxns,smallModel.rxns(mixIndexes(J)));
0081 <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>

0 comments on commit d95e03e

Please sign in to comment.