0001 function [vs, output, v0] = C13ConfidenceInterval(v0, expdata, model, max_score, directions, majorIterationLimit)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 if nargin < 5
0018 directions = ones(size(v0,1),1);
0019 end
0020 if isempty(directions)
0021 directions = ones(size(v0,1),1);
0022 end
0023 t_start = clock;
0024 printLevel = 3;
0025 if nargin < 6
0026 majorIterationLimit = 1000;
0027 end
0028 diffInterval = 1e-5;
0029 feasibilityTolerance = max_score/20;
0030 logdirectory = strcat('temp', filesep);
0031
0032
0033
0034
0035
0036 denom = zeros(size(directions));
0037 if isnumeric(directions)
0038 if max(directions) == 1
0039 isratio = find(directions);
0040 else
0041 isratio = directions;
0042 end
0043 else
0044 isratio = zeros(size(directions));
0045 for i = 1:length(directions)
0046 if findstr(directions{i}, '/')
0047 [rxn1,rest] = strtok(directions{i}, '/');
0048 rxnID = findRxnIDs(model,rxn1);
0049 if rxnID == 0
0050 display('unable to process rxn from list');
0051 display(directions{i});
0052 return;
0053 else
0054 isratio(i) = -rxnID;
0055 end
0056 rxn2 = rest(2:end);
0057 rxnID = findRxnIDs(model,rxn2);
0058 if rxnID == 0
0059 display('unable to process rxn from list');
0060 display(rest(2:end));
0061 return;
0062 else
0063 denom(i) = rxnID;
0064 end
0065 else
0066 rxnID = findRxnIDs(model,directions{i});
0067 if rxnID == 0
0068 display('unable to process rxn from list');
0069 display(directions{i});
0070 return;
0071 else
0072 isratio(i) = rxnID;
0073 end
0074 end
0075 end
0076 end
0077
0078 numdirections = length(isratio);
0079 numpoints = size(v0,2);
0080
0081 numiterations = numdirections*numpoints*2;
0082
0083 x0 = model.N\v0;
0084 scores = zeros(numpoints,1);
0085
0086 tProb.user.expdata = expdata;
0087 tProb.user.model = model;
0088 for i = 1:numpoints
0089 scores(i) = errorComputation2(x0(:,i),tProb);
0090 end
0091
0092
0093 v0(:,scores> max_score) = fitC13Data(v0(:,scores > max_score),expdata,model, majorIterationLimit);
0094
0095 if ~isfield(model, 'N')
0096 model.N = null(model.S);
0097 end
0098
0099 x0 = model.N\v0;
0100
0101
0102 if (max(abs(model.S*v0))> 1e-6)
0103 display('v0 not quite in null space');
0104 pause;
0105 end
0106 if(max(abs(model.N*x0 - v0)) > 1e-6)
0107 display('null basis is weird');
0108 pause;
0109 end
0110
0111 Name = 't2';
0112 nalpha = size(model.N, 2);
0113
0114 x_L = -1000*ones(nalpha,1);
0115 x_U = 1000*ones(nalpha,1);
0116 [A, b_L, b_U] = defineLinearConstraints(model);
0117
0118 scores = zeros(numpoints,1);
0119
0120 for i = 1:numpoints
0121 scores(i) = errorComputation2(x0(:,i),tProb);
0122 end
0123 valid_index = scores < max_score + feasibilityTolerance;
0124 fprintf('found %d valid points\n', sum(valid_index));
0125
0126 x0_valid = x0(:,valid_index);
0127 x0_invalid = x0;
0128 scores_valid = scores(valid_index);
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138 checkedbefore = zeros(length(isratio),1);
0139 for i = 2:length(isratio)
0140 if(isratio(i) < 0)
0141 continue
0142 end
0143 d = objective_coefficient(isratio(i), model);
0144 for j = 1:i-1
0145 if(isratio(j) < 0)
0146 continue
0147 end
0148 dj = objective_coefficient(isratio(j), model);
0149 if max(abs(dj - d))< 1e-4
0150 checkedbefore(i) = j;
0151 break;
0152 elseif max(abs(dj + d))< 1e-4
0153 checkedbefore(i) = -j;
0154 break;
0155 end
0156 end
0157 end
0158
0159
0160 outputv = 222*ones(numiterations,1);
0161 outputexitflag = -222*ones(numiterations,1);
0162 outputfinalscore = -222*ones(numiterations,1);
0163 outputstruct = cell(numiterations,1);
0164
0165 csense = '';
0166 for mm = 1:length(b_L),csense(mm,1) = 'L';end
0167 for mm = 1:length(b_L),csense(mm+length(b_L),1) = 'G';end
0168
0169 fLowBnds = zeros(length(isratio), 2);
0170 for rxn = 1:length(isratio)
0171 for direction = -1:2:1
0172 if isratio(rxn) < 0
0173 ration = objective_coefficient(-isratio(rxn),model);
0174 ratiod = objective_coefficient(denom(rxn),model);
0175
0176
0177 Result = solveCobraNLP(...
0178 struct('lb', x_L, 'ub', x_U,...
0179 'name', Name,...
0180 'A', A,...
0181 'b_L', b_L, 'b_U', b_U,...
0182 'objFunction', 'ratioScore', 'g', 'ratioScore_grad',...
0183 'userParams', struct(...
0184 'ration', direction*ration, 'ratiod', ratiod,...
0185 'diff_interval', diffInterval,'useparfor', true)...
0186 ),...
0187 'printLevel', 1, ...
0188 'iterationLimit', 1000);
0189 else
0190 d = objective_coefficient(isratio(rxn),model);
0191 Result = solveCobraLP(...
0192 struct('A', [A;A],'b',[b_U;b_L],'csense', csense, ...
0193 'c', direction*d, ...
0194 'lb', x_L,'ub', x_U, ...
0195 'osense', 1),...
0196 'feasTol',1e-7,'optTol',1e-7);
0197 if Result.stat ~= 1
0198 Result
0199 pause
0200 end
0201 end
0202 fLowBnds(rxn, (direction+3)/2) = direction*Result.obj;
0203 end
0204 end
0205
0206 if ~exist (logdirectory, 'dir')
0207 if ~mkdir(logdirectory)
0208 display('unable to create logdirectory');
0209 return;
0210 end
0211 end
0212 clear d
0213
0214 parfor itnum = 1:numiterations
0215 if exist('ttt.txt', 'file')
0216 fprintf('quitting due to file found\n');
0217 continue;
0218 end
0219 [rxn, direction, point] = getValues(itnum, numpoints);
0220
0221
0222
0223 if checkedbefore(rxn) ~= 0
0224 continue;
0225 end
0226
0227
0228 fLowBnd = fLowBnds(rxn, (direction+3)/2);
0229 fprintf('reaction %d of %d, direction %d, lowerbound %f point %d of %d\n', rxn ,length(isratio), direction, fLowBnd, point, numpoints);
0230
0231 if isratio(rxn) > 0
0232 di = objective_coefficient(isratio(rxn),model);
0233 obj1 = di'*x0_valid;
0234 else
0235 rationi = objective_coefficient(-isratio(rxn),model);
0236 ratiodi = objective_coefficient(denom(rxn),model);
0237 obj1 = (rationi'*x0_valid) ./ (ratiodi'*x0_valid);
0238 end
0239
0240 if(any(abs(obj1-fLowBnd)<.0001))
0241 display('short circuiting');
0242 [nil, min_index] = min(obj1);
0243 outputv(itnum,1) = fLowBnd;
0244 outputexitflag(itnum,1) = 111;
0245 outputfinalscore(itnum,1) = scores_valid(min_index);
0246
0247 else
0248 xinitial = x0_invalid(:,point);
0249 if isratio(rxn) > 0
0250 NLPsolution = solveCobraNLP(...
0251 struct('x0', xinitial, ...
0252 'lb', x_L, 'ub', x_U,...
0253 'name', Name,...
0254 'A', A, 'b_L', b_L, 'b_U', b_U,...
0255 'd', 'errorComputation2', 'dd', 'errorComputation2_grad',...
0256 'd_L', 0, 'd_U', max_score,...
0257 'c', di*direction, ...
0258 'userParams', struct(...
0259 'expdata', expdata,'model', model,'max_error', max_score,...
0260 'diff_interval', diffInterval,'useparfor', true)...
0261 ),...
0262 'printLevel', printLevel, ...
0263 ...
0264 'iterationLimit', majorIterationLimit, ...
0265 'logFile', strcat(logdirectory, 'ci_', num2str(rxn),'x',num2str(direction),'x', point, '.txt'));
0266 else
0267 NLPsolution = solveCobraNLP(...
0268 struct('x0', xinitial, ...
0269 'lb', x_L, 'ub', x_U,...
0270 'name', Name,...
0271 'A', A, 'b_L', b_L, 'b_U', b_U,...
0272 'd', 'errorComputation2', 'dd', 'errorComputation2_grad',...
0273 'd_L', 0, 'd_U', max_score,...
0274 'objFunction', 'ratioScore', 'g', 'ratioScore_grad',...
0275 'userParams', struct(...
0276 'expdata', expdata,'model', model,'max_error', max_score,...
0277 'ration', direction*rationi,...
0278 'ratiod', ratiodi,...
0279 'diff_interval', diffInterval,'useparfor', false)...
0280 ),...
0281 'printLevel', printLevel, ...
0282 ...
0283 'iterationLimit', majorIterationLimit, ...
0284 'logFile', strcat(logdirectory, 'ci_', num2str(rxn),'x',num2str(direction),'x', point, '.txt'));
0285 end
0286
0287 tscore = errorComputation2(NLPsolution.full, tProb);
0288 tbest = NLPsolution.obj;
0289
0290 fprintf('reaction %d (%d), x %d; x=%f (%f); score=%f (%f)\n', rxn, length(isratio),direction, tbest,fLowBnd, tscore, max_score)
0291
0292 outputv(itnum,1) = direction*tbest;
0293 outputexitflag(itnum,1) = NLPsolution.origStat;
0294 outputfinalscore(itnum,1) = tscore;
0295 outputstruct{itnum,1} = NLPsolution;
0296 end
0297 end
0298
0299 for itnum = 1:numiterations
0300 [rxn, direction, point] = getValues(itnum,numpoints);
0301 if direction == 1;
0302 output.minv(rxn,point) = outputv(itnum);
0303 output.minexitflag(rxn,point) = outputexitflag(itnum);
0304 output.minfinalscore(rxn,point) = outputfinalscore(itnum);
0305 output.minstruct(rxn,point) = outputstruct(itnum);
0306 else
0307 output.maxv(rxn,point) = outputv(itnum);
0308 output.maxexitflag(rxn,point) = outputexitflag(itnum);
0309 output.maxfinalscore(rxn,point) = outputfinalscore(itnum);
0310 output.maxstruct(rxn,point) = outputstruct(itnum);
0311 end
0312 end
0313
0314 for i = 1:length(isratio)
0315 if checkedbefore(i) > 0
0316 output.minv(i) = output.minv(checkedbefore(i));
0317 output.maxv(i) = output.maxv(checkedbefore(i));
0318 output.minexitflag(i) = output.minexitflag(checkedbefore(i));
0319 output.maxexitflag(i) = output.maxexitflag(checkedbefore(i));
0320 output.minfinalscore(i) = output.minfinalscore(checkedbefore(i));
0321 output.maxfinalscore(i) = output.maxfinalscore(checkedbefore(i));
0322 output.minstruct(i) = output.minstruct(checkedbefore(i));
0323 output.maxstruct(i) = output.maxstruct(checkedbefore(i));
0324 elseif checkedbefore(i) < 0
0325 output.minv(i) = -output.maxv(-checkedbefore(i));
0326 output.maxv(i) = -output.minv(-checkedbefore(i));
0327 output.minexitflag(i) = output.maxexitflag(-checkedbefore(i));
0328 output.maxexitflag(i) = output.minexitflag(-checkedbefore(i));
0329 output.minfinalscore(i) = output.maxfinalscore(-checkedbefore(i));
0330 output.maxfinalscore(i) = output.minfinalscore(-checkedbefore(i));
0331 output.minstruct(i) = output.maxstruct(-checkedbefore(i));
0332 output.maxstruct(i) = output.minstruct(-checkedbefore(i));
0333 end
0334 end
0335
0336 vs = zeros(length(isratio), 2);
0337 for i = 1:length(isratio)
0338 validindex = output.minfinalscore(i,:) < max_score + feasibilityTolerance;
0339 if any(validindex)
0340 vs(i,1) = min(output.minv(i,validindex));
0341 else
0342 vs(i,1) = 222;
0343 end
0344 validindex = output.maxfinalscore(i,:) < max_score + feasibilityTolerance;
0345 if any(validindex)
0346 vs(i,2) = max(output.maxv(i,validindex));
0347 else
0348 vs(i,2) = -222;
0349 end
0350 end
0351
0352 elapsed_time = etime(clock, t_start)
0353 return;
0354
0355
0356
0357
0358
0359
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372
0373 function [rxn, direction, point] = getValues(index, numpoints)
0374
0375
0376
0377
0378 index = index - 1;
0379 point = mod(index, numpoints);
0380 index = index - point;
0381 index = index/numpoints;
0382 direction = mod(index,2);
0383 index = index - direction;
0384 index = index/2;
0385 rxn = index;
0386
0387 point = point +1;
0388 direction = direction*2-1;
0389 rxn = rxn+1;
0390 return;
0391
0392
0393
0394 function [d] = objective_coefficient(i, model)
0395 d = zeros(length(model.lb),1);
0396 d(i) = 1;
0397 if (model.match(i))
0398 d(model.match(i)) = -1;
0399 end
0400 d = (d'*model.N)';
0401 return