0001 function [sampleStructOut, mixedFrac] = gpSampler(sampleStruct, nPoints, bias, maxTime, maxSteps, threads, nPointsCheck)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052 sampleStructOut = 0;
0053 mixedFrac = 1;
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;
0068 end
0069 if (nargin < 5) || isempty(maxSteps)
0070
0071 maxSteps = 1e10;
0072 else
0073
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
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
0113
0114 [A, b, csense, lb, ub] = deal(sampleStruct.A, sampleStruct.b, sampleStruct.csense, sampleStruct.lb, sampleStruct.ub);
0115
0116
0117
0118
0119
0120
0121
0122
0123
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
0155 Bnew = Bnew(1:rAnew,:);
0156 Dnew = Dnew(1:rCnew,:);
0157
0158
0159 if find(Bnew ~= 0)
0160 offset = Anew\Bnew;
0161 else
0162 offset = zeros(size(Anew,2), 1);
0163 end
0164
0165
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
0195 if (~ isfield(sampleStruct, 'warmupPts') )
0196 fprintf('Generating warmup points\n');
0197
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
0211 end
0212
0213
0214 fprintf('Sampling\n');
0215 if(maxTime > 0 && maxSteps > 0)
0216 if threads < 0
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
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
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
0249
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
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
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
0312
0313
0314
0315
0316 LPproblem.c = rand(dimX, 1)-.5;
0317 validFlag = true;
0318
0319 for maxMin = [1, -1]
0320
0321 if i <= dimX
0322
0323 LPproblem.c(i) = 5000;
0324 end
0325 LPproblem.osense = maxMin;
0326
0327
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
0344 x(x > LPproblem.ub) = LPproblem.ub(x > LPproblem.ub);
0345 x(x < LPproblem.lb) = LPproblem.lb(x < LPproblem.lb);
0346
0347
0348
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
0359 i = i+1;
0360 end
0361 end
0362 centerPoint = mean(warmupPts,2);
0363
0364
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
0372
0373
0374
0375
0376 return;
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410