0001 function ACBSampler(model,warmupPoints,fileName,nFiles,pointsPerFile,nMixPts,nWarmupNeeded,saveMatFlag,biasOpt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025 if (nargin < 7)
0026
0027
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
0042 maxMinTol = 1e-9;
0043
0044 uTol = 1e-9;
0045
0046 dTol = 1e-10;
0047
0048
0049 [nRxns,nWarmup] = size(warmupPoints);
0050
0051
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
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
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
0092 u = (randPoint2-randPoint1);
0093 u = u/norm(u);
0094
0095
0096
0097
0098
0099
0100 distUb = (model.ub - randPoint1);
0101 distLb = (randPoint1 - model.lb);
0102
0103
0104 validDir = ((distUb > dTol) & (distLb > dTol));
0105
0106
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
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
0121 maxStep = min(maxStepVec);
0122 minStep = max(minStepVec);
0123
0124
0125 boundaryPoint1 = randPoint1 + minStep*u;
0126 boundaryPoint2 = randPoint1 + maxStep*u;
0127
0128
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
0137 centerPoint = (boundaryPoint1 + boundaryPoint2)/2.0;
0138
0139
0140
0141
0142 if (biasFlag)
0143
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
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
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
0194 totalPointCount = totalPointCount + 1;
0195
0196
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
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
0219 if (saveMatFlag)
0220 file = [fileName '_' num2str(fileNo) '.mat'];
0221 save (file,'points');
0222 else
0223 fclose(fid);
0224 end
0225
0226 end
0227
0228 fclose(fidErr);
0229
0230
0231