ACHRSampler

PURPOSE ^

ACHRSampler Artificial Centering Hit-and-Run sampler

SYNOPSIS ^

function ACHRSampler(model,warmupPoints,fileName,nFiles,pointsPerFile,stepsPerPoint,initPoint,fileBaseNo,maxTime)

DESCRIPTION ^

 ACHRSampler Artificial Centering Hit-and-Run sampler

 ACHRSampler(model,warmupPoints,fileName,nFiles,pointsPerFile,stepsPerPoint,initPoint,fileBaseNo,maxTime)

INPUTS
 model    Model structure
 warmupPoints Warmup points
 fileName Base fileName for saving results
 nFiles   Number of files created
 pointsPerFile Number of points per file saved
 stepsPerPoint Number of sampler steps per point saved

OPTIONAL INPUTS
 initPoint     Initial point (Default = centerPoint)
 fileBaseNo    Base file number for continuing previous sampler run 
               (Default = 0)
 maxTime       Maximum time limit (Default = 36000)

 Markus Herrgard, Gregory Hannum, Ines Thiele, Nathan Price 4/14/06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function ACHRSampler(model,warmupPoints,fileName,nFiles,pointsPerFile,stepsPerPoint,initPoint,fileBaseNo,maxTime)
0002 % ACHRSampler Artificial Centering Hit-and-Run sampler
0003 %
0004 % ACHRSampler(model,warmupPoints,fileName,nFiles,pointsPerFile,stepsPerPoint,initPoint,fileBaseNo,maxTime)
0005 %
0006 %INPUTS
0007 % model    Model structure
0008 % warmupPoints Warmup points
0009 % fileName Base fileName for saving results
0010 % nFiles   Number of files created
0011 % pointsPerFile Number of points per file saved
0012 % stepsPerPoint Number of sampler steps per point saved
0013 %
0014 %OPTIONAL INPUTS
0015 % initPoint     Initial point (Default = centerPoint)
0016 % fileBaseNo    Base file number for continuing previous sampler run
0017 %               (Default = 0)
0018 % maxTime       Maximum time limit (Default = 36000)
0019 %
0020 % Markus Herrgard, Gregory Hannum, Ines Thiele, Nathan Price 4/14/06
0021 
0022 warning off MATLAB:divideByZero;
0023 
0024 if (nargin < 8)
0025   fileBaseNo = 0;
0026 else
0027     if (isempty(fileBaseNo))
0028         fileBaseNo = 0;
0029     end
0030 end
0031 if (nargin < 9)
0032     maxTime = 10*3600;
0033 end
0034 
0035 N = null(full(model.S));
0036 
0037 % Minimum allowed distance to the closest constraint
0038 maxMinTol = 1e-9;
0039 % Ignore directions where u is really small
0040 uTol = 1e-9; 
0041 % Project out of directions that are too close to the boundary
0042 dTol = 1e-14;
0043 
0044 % Number of warmup points
0045 [nRxns,nWrmup] = size(warmupPoints);
0046 
0047 % Find the center of the space
0048 centerPoint = mean(warmupPoints,2);
0049 
0050 % Set the start point
0051 if (nargin < 7)
0052     prevPoint = centerPoint;
0053 else
0054     if (~isempty(initPoint))
0055         prevPoint = initPoint;
0056     else
0057         prevPoint = centerPoint;
0058     end
0059 end
0060 
0061 fidErr = fopen('ACHRerror.txt','w');
0062 
0063 totalStepCount = 0;
0064 
0065 h = waitbar(0,'ACHR sampling in progress ...');
0066 totalCount = nFiles*pointsPerFile*stepsPerPoint;
0067 
0068 t0 = cputime;
0069 fprintf('File #\tPoint #\tStep #\tTime\t#Time left\n');
0070 for i = 1:nFiles
0071 
0072     % Allocate memory for all points
0073     points = zeros(nRxns,pointsPerFile); 
0074     
0075     pointCount = 1;
0076     while (pointCount <= pointsPerFile)
0077             
0078         % Create the random step size vector
0079         randVector = rand(stepsPerPoint,1);
0080         
0081         stepCount = 1;
0082         while (stepCount <= stepsPerPoint)
0083             
0084             % Pick a random warmup point
0085             randPointID = ceil(nWrmup*rand);
0086             randPoint = warmupPoints(:,randPointID);
0087             
0088             % Get a direction from the center point to the warmup point
0089             u = (randPoint-centerPoint);
0090             u = u/norm(u);
0091             
0092             % Figure out the distances to upper and lower bounds
0093             distUb = (model.ub - prevPoint);
0094             distLb = (prevPoint - model.lb);
0095             
0096             % Figure out if we are too close to a boundary
0097             validDir = ((distUb > dTol) & (distLb > dTol));
0098             %model.rxns(~validDir)
0099             
0100             % Zero out the directions that would bring us too close to the
0101             % boundary. This may cause problems.
0102             %u(~validDir) = 0;
0103             
0104             % Figure out positive and negative directions
0105             posDirn = find(u(validDir) > uTol);
0106             negDirn = find(u(validDir) < -uTol);
0107             
0108             % Figure out all the possible maximum and minimum step sizes
0109             maxStepTemp = distUb(validDir)./u(validDir);
0110             minStepTemp = -distLb(validDir)./u(validDir);
0111             maxStepVec = [maxStepTemp(posDirn);minStepTemp(negDirn)];
0112             minStepVec = [minStepTemp(posDirn);maxStepTemp(negDirn)];
0113             
0114             % Figure out the true max & min step sizes
0115             maxStep = min(maxStepVec);
0116             minStep = max(minStepVec);
0117             %fprintf('%f\t%f\n',minStep,maxStep);
0118             
0119             % Find new direction if we're getting too close to a constraint
0120             if (abs(minStep) < maxMinTol & abs(maxStep) < maxMinTol) | (minStep > maxStep)
0121                 fprintf('Warning %f %f\n',minStep,maxStep);
0122                 continue;
0123             end
0124             
0125             % Pick a rand out of list_of_rands and use it to get a random
0126             % step distance
0127             stepDist = randVector(stepCount)*(maxStep-minStep)+minStep;
0128             
0129             % Advance to the next point
0130             curPoint = prevPoint + stepDist*u;
0131             
0132             % Reproject the current point and go to the next step
0133             if mod(totalStepCount,10) == 0
0134                 if (full(max(max(abs(model.S*curPoint)))) > 1e-9)
0135                   curPoint = N*(N'*curPoint);
0136                 end
0137             end
0138             
0139             % Print out errors
0140             if (mod(totalStepCount,2000)==0)
0141               fprintf(fidErr,'%10.8f\t%10.8f\t',max(curPoint-model.ub),max(model.lb-curPoint));
0142             end
0143 
0144             timeElapsed = cputime-t0;
0145             
0146             % Print step information
0147             if (mod(totalStepCount,5000)==0)  
0148               timePerStep = timeElapsed/totalStepCount;
0149               fprintf('%d\t%d\t%d\t%8.2f\t%8.2f\n',i,pointCount,totalStepCount,timeElapsed/60,(totalCount-totalStepCount)*timePerStep/60);
0150             end
0151             
0152             overInd = find(curPoint > model.ub);
0153             underInd = find(curPoint < model.lb);
0154             
0155             if (any((model.ub-curPoint) < 0) | any((curPoint-model.lb) ...
0156                                                    < 0))
0157               curPoint(overInd) = model.ub(overInd);
0158               curPoint(underInd) = model.lb(underInd);
0159             
0160             end
0161             
0162             if (mod(totalStepCount,2000)==0)
0163               fprintf(fidErr,'%10.8f\n',full(max(max(abs(model.S*curPoint)))));
0164             end
0165             
0166             prevPoint = curPoint;
0167             stepCount = stepCount + 1;
0168             
0169             % Count the total number of steps
0170             totalStepCount = totalStepCount + 1;
0171             if mod(totalStepCount,50) == 0
0172                 waitbar(totalStepCount/totalCount,h);
0173             end        
0174             %recalculate the center point
0175             centerPoint = ((nWrmup+totalStepCount)*centerPoint + curPoint)/(nWrmup+totalStepCount+1);
0176 
0177             % Exit if time exceeded
0178             if (timeElapsed > maxTime)
0179                 points(:,pointCount) = curPoint;
0180                 file = [fileName '_' num2str(fileBaseNo+i) '.mat'];
0181                 save (file,'points');
0182                 save ACHR_last_point.mat curPoint
0183                 close(h);
0184                 return;
0185             end
0186             
0187         end % Steps per point
0188         
0189         % Add the current point to points
0190         points(:,pointCount) = curPoint;
0191        
0192         pointCount = pointCount + 1;
0193          
0194     end % Points per cycle
0195     
0196     % Save current points to a file
0197     file = [fileName '_' num2str(fileBaseNo+i) '.mat'];
0198     save (file,'points');
0199     
0200 end
0201 close(h);
0202 
0203 % Save last point for restart
0204 save ACHR_last_point.mat curPoint
0205 
0206 fclose(fidErr);

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