ACBSampler

PURPOSE ^

ACBSampler Artificial centering boundary sampler

SYNOPSIS ^

function ACBSampler(model,warmupPoints,fileName,nFiles,pointsPerFile,nMixPts,nWarmupNeeded,saveMatFlag,biasOpt)

DESCRIPTION ^

 ACBSampler Artificial centering boundary sampler

 ACBSampler(model,warmupPoints,fileName,nFiles,pointsPerFile,nMixPts,nWarmupNeeded,saveMatFlag,biasOpt)

INPUTS
 model         Model structure
 warmupPoints  Warmup points
 fileName      Base fileName for saving results
 nFiles        Number of sample point files created
 pointsPerFile Number of points per file saved
 nMixPts       Number of steps initially used for mixing (not saved)

OPTIONAL INPUTS
 nWarmupNeeded Number of warmup points needed (Default = 20000)
 saveMatFlag   Save points in mat format vs txt format (Default = true)
 biasOpt       Options for biasing sampler (Default = no bias)
   rxns
   values
   stds
   percThr

 Christian Barrett and Markus Herrgard 8/24/06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function ACBSampler(model,warmupPoints,fileName,nFiles,pointsPerFile,nMixPts,nWarmupNeeded,saveMatFlag,biasOpt)
0002 % ACBSampler Artificial centering boundary sampler
0003 %
0004 % ACBSampler(model,warmupPoints,fileName,nFiles,pointsPerFile,nMixPts,nWarmupNeeded,saveMatFlag,biasOpt)
0005 %
0006 %INPUTS
0007 % model         Model structure
0008 % warmupPoints  Warmup points
0009 % fileName      Base fileName for saving results
0010 % nFiles        Number of sample point files created
0011 % pointsPerFile Number of points per file saved
0012 % nMixPts       Number of steps initially used for mixing (not saved)
0013 %
0014 %OPTIONAL INPUTS
0015 % nWarmupNeeded Number of warmup points needed (Default = 20000)
0016 % saveMatFlag   Save points in mat format vs txt format (Default = true)
0017 % biasOpt       Options for biasing sampler (Default = no bias)
0018 %   rxns
0019 %   values
0020 %   stds
0021 %   percThr
0022 %
0023 % Christian Barrett and Markus Herrgard 8/24/06
0024 
0025 if (nargin < 7)
0026     % Expands the number of warmup points if the initially provided set is too
0027     % small
0028     nWarmupNeeded = 20000;
0029 end
0030 if (nargin < 8)
0031     saveMatFlag = true;
0032 end
0033 if (nargin < 9)
0034     biasFlag = false;
0035 else
0036     biasFlag = true;
0037 end
0038 
0039 warning off MATLAB:divideByZero;
0040 
0041 % Minimum allowed distance to the closest constraint
0042 maxMinTol = 1e-9;
0043 % Ignore directions where u is really small
0044 uTol = 1e-9;
0045 % Project out of directions that are too close to the boundary
0046 dTol = 1e-10;
0047 
0048 % Number of initial warmup points
0049 [nRxns,nWarmup] = size(warmupPoints);
0050 
0051 % Find indices for fluxes to be biased
0052 if (biasFlag)
0053     [tmp,biasInd] = ismember(biasOpt.rxns,model.rxns);
0054     if (length(biasInd) > 1)
0055         biasDist = sum((warmupPoints(biasInd,:)-repmat(biasOpt.values,1,nWarmup)).^2./repmat(biasOpt.stds.^2,1,nWarmup));
0056     else
0057         biasDist = (warmupPoints(biasInd,:)-repmat(biasOpt.values,1,nWarmup)).^2./repmat(biasOpt.stds.^2,1,nWarmup);
0058     end
0059 end
0060 
0061 totalPointCount = 0;
0062 
0063 fidErr = fopen([fileName '.err'],'w');
0064 
0065 for fileNo = 1:nFiles
0066 
0067     if (~saveMatFlag)
0068         fid = fopen([fileName '_' num2str(fileNo) '.fls'],'w');
0069     else
0070         % Allocate memory for all points
0071         points = zeros(nRxns,pointsPerFile);
0072     end
0073 
0074     pointCount = 1;
0075     while (pointCount <= pointsPerFile)
0076 
0077         if (mod(totalPointCount,500) == 0)
0078             fprintf('%d %d %d\n',fileNo,pointCount,totalPointCount);
0079         end
0080 
0081         % Pick two random points
0082         if (biasFlag)
0083            [tmp,sortInd] = sort(biasDist);
0084             randPoint1 = warmupPoints(:,sortInd(ceil(biasOpt.percThr*nWarmup*rand)));
0085             randPoint2 = warmupPoints(:,sortInd(ceil(biasOpt.percThr*nWarmup*rand)));
0086         else
0087             randPoint1 = warmupPoints(:,ceil(nWarmup*rand));
0088             randPoint2 = warmupPoints(:,ceil(nWarmup*rand));
0089         end
0090 
0091         % Get a line formed by the two
0092         u = (randPoint2-randPoint1);
0093         u = u/norm(u);
0094 
0095         % This gets randPoint1 away from the boundary, which causes problems
0096         % when trying to move from randPoint1
0097         % randPoint1 = randPoint1 + tol*u;
0098 
0099         % Figure out the distances to upper and lower bounds
0100         distUb = (model.ub - randPoint1);
0101         distLb = (randPoint1 - model.lb);
0102 
0103         % Figure out if we are too close to a boundary
0104         validDir = ((distUb > dTol) & (distLb > dTol));
0105 
0106         % Figure out positive and negative directions
0107         posDirn = find(u(validDir) > uTol);
0108         negDirn = find(u(validDir) < -uTol);
0109 
0110         if (isempty(posDirn) & isempty(negDirn))
0111             continue
0112         end
0113 
0114         % Figure out all the possible maximum and minimum step sizes
0115         maxStepTemp = distUb(validDir)./u(validDir);
0116         minStepTemp = -distLb(validDir)./u(validDir);
0117         maxStepVec = [maxStepTemp(posDirn);minStepTemp(negDirn)];
0118         minStepVec = [minStepTemp(posDirn);maxStepTemp(negDirn)];
0119 
0120         % Figure out the true max & min step sizes
0121         maxStep = min(maxStepVec);
0122         minStep = max(minStepVec);
0123 
0124         % Compute the boundary points along this line
0125         boundaryPoint1 = randPoint1 + minStep*u;
0126         boundaryPoint2 = randPoint1 + maxStep*u;
0127 
0128         % Move on if the points left the null space for some reason
0129         error1 = full(max(max(abs(model.S*boundaryPoint1))));
0130         error2 = full(max(max(abs(model.S*boundaryPoint2))));
0131         if (error1 > 1e-7 | error2 > 1e-7)
0132           fprintf('Point out of N-space: %g %g\n',error1,error2);
0133           continue;
0134         end
0135 
0136         % Get the center point
0137         centerPoint = (boundaryPoint1 + boundaryPoint2)/2.0;
0138 
0139         % Check if we want to add this point to warmup points or replace a
0140         % previous point
0141 
0142         if (biasFlag)
0143             % Calculate distance to bias flux values
0144             [maxBiasDist,maxInd] = max(biasDist);
0145 
0146             if (length(biasInd) > 1)
0147                 biasDistThisStep = sum((centerPoint(biasInd)-biasOpt.values).^2./biasOpt.stds.^2);
0148             else
0149                 biasDistThisStep = (centerPoint(biasInd)-biasOpt.values).^2./biasOpt.stds.^2;
0150             end
0151 
0152             if (mod(totalPointCount,100) == 0)
0153               minBiasDist = min(biasDist);
0154               fprintf('%d\t%f\t%f\t%f\t%f\t%d\n',totalPointCount, ...
0155                       biasDistThisStep,minBiasDist,maxBiasDist,mean(biasDist),nWarmup);
0156               fprintf(fidErr,'%d\t%f\t%f\t%f\t%f\t%d\n',totalPointCount, ...
0157                       biasDistThisStep,minBiasDist,maxBiasDist,mean(biasDist),nWarmup);
0158               
0159             end
0160 
0161             if (biasDistThisStep < maxBiasDist)
0162                 replaceWarmup = true;
0163             else
0164                 replaceWarmup = false;
0165             end
0166         else
0167             replaceWarmup = true;
0168         end
0169         
0170         nWarmup = size(warmupPoints,2);
0171         if (replaceWarmup)
0172             if (biasFlag)
0173                 % Add warmup point
0174                 replaceInd = max(find(biasDistThisStep > ...
0175                                       biasDist));
0176                 if (nWarmup < nWarmupNeeded)
0177                     biasDist(end+1) = biasDistThisStep;
0178                     warmupPoints(:,end+1) = centerPoint;
0179                     nWarmup = nWarmup+1;
0180                 else
0181                     % Replace warmup point
0182                     biasDist(maxInd) = biasDistThisStep;
0183                     warmupPoints(:,maxInd) = centerPoint;
0184                 end
0185             else
0186                 if (nWarmup < nWarmupNeeded)
0187                     warmupPoints(:,end+1) = centerPoint;
0188                     nWarmup = nWarmup+1;
0189                 end
0190             end
0191         end
0192         
0193         % Count the total number of points generated
0194         totalPointCount = totalPointCount + 1;
0195 
0196         % Check if it is time to start saving points
0197         if (totalPointCount > nMixPts)
0198             if (rand < 0.5)
0199                 pointToSave = boundaryPoint2;
0200             else
0201                 pointToSave = boundaryPoint1;
0202             end
0203 
0204             if (saveMatFlag)
0205                 points(:,pointCount) = pointToSave;
0206             else
0207                 % Print out point
0208                 for i = 1:length(pointToSave)
0209                     fprintf(fid,'%g ',pointToSave(i));
0210                 end
0211                 fprintf(fid,'\n');
0212             end
0213 
0214             pointCount = pointCount + 1;
0215         end
0216     end
0217 
0218     % Save current points within the cycle to a file
0219     if (saveMatFlag)
0220         file = [fileName '_' num2str(fileNo) '.mat'];
0221         save (file,'points');
0222     else
0223         fclose(fid);
0224     end
0225 
0226 end % Files
0227 
0228 fclose(fidErr);
0229 
0230 
0231

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