Home > RAVEN > randomSampling.m

randomSampling

PURPOSE ^

randomSampling

SYNOPSIS ^

function solutions=randomSampling(model, nSamples, replaceBoundsWithInf,supressErrors)

DESCRIPTION ^

 randomSampling
   Returns a number of random solutions

   model                   a model structure
   nSamples                the number of solutions to return
                           (opt, default 1000)
   replaceBoundsWithInf    replace the largest upper bounds with Inf and
                           the smallest lower bounds with -Inf. This is
                           needed in order to get solutions without loops
                           if your model has for example 1000/-1000 as 
                           arbitary large bounds. If your model only has
                           "biologically relevant" bounds, then set this
                           to false (opt, default true)
   supressErrors           the program will halt if it has problems
                           finding non-zero solutions which are not
                           involved in loops. This could be because the
                           constraints on the model are too relaxed (such
                           as unlimited glucose uptake) or too strict
                           (such as too many and too narrow constraints)
                           (opt, default false)

   solutions               matrix with the solutions

   The solutions are generated by maximizing (with random weights) for a
   random set of three reactions. For reversible reactions it randomly
   chooses between maximizing and minimizing.

   Usage: solutions=randomSampling(model, nSamples, replaceBoundsWithInf)

   Rasmus Agren, 2013-04-16

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function solutions=randomSampling(model, nSamples, replaceBoundsWithInf,supressErrors)
0002 % randomSampling
0003 %   Returns a number of random solutions
0004 %
0005 %   model                   a model structure
0006 %   nSamples                the number of solutions to return
0007 %                           (opt, default 1000)
0008 %   replaceBoundsWithInf    replace the largest upper bounds with Inf and
0009 %                           the smallest lower bounds with -Inf. This is
0010 %                           needed in order to get solutions without loops
0011 %                           if your model has for example 1000/-1000 as
0012 %                           arbitary large bounds. If your model only has
0013 %                           "biologically relevant" bounds, then set this
0014 %                           to false (opt, default true)
0015 %   supressErrors           the program will halt if it has problems
0016 %                           finding non-zero solutions which are not
0017 %                           involved in loops. This could be because the
0018 %                           constraints on the model are too relaxed (such
0019 %                           as unlimited glucose uptake) or too strict
0020 %                           (such as too many and too narrow constraints)
0021 %                           (opt, default false)
0022 %
0023 %   solutions               matrix with the solutions
0024 %
0025 %   The solutions are generated by maximizing (with random weights) for a
0026 %   random set of three reactions. For reversible reactions it randomly
0027 %   chooses between maximizing and minimizing.
0028 %
0029 %   Usage: solutions=randomSampling(model, nSamples, replaceBoundsWithInf)
0030 %
0031 %   Rasmus Agren, 2013-04-16
0032 %
0033 
0034 if nargin<2
0035     nSamples=1000;
0036 end
0037 if nargin<3
0038     replaceBoundsWithInf=true;
0039 end
0040 if nargin<4
0041     supressErrors=false;
0042 end
0043 
0044 nRxns=2; %Number of reactions in the objective function in each iteration
0045 
0046 %Simplify the model to speed stuff up a little. Keep original mapping
0047 originalRxns=model.rxns;
0048 model=simplifyModel(model,false,false,true,true);
0049 
0050 %First check that the model is feasible given the constraints
0051 sol=solveLP(model);
0052 if isempty(sol.x)
0053    throw(MException('','The model has no feasible solution')); 
0054 end 
0055 
0056 %Then change the bounds to +/- Inf. This is needed in order to not have
0057 %loops in the solutions
0058 if replaceBoundsWithInf==true
0059    model.ub(model.ub==max(model.ub))=Inf;
0060    model.lb(model.lb==min(model.lb))=-Inf;
0061 end
0062 
0063 %Reactions which can be involved in loops should not be optimized
0064 %for. Check which reactions reach an arbitary high upper bound
0065 goodRxns=true(numel(model.rxns),1);
0066 for i=1:numel(model.rxns)
0067     if goodRxns(i)==true
0068         testModel=setParam(model,'eq',model.rxns(i),1000);
0069         sol=solveLP(testModel);
0070         if ~isempty(sol.f)
0071             goodRxns(abs(sol.x)>999)=false;
0072         else
0073             %If the reaction is reversible, also check in that direction
0074             if model.rev(i)
0075                 testModel=setParam(model,'eq',model.rxns(i),-1000);
0076                 sol=solveLP(testModel);
0077                 if ~isempty(sol.f)
0078                     goodRxns(abs(sol.x)>999)=false;
0079                 end
0080             end
0081         end
0082     end
0083 end
0084 
0085 %Reserve space for a solution matrix
0086 sols=zeros(numel(model.rxns),nSamples);
0087 
0088 %Main loop
0089 counter=1;
0090 badSolutions=0;
0091 goodRxns=find(goodRxns);
0092 while counter<=nSamples
0093    rxns=randsample(numel(goodRxns),nRxns);
0094    model.c=zeros(numel(model.rxns),1);
0095    multipliers=randsample([-1 1],nRxns,true);
0096    multipliers(model.rev(goodRxns(rxns))==0)=1;
0097    model.c(goodRxns(rxns))=rand(nRxns,1).*multipliers;
0098    sol=solveLP(model);
0099    if any(sol.x)
0100        if abs(sol.f)>10^-8
0101             sols(:,counter)=sol.x;
0102             counter=counter+1;
0103             badSolutions=0;
0104        else
0105             badSolutions=badSolutions+1;
0106             %If it only finds bad solutions then throw an error.
0107             if badSolutions==50 && supressErrors==false
0108                 throw(MException('','The program is having problems finding non-zero solutions that are not involved in loops. Review the constraints on your model. Set supressErrors to true to ignore this error.')); 
0109             end
0110        end
0111    end
0112 end
0113 
0114 %Map to original model
0115 [crap I]=ismember(model.rxns,originalRxns);
0116 solutions=zeros(numel(originalRxns),nSamples);
0117 solutions(I,:)=sols;
0118 solutions=sparse(solutions);
0119 end
0120 
0121 %To use instead of the normal Matlab randsample function. This is in order
0122 %to not depend on the Matlab statistical toolbox.
0123 function I=randsample(n,k,replacement)
0124     if nargin<3
0125         replacement=false;
0126     end
0127     %n can be a integer, which leads to I being sampled from 1:n, or it can
0128     %be a population to sample from.
0129     if numel(n)==1 && isnumeric(n)
0130        n=1:n;
0131     end
0132     %Loop and get random numbers until the list is unique. This is only a
0133     %good option is the number of samples is small compared to the
0134     %population. There are several checks that should be made here, for
0135     %example regarding size and that the number of samples is <=population
0136     %size if replacement==false. This is not the case in randomSampling, so
0137     %such checks are ignored
0138     while true
0139         J=randi(numel(n),[k,1]);
0140         if replacement==true || numel(J)==numel(unique(J))
0141             I=n(J);
0142             break;
0143         end
0144     end
0145     I=I(:);
0146 end

Generated on Tue 16-Jul-2013 21:50:02 by m2html © 2005