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
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