feat: multiple small tweaks of functions

core/addMets.m

@@ -69,19 +69,21 @@
     metsToAdd.mets=generateNewIds(newModel,'mets',prefix,numel(metsToAdd.metNames));
 end
 if ~isfield(metsToAdd,'metNames')
-    EM='metNames is a required field in metsToAdd';
-    dispEM(EM);
+    metsToAdd.metNames=metsToAdd.mets;
 end
 if ~isfield(metsToAdd,'compartments')
     EM='compartments is a required field in metsToAdd';
     dispEM(EM);
 end
-
-if ~iscellstr(metsToAdd.mets)
+if ischar(metsToAdd.mets)
+    metsToAdd.mets={metsToAdd.mets};
+elseif ~iscellstr(metsToAdd.mets)
     EM='metsToAdd.mets must be a cell array of strings';
     dispEM(EM);
 end
-if ~iscellstr(metsToAdd.metNames)
+if ischar(metsToAdd.metNames)
+    metsToAdd.metNames={metsToAdd.metNames};
+elseif ~iscellstr(metsToAdd.metNames)
     EM='metsToAdd.metNames must be a cell array of strings';
     dispEM(EM);
 end

core/getIndexes.m

@@ -8,7 +8,7 @@
 %                   same number of elements as metabolites in the model,
 %                   of a vector of indexes
 %   type            'rxns', 'mets', or 'genes' depending on what to retrieve
-%                   'metnames' queries metabolite names, while 'metscomps'
+%                   'metnames' queries metabolite names, while 'metcomps'
 %                   allows to provide specific metabolites and their
 %                   compartments in the format metaboliteName[comp]
 %   returnLogical   Sets whether to return a logical array or an array with
@@ -32,42 +32,43 @@
 
 indexes=[];
 
-if strcmpi(type,'rxns')
-    searchIn=model.rxns;
-elseif strcmpi(type,'mets')
-    searchIn=model.mets;
-elseif strcmpi(type,'genes')
-    searchIn=model.genes;
-elseif strcmpi(type,'metnames')
-    searchIn=model.metNames;
-elseif strcmpi(type,'metscomps')
-    % If provided as metaboliteName[comp], then index search is quite
-    % different from general approach, and therefore coded separately.
-    if isstr(objects)
-        objects={objects};
-    end
-    mets=regexprep(objects,'(^.+)\[(.+)\]$','$1');
-    comps=regexprep(objects,'(^.+)\[(.+)\]$','$2');
-    for i=1:numel(objects)
-        index=find(strcmp(mets(i),model.metNames));
-        index=index(strcmp(comps(i),model.comps(model.metComps(index))));
-        if ~isempty(index)
-            indexes(i)=index;
+switch type
+    case 'rxns'
+        searchIn=model.rxns;
+    case 'mets'
+        searchIn=model.mets;
+    case 'genes'
+        searchIn=model.genes;
+    case 'metnames'
+        searchIn=model.metNames;
+    case 'metcomps'
+        % If provided as metaboliteName[comp], then index search is quite
+        % different from general approach, and therefore coded separately.
+        mets=regexprep(objects,'(^.+)\[(.+)\]$','$1');
+        comps=regexprep(objects,'(^.+)\[(.+)\]$','$2');
+        indexes=zeros(1,numel(objects));
+        for i=1:numel(objects)
+            index=find(strcmp(mets(i),model.metNames));
+            index=index(strcmp(comps(i),model.comps(model.metComps(index))));
+            if ~isempty(index)
+                indexes(i)=index;
+            else
+                error(['Could not find object ' objects{i} ' in the model']);
+            end
+        end
+        indexes=indexes(:);
+        if returnLogical==true
+            tempIndexes=false(numel(model.mets),1);
+            tempIndexes(indexes)=true;
+            indexes=tempIndexes;
+        end
+        return %None of the remaining function needs to run if metcomps
+    otherwise
+        if contains(objects,{'rxns','mets','metnames','metcomps','genes'})
+            error('The second and third parameter provided to run getIndexes are likely switched. Note that the second parameter is the object to find the index for, while the third parameter is the object type (''rxns'', ''mets'', ''metnames'', ''metcomps'' or ''genes'').')
         else
-            EM=['Could not find object ' objects{i} ' in the model'];
-            dispEM(EM);
+            error('Incorrect value of the ''type'' parameter. Allowed values are ''rxns'', ''mets'', ''metnames'', ''metcomps'' or ''genes''.');
         end
-    end
-    indexes=indexes(:);
-    if returnLogical==true
-        tempIndexes=false(numel(model.mets),1);
-        tempIndexes(indexes)=true;
-        indexes=tempIndexes;
-    end
-    return; % If metscomps is queried, remaining codes doesn't need executing
-else
-    EM='Incorrect value of the "type" parameter. Allowed values are "rxns", "mets" or "genes"';
-    dispEM(EM);
 end
 
 if iscell(objects)
@@ -78,8 +79,7 @@
         elseif ~isempty(index)
             indexes(i)=index;
         else
-            EM=['Could not find object ' objects{i} ' in the model'];
-            dispEM(EM);
+            error(['Could not find object ' objects{i} ' in the model']);
         end
     end
 else

core/randomSampling.m

@@ -1,4 +1,4 @@
-function [solutions, goodRxns]=randomSampling(model, nSamples, replaceBoundsWithInf,supressErrors,showProgress,goodRxns)
+function [solutions, goodRxns]=randomSampling(model,nSamples,replaceBoundsWithInf,supressErrors,showProgress,goodRxns,minFlux)
 % randomSampling
 %   Returns a number of random solutions
 %
@@ -27,6 +27,12 @@
 %                           as random objective functions, as generated by
 %                           a previous run of randomSampling on the same
 %                           model (opt, default empty)
+%   minFlux                 determines if a second optimization should be
+%                           performed for each random sample, to minimize
+%                           the number of fluxes and thereby preventing
+%                           loops. Typically, loops are averaged out when a
+%                           large number of samples are taken, but this is
+%                           not always the case (opt, default false)
 %
 %   solutions               matrix with the solutions
 %   goodRxns                double vector of indexes of those reactions
@@ -37,7 +43,7 @@
 %   random set of three reactions. For reversible reactions it randomly
 %   chooses between maximizing and minimizing.
 %
-%   Usage: solutions=randomSampling(model, nSamples, replaceBoundsWithInf)
+%   Usage: solutions=randomSampling(model,nSamples,replaceBoundsWithInf,supressErrors,showProgress,goodRxns,minFlux)
 
 if nargin<2
     nSamples=1000;
@@ -51,6 +57,9 @@
 if nargin<5
     showProgress=false;
 end
+if nargin<6
+    minFlux=false;
+end
 
 nRxns=2; %Number of reactions in the objective function in each iteration
 
@@ -76,7 +85,7 @@
 
 %Reactions which can be involved in loops should not be optimized for.
 %Check which reactions reach an arbitary high upper bound
-if nargin<6
+if nargin<6 | isempty(goodRxns)
     goodRxns=true(numel(model.rxns),1);
     for i=1:numel(model.rxns)
         if showProgress && rem(i,100) == 0
@@ -117,7 +126,11 @@
     multipliers=randsample([-1 1],nRxns,true);
     multipliers(model.rev(goodRxns(rxns))==0)=1;
     model.c(goodRxns(rxns))=rand(nRxns,1).*multipliers;
-    sol=solveLP(model);
+    if true(minFlux)
+        sol=solveLP(model,1);
+    else
+        sol=solveLP(model);
+    end
     if any(sol.x)
         if abs(sol.f)>10^-8
             sols(:,counter)=sol.x;