fix: various minor fixes

core/getAllowedBounds.m

@@ -3,25 +3,30 @@
 %   Returns the minimal and maximal fluxes through each reaction.
 %
 % Input:
-%   model       a model structure
-%   rxns        either a cell array of reaction IDs, a logical vector with the
-%               same number of elements as reactions in the model, or a vector
-%               of reaction indexes (opt, default model.rxns)
-%   runParallel make use of MATLAB parallel pool to speed up calculations. Not
-%               beneficial if only a limited number of reactions are simulated.
-%               (opt, default true)
+%   model           a model structure
+%   rxns            either a cell array of reaction IDs, a logical vector
+%                   with the same number of elements as reactions in the
+%                   model, or a vector of reaction indexes (opt, default
+%                   model.rxns)
+%   runParallel     speed up calculations by parallel processing. This is
+%                   not beneficial if allowed bounds are calculated for
+%                   only a few reactions, as the overhead of parallel
+%                   processing will take longer. It requires MATLAB
+%                   Parallel Computing Toolbox. If this is not installed,
+%                   the calculations will not be parallelized, regardless
+%                   what is indicated as runParallel. (opt, default true)
 %
 % Output:
-%   minFluxes   minimal allowed fluxes
-%   maxFluxes   maximal allowed fluxes
-%   exitFlags   exit flags for min/max for each of the reactions. True if it was
-%               possible to calculate a flux
+%   minFluxes       minimal allowed fluxes
+%   maxFluxes       maximal allowed fluxes
+%   exitFlags       exit flags for min/max for each of the reactions. True
+%                   if it was possible to calculate a flux
 %
 %   NOTE: In cases where no solution can be calculated, NaN is returned.
 %
 % Usage: [minFluxes, maxFluxes, exitFlags] = getAllowedBounds(model, rxns, runParallel)
 
-if nargin<2
+if nargin<2 || isempty(rxns)
     rxns = 1:numel(model.rxns);
 elseif ~islogical(rxns) && ~isnumeric(rxns)
     rxns = convertCharArray(rxns);
@@ -30,78 +35,39 @@
 if nargin<3
     runParallel = true;
 end
-if runParallel
-    addonList = matlab.addons.installedAddons;
-    if ~any(strcmpi(addonList.Name,'Parallel Computing Toolbox'))
-        disp('Cannot find MATLAB Parallel Computing Toolbox, process is not parallelized.')
-        runParallel = false;
-    end
-end
+
+[ps, oldPoolAutoCreateSetting] = parallelPoolRAVEN(runParallel);
 
 minFluxes = zeros(numel(rxns),1);
 maxFluxes = zeros(numel(rxns),1);
 exitFlags = zeros(numel(rxns),2);
 c = zeros(numel(model.rxns),1);
 
-p = 1;
-progressbar('Calculating the minimal and maximal fluxes')
-if runParallel
-    D = parallel.pool.DataQueue;
-    afterEach(D, @updateProgressParallel);
-    parfor i = 1:numel(rxns)
-        tmpModel = model;
-        tmpModel.c = c;
-
-        % Get minimal flux
-        tmpModel.c(rxns(i)) = -1;
-        solMin = solveLP(tmpModel);
-        if ~isempty(solMin.f)
-            minFluxes(i) = solMin.x(rxns(i));
-        else
-            minFluxes(i) = NaN;
-        end
+PB = ProgressBar2(numel(rxns),'Running getAllowedBounds','cli');
+parfor i = 1:numel(rxns)
+    count(PB)
+    tmpModel = model;
+    tmpModel.c = c;
 
-        % 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
-        send(D, i);
+    % Get minimal flux
+    tmpModel.c(rxns(i)) = -1;
+    solMin = solveLP(tmpModel);
+    if ~isempty(solMin.f)
+        minFluxes(i) = solMin.x(rxns(i));
+    else
+        minFluxes(i) = NaN;
     end
-else
-    for i = 1:numel(rxns)
-        tmpModel = model;
-        tmpModel.c = c;
 
-        % Get minimal flux
-        tmpModel.c(rxns(i)) = -1;
-        solMin=solveLP(tmpModel);
-        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
-        progressbar(i/numel(rxns))
+    % 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
-progressbar(1) % Make sure it closes
-
-    function updateProgressParallel(~)
-        progressbar(p/numel(rxns))
-        p = p + 1;
-    end
+% Reset original Parallel setting
+ps.Pool.AutoCreate = oldPoolAutoCreateSetting;
 end

core/getGenesFromGrRules.m

@@ -31,7 +31,7 @@
 end
 
 % check if the grRules use written or symbolic boolean operators
-if any(contains(grRules,{'&','|'}))
+if any(contains(grRules,{' & ',' | '}))
     % fix some potential missing spaces between parentheses and &/|
     grRules = regexprep(grRules,'\)&',') &');   % ")&"  ->  ") &"
     grRules = regexprep(grRules,'&\(','& (');   % "&("  ->  "& ("
@@ -50,11 +50,12 @@
 end
 
 % extract list of genes from each reaction
-rxnGenes = cellfun(@(r) unique(regexp(r,'[^&|\(\) ]+','match')),grRules,'UniformOutput',false);
+rxnGenes = cellfun(@(r) regexprep(unique(strsplit(r,{' | ',' & '})),'[\(\) ]+',''),grRules,'UniformOutput',false);
 
 % construct new gene list
 nonEmpty = ~cellfun(@isempty,rxnGenes);
 genes = unique([rxnGenes{nonEmpty}]');
+genes(cellfun(@isempty,genes)) = [];
 
 if ~isempty(originalGenes)
     if ~isequal(sort(originalGenes), sort(genes))

core/getMinNrFluxes.m

@@ -135,7 +135,7 @@
 prob=rmfield(prob,{'blx','bux','blc','buc'});
 
 % Optimize the problem
-res = optimizeProb(prob,params);
+res = optimizeProb(prob,params,false);
 isFeasible=checkSolution(res);
 
 if ~isFeasible

core/getModelFromHomology.m

@@ -115,6 +115,11 @@
     if isfield(models{i},'geneMiriams')
         models{i}=rmfield(models{i},'geneMiriams');
     end
+    %The geneFrom field also loses meaning if the genes are replaced by
+    %orthologs
+    if isfield(models{i},'geneFrom')
+        models{i}=rmfield(models{i},'geneFrom');
+    end
 end
 
 %Check that genes do not begin with ( or end with ), as this makes problematic grRules

core/mapCompartments.m

@@ -69,7 +69,7 @@
 %Check that there are no compartments in the rules that are not in the
 %geneScoreStructure.
 uComps=upper(geneScoreStructure.compartments);
-J=[uComps;{'OTHER'}];
+J=[uComps,{'OTHER'}];
 
 if ~isempty(setdiff([toKeep;toMerge],J))
     EM='There are compartment in the rules that are not in geneScoreStructure.compartments';

core/parallelPoolRAVEN.m

@@ -0,0 +1,53 @@
+function [ps, oldPoolAutoCreate] = parallelPoolRAVEN(runParallel)
+% handleParallelRAVEN
+%   Called by RAVEN functions that support parallel processing, to confirm
+%   whether the MATLAB Parallel Computing Toolbox is installed.
+%   - The toolbox is installed, and runParallel == true, ==> a parallel
+%     pool is started.
+%   - The toolbox is installed, but runParallel == false, ==> the auto-
+%     creation of a parallel pool is disabled, to prevent "parfor" in
+%     the target function to start a pool anyway.
+%   - The toolbox is not installed, and runParallel == true, ==> a warning
+%     is displayed that parallel computer is not possible.
+%   - The toolbox is not installed, and runParallel == false, ==> the
+%     target runs as intended, as "parfor" will automatically run in serial
+%     mode instead.
+%
+% Input:
+%   runParallel         logical, whether the target function (which calls
+%                       parallelPoolRAVEN) should be run in parallel (opt,
+%                       default true)
+%
+% Output:
+%   ps                  parallel settings structure that will be used by
+%                       the target function    
+%   oldPoolAutoCreate   logical, to reset the original ps.Pool.AutoCreate
+%                       setting once the target function has finished
+%
+% Use: [ps, oldPoolAutoCreate] = parallelPoolRAVEN(runParallel)
+
+if nargin<1 || isempty(runParallel)
+    runParallel = true;
+end
+
+addonList = matlab.addons.installedAddons;
+ps = []; oldPoolAutoCreate = [];
+if ~any(strcmpi(addonList.Name,'Parallel Computing Toolbox'))
+    if runParallel % User wants parallel, but will not be possible
+        disp('Cannot find MATLAB Parallel Computing Toolbox, process is not parallelized.')
+    end
+else
+    if ~runParallel % User has Parallel toolbox, but does not want pool to start.
+        % If pool is already running, then it will use the max workers
+        % anyway, but this is probably okay.
+        ps = parallel.Settings;
+        oldPoolAutoCreate = ps.Pool.AutoCreate;
+        ps.Pool.AutoCreate = false;
+    else
+        pool = gcp('nocreate');
+        if isempty(pool)
+            parpool('IdleTimeout',120)
+        end
+    end
+end
+end

core/predictLocalization.m

@@ -1005,3 +1005,122 @@
 transportCost=sum(transportCost(I));
 score=geneScore-transportCost;
 end
+
+% To avoid dependency on stats toolbox, use this alternative implementation
+% of randsample, source:
+% https://github.com/gpeyre/numerical-tours/blob/dacee30081c04ef5f67b26b387ead85f2b193af9/matlab/toolbox_signal/randsample.m
+function y = randsample(n, k, replace, w)
+%RANDSAMPLE Random sample, with or without replacement.
+%   Y = RANDSAMPLE(N,K) returns Y as a vector of K values sampled uniformly
+%   at random, without replacement, from the integers 1:N.
+%
+%   Y = RANDSAMPLE(POPULATION,K) returns K values sampled uniformly at
+%   random, without replacement, from the values in the vector POPULATION.
+%
+%   Y = RANDSAMPLE(...,REPLACE) returns a sample taken with replacement if
+%   REPLACE is true, or without replacement if REPLACE is false (the default).
+%
+%   Y = RANDSAMPLE(...,true,W) returns a weighted sample, using positive
+%   weights W, taken with replacement.  W is often a vector of probabilities.
+%   This function does not support weighted sampling without replacement.
+%
+%   Example:  Generate a random sequence of the characters ACGT, with
+%   replacement, according to specified probabilities.
+%
+%      R = randsample('ACGT',48,true,[0.15 0.35 0.35 0.15])
+%
+%   See also RAND, RANDPERM.
+
+%   Copyright 1993-2008 The MathWorks, Inc.
+%   $Revision: 1.1.4.3 $  $Date: 2008/12/01 08:09:34 $
+
+if nargin < 2
+    error('stats:randsample:TooFewInputs','Requires two input arguments.');
+elseif numel(n) == 1
+    population = [];
+else
+    population = n;
+    n = numel(population);
+    if length(population)~=n
+       error('stats:randsample:BadPopulation','POPULATION must be a vector.');
+    end
+end
+
+if nargin < 3
+    replace = false;
+end
+
+if nargin < 4
+    w = [];
+elseif ~isempty(w)
+    if length(w) ~= n
+        if isempty(population)
+            error('stats:randsample:InputSizeMismatch',...
+                  'W must have length equal to N.');
+        else
+            error('stats:randsample:InputSizeMismatch',...
+                  'W must have the same length as the population.');
+        end
+    else
+        p = w(:)' / sum(w);
+    end
+end
+
+switch replace
+
+% Sample with replacement
+case {true, 'true', 1}
+    if isempty(w)
+        y = ceil(n .* rand(k,1));
+    else
+        [dum, y] = histc(rand(k,1),[0 cumsum(p)]);
+    end
+
+% Sample without replacement
+case {false, 'false', 0}
+    if k > n
+        if isempty(population)
+            error('stats:randsample:SampleTooLarge',...
+        'K must be less than or equal to N for sampling without replacement.');
+        else
+            error('stats:randsample:SampleTooLarge',...
+                  'K must be less than or equal to the population size.');
+        end
+    end
+
+    if isempty(w)
+        % If the sample is a sizeable fraction of the population,
+        % just randomize the whole population (which involves a full
+        % sort of n random values), and take the first k.
+        if 4*k > n
+            rp = randperm(n);
+            y = rp(1:k);
+
+        % If the sample is a small fraction of the population, a full sort
+        % is wasteful.  Repeatedly sample with replacement until there are
+        % k unique values.
+        else
+            x = zeros(1,n); % flags
+            sumx = 0;
+            while sumx < k
+                x(ceil(n * rand(1,k-sumx))) = 1; % sample w/replacement
+                sumx = sum(x); % count how many unique elements so far
+            end
+            y = find(x > 0);
+            y = y(randperm(k));
+        end
+    else
+        error('stats:randsample:NoWeighting',...
+              'Weighted sampling without replacement is not supported.');
+    end
+otherwise
+    error('stats:randsample:BadReplaceValue',...
+          'REPLACE must be either true or false.');
+end
+
+if ~isempty(population)
+    y = population(y);
+else
+    y = y(:);
+end
+end

core/printFluxes.m

@@ -63,7 +63,10 @@ function printFluxes(model, fluxes, onlyExchange, cutOffFlux, outputFile,outputS
 else
     metaboliteList=convertCharArray(metaboliteList);
 end
-if size(fluxes,1)~=numel(model.rxns)
+if isempty(fluxes)
+    EM='Empty vector of fluxes, solveLP possibly returned infeasible';
+    dispEM(EM);
+elseif size(fluxes,1)~=numel(model.rxns)
     EM='The number of fluxes and the number of reactions must be the same';
     dispEM(EM);
 end