createHRWarmup

PURPOSE ^

createHRWarmup Create a warmup point set for hit-and-run sampling by

SYNOPSIS ^

function warmupPts= createHRWarmup(model,nPoints,verbFlag,bias,nPointsCheck)

DESCRIPTION ^

 createHRWarmup Create a warmup point set for hit-and-run sampling by
 combining orthogonal and random points

 warmupPts= createHRWarmup(model,nPoints,verbFlag)

INPUTS
 model     Model structure

OPTIONAL INPUTS
 nPoints   Number of warmup points (Default = 5000);
 verbFlag  Verbose flag (Default = false)
 bias
   method          Biasing distribution: 'uniform', 'normal'
   index           The reaction indexes which to bias (nBias total)
   param           nBias x 2 matrix of parameters (for uniform it's min
   max, for normal it's mu, sigma).

OUTPUT
 warmupPts Set of warmup points

 Markus Herrgard 4/21/06

 Richard Que (11/23/09) Integrated subfunctions into script.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function warmupPts= createHRWarmup(model,nPoints,verbFlag,bias,nPointsCheck)
0002 % createHRWarmup Create a warmup point set for hit-and-run sampling by
0003 % combining orthogonal and random points
0004 %
0005 % warmupPts= createHRWarmup(model,nPoints,verbFlag)
0006 %
0007 %INPUTS
0008 % model     Model structure
0009 %
0010 %OPTIONAL INPUTS
0011 % nPoints   Number of warmup points (Default = 5000);
0012 % verbFlag  Verbose flag (Default = false)
0013 % bias
0014 %   method          Biasing distribution: 'uniform', 'normal'
0015 %   index           The reaction indexes which to bias (nBias total)
0016 %   param           nBias x 2 matrix of parameters (for uniform it's min
0017 %   max, for normal it's mu, sigma).
0018 %
0019 %OUTPUT
0020 % warmupPts Set of warmup points
0021 %
0022 % Markus Herrgard 4/21/06
0023 %
0024 % Richard Que (11/23/09) Integrated subfunctions into script.
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 % Generate the correct parameters for the biasing reactions
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         % Find upper & lower bounds for bias rxns to ensure that no
0055         % problems arise with values out of bounds
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 %Generate the points
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                 % Try until get points inside the space
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     % Create random objective function
0118     model.c = rand(nRxns,1)-0.5;
0119     
0120     for maxMin = [1, -1]
0121         % Set the objective function
0122         if i <= nRxns
0123             model.c = zeros(nRxns,1);
0124             model.c(i) = 1;
0125         end
0126         model.osense = maxMin;
0127         
0128         % Determine the max or min for the rxn
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         % Continue if optimal solution is found
0142         
0143         % Move points to within bounds
0144         x(x > model.ub) = model.ub(x > model.ub);
0145         x(x < model.lb) = model.lb(x < model.lb);
0146         
0147         % Store point
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 % Move points in
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 % % Create orthogonal warmup points
0176 % warmupPts1= createHRWarmupOrth(model,true,verbFlag);
0177 %
0178 % % Create warmup points using random directions
0179 % warmupPts2 = createHRWarmupRand(model,1000,verbFlag);
0180 %
0181 % [nRxns,nPts1] = size(warmupPts1);
0182 % [nRxns,nPts2] = size(warmupPts2);
0183 % warmupPts = zeros(nRxns,nPoints);
0184 %
0185 % pointOrder = randperm(nPts1);
0186 %
0187 % if (nPoints < nPts1)
0188 %     warning(['Need a minimum of ' num2str(nPts1) ' warmup points']);
0189 %     nPoints = nPts1;
0190 % end
0191 %
0192 % % Combine point sets
0193 % for i = 1:nPoints
0194 %   if (i <= nPts1)
0195 %       % Ensure that each direction is used at least once
0196 %       x1 = warmupPts1(:,pointOrder(i));
0197 %   else
0198 %       % All direction already used
0199 %       x1 = warmupPts1(:,ceil(rand*nPts1));
0200 %   end
0201 %   x2 = warmupPts2(:,ceil(rand*nPts2));
0202 %   r = rand;
0203 %   thisPoint = sparse((x1*r + x2*(1-r)));
0204 %   % Make sure points stay within bounds
0205 %   thisPoint(thisPoint > model.ub) = model.ub(thisPoint ...
0206 %                                              > ...
0207 %                                              model.ub);
0208 %   thisPoint(thisPoint < model.lb) = model.lb(thisPoint ...
0209 %                                              < model.lb);
0210 %   warmupPts(:,i) = thisPoint;
0211 %
0212 % end
0213

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003