0001 function warmupPts= createHRWarmup(model,nPoints,verbFlag,bias,nPointsCheck)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026 if (nargin < 2)||isempty(nPoints), nPoints = 5000; end
0027 if (nargin < 3)||isempty(verbFlag), verbFlag = false; end
0028 if (nargin < 4), bias = []; end
0029 if (nargin < 5)||isempty(nPointsCheck), nPointsCheck = true; end
0030
0031 if isfield(model,'A')
0032 [nMets,nRxns] = size(model.A);
0033 else
0034 [nMets,nRxns] = size(model.S);
0035 model.A=model.S;
0036 end
0037 if ~isfield(model,'csense')
0038 model.csense(1:size(model.S,1)) = 'E';
0039 end
0040
0041 if nPointsCheck && (nPoints < nRxns*2)
0042 warning(['Need a minimum of ' num2str(nRxns*2) ' warmup points']);
0043 nPoints = nRxns*2;
0044 end
0045 warmupPts = sparse(nRxns,nPoints);
0046
0047
0048 if ~isempty(bias)
0049 if (~ismember(bias.method,{'uniform','normal'}))
0050 error('Biasing method not implemented');
0051 end
0052 for k = 1:size(bias.index)
0053 ind = bias.index(k);
0054
0055
0056 model.c = zeros(nRxns,1);
0057 model.c(ind) = 1;
0058 model.osense = -1;
0059 sol = solveCobraLP(model);
0060 maxFlux = sol.obj;
0061 model.osense = 1;
0062 sol = solveCobraLP(model);
0063 minFlux = sol.obj;
0064
0065 if strcmp(bias.method, 'uniform')
0066 upperBias = bias.param(k,2);
0067 lowerBias = bias.param(k,1);
0068 if (upperBias > maxFlux || upperBias < minFlux)
0069 upperBias = maxFlux;
0070 disp('Invalid bias bounds - using default bounds instead');
0071 end
0072 if (lowerBias < minFlux || lowerBias > maxFlux)
0073 lowerBias = minFlux;
0074 disp('Invalid bias bounds - using default bounds instead');
0075 end
0076 bias.param(k,1) = lowerBias;
0077 bias.param(k,2) = upperBias;
0078 elseif strcmp(bias.method, 'normal')
0079 biasMean = bias.param(k,1);
0080 if (biasMean > maxFlux || biasMean < minFlux)
0081 bias.param(k,1) = (minFlux + maxFlux)/2;
0082 disp('Invalid bias mean - using default mean instead');
0083 end
0084 biasFluxMin(k) = minFlux;
0085 biasFluxMax(k) = maxFlux;
0086 end
0087 end
0088 end
0089
0090 i = 1;
0091 h = waitbar(0,'Creating warmup points ...');
0092
0093 while i <= nPoints/2
0094 if mod(i,10) == 0
0095 waitbar(2*i/nPoints,h);
0096 end
0097 if ~isempty(bias)
0098 for k = 1:size(bias.index)
0099 ind = bias.index(k);
0100 if strcmp(bias.method, 'uniform')
0101 diff = bias.param(k,2) - bias.param(k,1);
0102 fluxVal = diff*rand() + bias.param(k,1);
0103 elseif strcmp(bias.method, 'normal')
0104 valOK = false;
0105
0106 while (~valOK)
0107 fluxVal = randn()*bias.param(k,2)+bias.param(k,1);
0108 if (fluxVal <= biasFluxMax(k) && fluxVal >= biasFluxMin(k))
0109 valOK = true;
0110 end
0111 end
0112 end
0113 model.lb(ind) = 0.99999999*fluxVal;
0114 model.ub(ind) = 1.00000001*fluxVal;
0115 end
0116 end
0117
0118 model.c = rand(nRxns,1)-0.5;
0119
0120 for maxMin = [1, -1]
0121
0122 if i <= nRxns
0123 model.c = zeros(nRxns,1);
0124 model.c(i) = 1;
0125 end
0126 model.osense = maxMin;
0127
0128
0129 sol = solveCobraLP(model);
0130 x = sol.full;
0131 status = sol.stat;
0132 if status == 1
0133 validFlag = true;
0134 else
0135 display ('invalid solution')
0136 validFlag = false;
0137 display(status)
0138 pause;
0139 end
0140
0141
0142
0143
0144 x(x > model.ub) = model.ub(x > model.ub);
0145 x(x < model.lb) = model.lb(x < model.lb);
0146
0147
0148 if (maxMin == 1)
0149 warmupPts(:,2*i-1) = x;
0150 else
0151 warmupPts(:,2*i) = x;
0152 end
0153
0154 if (verbFlag)
0155 if mod(i,100)==0
0156 fprintf('%4.1f\n',i/nPoints*100);
0157 end
0158 end
0159
0160
0161 end
0162 if validFlag
0163 i = i+1;
0164 end
0165 end
0166 centerPoint = mean(warmupPts,2);
0167
0168 if isempty(bias)
0169 warmupPts = warmupPts*.33 + .67*centerPoint*ones(1,nPoints);
0170 else
0171 warmupPts = warmupPts*.99 + .01*centerPoint*ones(1,nPoints);
0172 end
0173 close(h);
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213