0001 function ACHRSampler(model,warmupPoints,fileName,nFiles,pointsPerFile,stepsPerPoint,initPoint,fileBaseNo,maxTime)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
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
0038 maxMinTol = 1e-9;
0039
0040 uTol = 1e-9;
0041
0042 dTol = 1e-14;
0043
0044
0045 [nRxns,nWrmup] = size(warmupPoints);
0046
0047
0048 centerPoint = mean(warmupPoints,2);
0049
0050
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
0073 points = zeros(nRxns,pointsPerFile);
0074
0075 pointCount = 1;
0076 while (pointCount <= pointsPerFile)
0077
0078
0079 randVector = rand(stepsPerPoint,1);
0080
0081 stepCount = 1;
0082 while (stepCount <= stepsPerPoint)
0083
0084
0085 randPointID = ceil(nWrmup*rand);
0086 randPoint = warmupPoints(:,randPointID);
0087
0088
0089 u = (randPoint-centerPoint);
0090 u = u/norm(u);
0091
0092
0093 distUb = (model.ub - prevPoint);
0094 distLb = (prevPoint - model.lb);
0095
0096
0097 validDir = ((distUb > dTol) & (distLb > dTol));
0098
0099
0100
0101
0102
0103
0104
0105 posDirn = find(u(validDir) > uTol);
0106 negDirn = find(u(validDir) < -uTol);
0107
0108
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
0115 maxStep = min(maxStepVec);
0116 minStep = max(minStepVec);
0117
0118
0119
0120 if (abs(minStep) < maxMinTol & abs(maxStep) < maxMinTol) | (minStep > maxStep)
0121 fprintf('Warning %f %f\n',minStep,maxStep);
0122 continue;
0123 end
0124
0125
0126
0127 stepDist = randVector(stepCount)*(maxStep-minStep)+minStep;
0128
0129
0130 curPoint = prevPoint + stepDist*u;
0131
0132
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
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
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
0170 totalStepCount = totalStepCount + 1;
0171 if mod(totalStepCount,50) == 0
0172 waitbar(totalStepCount/totalCount,h);
0173 end
0174
0175 centerPoint = ((nWrmup+totalStepCount)*centerPoint + curPoint)/(nWrmup+totalStepCount+1);
0176
0177
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
0188
0189
0190 points(:,pointCount) = curPoint;
0191
0192 pointCount = pointCount + 1;
0193
0194 end
0195
0196
0197 file = [fileName '_' num2str(fileBaseNo+i) '.mat'];
0198 save (file,'points');
0199
0200 end
0201 close(h);
0202
0203
0204 save ACHR_last_point.mat curPoint
0205
0206 fclose(fidErr);