gpSampler

PURPOSE ^

gpSampler Samples an arbitrary linearly constrained space using a fixed

SYNOPSIS ^

function [sampleStructOut, mixedFrac] = gpSampler(sampleStruct, nPoints, bias, maxTime, maxSteps, threads, nPointsCheck)

DESCRIPTION ^

gpSampler Samples an arbitrary linearly constrained space using a fixed
number of points that are moved in parallel

 [sampleStructOut, mixedFraction] = gpSampler(sampleStruct, nPoints, bias, maxTime, maxSteps)

 The space is defined by
   A x <=,==,>= b
   lb <= x <= ub

INPUTS
 sampleStruct      Structure describing the space to be sampled and
                   previous point sets
   A               LHS matrix (optionally, if not A script checks for S)
   b               RHS vector
   lb              Lower bound
   ub              Upper bound
   csense          Constraint type for each row in A ('G', 'L', 'E')
   warmupPoints     Set of warmup points (optional, generated by default)
     points          Currently sampled points (optional)

OPTIONAL INPUTS
 nPoints           Number of points used in sampling 
                   (default = 2*nRxns or 5000 whichever is greater)
 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).
 maxTime           Maximum time alloted for the sampling in seconds
                   (default 600 s, pass an empty number [] to set maxSteps instead)
 maxSteps          Maximum number of steps to take (default 1e10). Sampler
                   will run until either maxStep or maxTime is reached.  
                   Set maxStep or maxTime to 0 and no sampling will occur 
                   (only warmup points generated).  
 threads           number of threads the sampler will use.  If you have a
                   dual core machine, you can set it to 2 etc.  The speed
                   up is almost linear w/ the number of cores.
                   If using this feature and 2009a or newer, a futher 
                   speedup can be obtained by starting matlab from the 
                   command line by "typing matlab -singleCompThread"
                   New feature:  if threads < 0, use distributed toolbox.
 nPointsCheck      Checks that minimum number of points (2*nRxns) are
                   used. (Default = true).

OUTPUT
 sampleStructOut   The sampling structure with some extra fields.
 mixedFract        The fraction mixed (relative to the warmupPts).  A value of 1
                   means not mixed at all.  A value of .5 means completely mixed.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [sampleStructOut, mixedFrac] = gpSampler(sampleStruct, nPoints, bias, maxTime, maxSteps, threads, nPointsCheck)
0002 %gpSampler Samples an arbitrary linearly constrained space using a fixed
0003 %number of points that are moved in parallel
0004 %
0005 % [sampleStructOut, mixedFraction] = gpSampler(sampleStruct, nPoints, bias, maxTime, maxSteps)
0006 %
0007 % The space is defined by
0008 %   A x <=,==,>= b
0009 %   lb <= x <= ub
0010 %
0011 %INPUTS
0012 % sampleStruct      Structure describing the space to be sampled and
0013 %                   previous point sets
0014 %   A               LHS matrix (optionally, if not A script checks for S)
0015 %   b               RHS vector
0016 %   lb              Lower bound
0017 %   ub              Upper bound
0018 %   csense          Constraint type for each row in A ('G', 'L', 'E')
0019 %   warmupPoints     Set of warmup points (optional, generated by default)
0020 %     points          Currently sampled points (optional)
0021 %
0022 %OPTIONAL INPUTS
0023 % nPoints           Number of points used in sampling
0024 %                   (default = 2*nRxns or 5000 whichever is greater)
0025 % bias
0026 %   method          Biasing distribution: 'uniform', 'normal'
0027 %   index           The reaction indexes which to bias (nBias total)
0028 %   param           nBias x 2 matrix of parameters (for uniform it's min
0029 %                   max, for normal it's mu, sigma).
0030 % maxTime           Maximum time alloted for the sampling in seconds
0031 %                   (default 600 s, pass an empty number [] to set maxSteps instead)
0032 % maxSteps          Maximum number of steps to take (default 1e10). Sampler
0033 %                   will run until either maxStep or maxTime is reached.
0034 %                   Set maxStep or maxTime to 0 and no sampling will occur
0035 %                   (only warmup points generated).
0036 % threads           number of threads the sampler will use.  If you have a
0037 %                   dual core machine, you can set it to 2 etc.  The speed
0038 %                   up is almost linear w/ the number of cores.
0039 %                   If using this feature and 2009a or newer, a futher
0040 %                   speedup can be obtained by starting matlab from the
0041 %                   command line by "typing matlab -singleCompThread"
0042 %                   New feature:  if threads < 0, use distributed toolbox.
0043 % nPointsCheck      Checks that minimum number of points (2*nRxns) are
0044 %                   used. (Default = true).
0045 %
0046 %OUTPUT
0047 % sampleStructOut   The sampling structure with some extra fields.
0048 % mixedFract        The fraction mixed (relative to the warmupPts).  A value of 1
0049 %                   means not mixed at all.  A value of .5 means completely mixed.
0050 
0051 %% Parameter Processing / error checking
0052 sampleStructOut = 0; % in case of returning early
0053 mixedFrac = 1; % in case of returning early
0054 if nargin < 2
0055     nPoints = 5000;
0056 end
0057 if nargin < 3
0058     bias = [];
0059 end
0060 if ~isempty(bias)
0061     if ~isfield (bias, 'method')
0062         display('bias does not have a method set');
0063         return;
0064     end
0065 end
0066 if nargin < 4 || (isempty(maxTime) && isempty(maxSteps))
0067     maxTime = 10*60; % 10 minutes
0068 end
0069 if (nargin < 5) || isempty(maxSteps)
0070     % Max time takes precedence
0071     maxSteps = 1e10;
0072 else
0073     % Set max steps instead of max time
0074     if (isempty(maxTime))
0075         maxTime = 1e10;
0076     end
0077 end
0078 
0079 if nargin < 6 || isempty(threads)
0080     threads = 1;
0081 end
0082 
0083 if nargin < 7, nPointsCheck = true; end
0084 
0085 % Sanity checking
0086 if (~ isfield (sampleStruct, 'A'))
0087     if isfield(sampleStruct, 'S')
0088         display('A set to S');
0089         sampleStruct.A = sampleStruct.S;
0090     else
0091         display('A and/or S not set');
0092         return;
0093     end
0094 end
0095 if (~ isfield (sampleStruct, 'b'))
0096     sampleStruct.b = zeros(size(sampleStruct.A,1), 1);
0097     display('Warning:  b not set.  Defaulting to zeros');
0098 end
0099 if (~ isfield (sampleStruct, 'csense'))
0100     sampleStruct.csense(1:size(sampleStruct.A,1)) = 'E';
0101     display('Warning:  csense not set.  Defaulting to all Equality constraints');
0102 end
0103 if (~isfield (sampleStruct, 'lb'))
0104     display('lb not set');
0105     return;
0106 end
0107 if (~isfield (sampleStruct, 'ub'))
0108     display('ub not set');
0109     return;
0110 end
0111 
0112 %% internal data generation
0113 % make internal structure
0114 [A, b, csense, lb, ub] = deal(sampleStruct.A, sampleStruct.b, sampleStruct.csense, sampleStruct.lb, sampleStruct.ub);
0115 
0116 % constInd = (lb == ub);
0117 % constVal = lb(constInd);
0118 % Aconst = A(:,constInd);
0119 % b = b - Aconst*constVal;
0120 % A = A(:,~constInd);
0121 % lb = lb(~constInd);
0122 % ub = ub(~constInd);
0123 % [sampleStruct.A, sampleStruct.b, sampleStruct.csense, sampleStruct.lb, sampleStruct.ub] = deal(A, b, csense, lb, ub);
0124 
0125 [rA, dimx] = size(A);
0126 if (~ isfield(sampleStruct, 'internal'))
0127     Anew = sparse(0, dimx);
0128     Cnew = sparse(0, dimx);
0129     Bnew = zeros(rA, 1);
0130     Dnew = zeros(rA, 1);
0131     rAnew = 0;
0132     rCnew = 0 ;
0133     for i = 1:size(A, 1)
0134         switch csense(i)
0135             case 'E'
0136                 rAnew = rAnew+1;
0137                 Anew(rAnew,:) = A(i,:);
0138                 Bnew(rAnew,:) = b(i);
0139 
0140             case 'G'
0141                 rCnew=rCnew+1;
0142                 Cnew(rCnew,:) = -A(i,:);
0143                 Dnew(rCnew,:) = -b(i);
0144             case 'L'
0145                 rCnew=rCnew+1;
0146                 Cnew(rCnew,:) = A(i,:);
0147                 Dnew(rCnew,:) = b(i);
0148             otherwise
0149                 display ('whoops.  csense can only contain E, G, or L')
0150                 return;
0151         end
0152     end
0153     
0154     %Anew = Anew(1:rAnew,:);
0155     Bnew = Bnew(1:rAnew,:);
0156     Dnew = Dnew(1:rCnew,:);
0157 
0158     % calculate offset
0159     if find(Bnew ~= 0)
0160         offset = Anew\Bnew;
0161     else
0162         offset = zeros(size(Anew,2), 1);
0163     end
0164 
0165     % rescale Bnew, Dnew
0166     Boffset = Bnew - Anew*offset;
0167     if (max(abs(Boffset)) > .0000000001)
0168         display('whoops.  It looks like the offset calculation made a mistake.  this should be zero.');
0169         max(abs(Boffset))
0170         return;
0171     end
0172     Doffset = Dnew - Cnew*offset;
0173 
0174     lbnew = lb - offset;
0175     ubnew = ub - offset;
0176     
0177     sampleStruct.internal.offset = offset;
0178     sampleStruct.internal.Anew = Anew;
0179     sampleStruct.internal.Cnew = Cnew;
0180     sampleStruct.internal.Dnew = Doffset;
0181     sampleStruct.internal.lbnew = lbnew;
0182     sampleStruct.internal.ubnew = ubnew;
0183     
0184     if (isfield(sampleStruct, 'warmupPts'))
0185         sampleStruct = rmfield(sampleStruct, 'warmupPts');
0186     end
0187     if ~isempty(bias)
0188         sampleStruct.internal.fixed = bias.index;
0189     else
0190         sampleStruct.internal.fixed = [];
0191     end
0192 end
0193 
0194 %% Generate warmup points
0195 if (~ isfield(sampleStruct, 'warmupPts') )
0196     fprintf('Generating warmup points\n');
0197 %     warmupPts = warmup(sampleStruct, nPoints, bias);
0198     warmupPts = createHRWarmup(sampleStruct, nPoints, false, bias, nPointsCheck);
0199     sampleStruct.warmupPts = warmupPts;
0200     if (isfield(sampleStruct, 'points'))
0201         sampleStruct = rmfield(sampleStruct, 'points');
0202     end
0203     if (isfield(sampleStruct, 'bias'))
0204         sampleStruct = rmfield(sampleStruct, 'bias');
0205     end
0206     sampleStruct.steps = 0;
0207     save sampleStructTmp sampleStruct
0208 else
0209     fprintf('Warmup points already present.\n');
0210 %    save sampleStructTmp sampleStruct
0211 end
0212 
0213 %% Do actual sampling
0214 fprintf('Sampling\n');
0215 if(maxTime > 0 && maxSteps > 0)
0216     if threads < 0  %uses distributed toolbox.
0217         sampleStruct = ACHRSamplerDistributedGeneral(sampleStruct, ceil(maxSteps/50), 50, maxTime);
0218     else
0219         sampleStruct = ACHRSamplerParallelGeneral(sampleStruct, ceil(maxSteps/50), 50, maxTime, threads);
0220     end
0221     mixedFrac = mixFraction(sampleStruct.points, sampleStruct.warmupPts, sampleStruct.internal.fixed);
0222 else
0223     mixedFrac = 1;
0224 end
0225 
0226 sampleStructOut = sampleStruct;
0227 
0228 return;
0229 
0230 
0231 
0232 %% warmup Point generator
0233 function warmupPts = warmup(sampleProblem, nPoints, bias)
0234 
0235 dimX = size(sampleProblem.A, 2); 
0236 warmupPts = zeros(dimX, nPoints);
0237 
0238 [LPproblem.A,LPproblem.b,LPproblem.lb,LPproblem.ub,LPproblem.csense,LPproblem.osense] = ...
0239     deal(sampleProblem.A,sampleProblem.b,sampleProblem.lb,sampleProblem.ub,sampleProblem.csense,1);
0240 
0241 % Generate the correct parameters for the biasing reactions
0242 if ~isempty(bias)
0243     if (~ismember(bias.method,{'uniform','normal'}))
0244         error('Biasing method not implemented');
0245     end
0246     for k = 1:size(bias.index)
0247         ind = bias.index(k);
0248         % Find upper & lower bounds for bias rxns to ensure that no
0249         % problems arise with values out of bounds
0250         LPproblem.c = zeros(size(LPproblem.A,2),1);
0251         LPproblem.c(ind) = 1;
0252         LPproblem.osense = -1;
0253         sol = solveCobraLP(LPproblem);
0254         maxFlux = sol.obj;
0255         LPproblem.osense = 1;
0256         sol = solveCobraLP(LPproblem);
0257         minFlux = sol.obj;
0258 
0259         if strcmp(bias.method, 'uniform')
0260             upperBias = bias.param(k,2);
0261             lowerBias = bias.param(k,1);
0262             if (upperBias > maxFlux || upperBias < minFlux)
0263                 upperBias = maxFlux;
0264                 disp('Invalid bias bounds - using default bounds instead');
0265             end             
0266             if (lowerBias < minFlux || lowerBias > maxFlux)
0267                 lowerBias = minFlux;
0268                 disp('Invalid bias bounds - using default bounds instead');
0269             end
0270             bias.param(k,1) = lowerBias;
0271             bias.param(k,2) = upperBias;
0272         elseif strcmp(bias.method, 'normal')
0273             biasMean = bias.param(k,1);
0274             if (biasMean > maxFlux || biasMean < minFlux)
0275                  bias.param(k,1) = (minFlux + maxFlux)/2;
0276                 disp('Invalid bias mean - using default mean instead');
0277             end
0278             biasFluxMin(k) = minFlux;
0279             biasFluxMax(k) = maxFlux;
0280         end
0281     end
0282 end
0283 
0284 %Generate the points
0285 i = 1;
0286 while i <= nPoints/2
0287     if mod(i,10) ==0
0288         fprintf('%d\n',2*i);
0289     end
0290     if ~isempty(bias)
0291         for k = 1:size(bias.index)
0292             ind = bias.index(k);
0293             if strcmp(bias.method, 'uniform')
0294                 diff = bias.param(k,2) - bias.param(k,1);
0295                 fluxVal = diff*rand() + bias.param(k,1);
0296             elseif strcmp(bias.method, 'normal')
0297                 valOK = false;
0298                 % Try until get points inside the space
0299                 while (~valOK)
0300                     fluxVal = randn()*bias.param(k,2)+bias.param(k,1);
0301                     if (fluxVal <= biasFluxMax(k) && fluxVal >= biasFluxMin(k))
0302                         valOK = true;
0303                     end
0304                 end
0305             end
0306             LPproblem.lb(ind) = 0.99999999*fluxVal;
0307             LPproblem.ub(ind) = 1.00000001*fluxVal;
0308         end
0309     end
0310     
0311     % Pick the next flux to optimize, cycles though each reaction
0312     % alternates minimization and maximization for each cycle
0313 
0314     
0315     
0316     LPproblem.c = rand(dimX, 1)-.5;
0317     validFlag = true;
0318     
0319     for maxMin = [1, -1]  
0320       % Set the objective function
0321       if i <= dimX
0322 %        LPproblem.c = rand(dimX, 1)-.5;
0323         LPproblem.c(i) = 5000;
0324       end
0325       LPproblem.osense = maxMin;
0326       
0327       % Determine the max or min for the rxn
0328       sol = solveCobraLP(LPproblem);
0329       x = sol.full;
0330       if maxMin == 1
0331           sol1 = sol;
0332       else
0333           sol2 = sol;
0334       end
0335       status = sol.stat;
0336       if status ~= 1
0337           display ('invalid solution')
0338           validFlag = false;
0339           display(status)
0340           pause;
0341       end
0342         
0343       % Move points to within bounds
0344       x(x > LPproblem.ub) = LPproblem.ub(x > LPproblem.ub);
0345       x(x < LPproblem.lb) = LPproblem.lb(x < LPproblem.lb); 
0346       
0347       % Store point
0348       % For non-random points just store a min/max point
0349       
0350       if (maxMin == 1)
0351           warmupPts(:,2*i-1) = x;
0352       else 
0353           warmupPts(:,2*i) = x;
0354       end
0355     end
0356     
0357     if validFlag
0358         %postprocess(LPproblem, sol1, sol2)
0359         i = i+1;
0360     end
0361 end
0362 centerPoint = mean(warmupPts,2);
0363 
0364 % Move points in
0365 if isempty(bias)
0366     warmupPts = warmupPts*.33 + .67*centerPoint*ones(1,nPoints);
0367 else
0368     warmupPts = warmupPts*.99 + .01*centerPoint*ones(1,nPoints);
0369 end
0370 
0371 % Permute point order
0372 % if (permFlag)
0373 %     [nRxns,nPoints] = size(warmupPts);
0374 %     warmupPts = warmupPts(:,randperm(nPoints));
0375 % end
0376 return;
0377 
0378 %% post processing for better warmup point generation.
0379 % function out = postprocess(LPproblem, sol1, sol2)
0380 % x1 = sol1.full;
0381 % x2 = sol2.full;
0382 % closetoboundary(LPproblem, x1)
0383 % closetoboundary(LPproblem, x2)
0384 %
0385 % separation = sol2.obj - sol1.obj;
0386 % if separation < .00001
0387 %     disp('low separation');
0388 %     pause;
0389 % end
0390 % LPproblem.A = [LPproblem.A; LPproblem.c];
0391 % LPproblem.b(end+1) = sol1.obj + separation*.1;
0392 % Lpproblem.csense(end+1) = 'L';
0393 %
0394 %
0395 %
0396 % pause;
0397 % out = 2;
0398 %
0399 % return;
0400 
0401 % function counter = closetoboundary(LPproblem, sol)
0402 % etol = 1e-5;
0403 % counter = 0;
0404 % counter = counter + sum(LPproblem.A((LPproblem.csense == 'G'),:)*sol - LPproblem.b(LPproblem.csense =='G') < etol)
0405 % counter = counter + sum(LPproblem.b(LPproblem.csense =='L') - LPproblem.A((LPproblem.csense == 'L'),:)*sol < etol)
0406 % counter = counter + sum(LPproblem.ub - sol < etol)
0407 % counter = counter + sum(sol - LPproblem.lb < etol)
0408 % return;
0409 
0410

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