C13ConfidenceInterval

PURPOSE ^

v0 - set of flux vectors to be used as initial guesses. They may be

SYNOPSIS ^

function [vs, output, v0] = C13ConfidenceInterval(v0, expdata, model, max_score, directions, majorIterationLimit)

DESCRIPTION ^

 v0 - set of flux vectors to be used as initial guesses.  They may be
 valid or not.
 expdata - experimental data.
 model - The standard model.  Additional field .N (= null(S)) should also
 be provided.  This is a basis of the flux space.
 max score - maximum allowable data fit error.  
 directions (optional) - ones and zeros of which reactions to compute (size = n
 x 1).  
   OR
 numbers of reactions to use  aka.  [1;5;7;8;200]
   OR
 reaction strings  aka.  {'GPK', 'PGL'}.  Ratios are possible with this
 input only.  Default = [] meaning do FVA with no ratios.
 majorIterationLimit (optional) - default = 10000

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [vs, output, v0] = C13ConfidenceInterval(v0, expdata, model, max_score, directions, majorIterationLimit)
0002 
0003 % v0 - set of flux vectors to be used as initial guesses.  They may be
0004 % valid or not.
0005 % expdata - experimental data.
0006 % model - The standard model.  Additional field .N (= null(S)) should also
0007 % be provided.  This is a basis of the flux space.
0008 % max score - maximum allowable data fit error.
0009 % directions (optional) - ones and zeros of which reactions to compute (size = n
0010 % x 1).
0011 %   OR
0012 % numbers of reactions to use  aka.  [1;5;7;8;200]
0013 %   OR
0014 % reaction strings  aka.  {'GPK', 'PGL'}.  Ratios are possible with this
0015 % input only.  Default = [] meaning do FVA with no ratios.
0016 % majorIterationLimit (optional) - default = 10000
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;  %max number of iterations
0027 end
0028 diffInterval = 1e-5;         %gradient step size.
0029 feasibilityTolerance = max_score/20; % how close you need to be to the max score.
0030 logdirectory = strcat('temp', filesep);
0031 
0032 % isratio(i):  >0 indicates not a ratio and points to reaction#
0033 %              <0 indicates -numerator of ratio fraction.  In this case,
0034 %              denom(i) stores the denominator
0035 %
0036 denom = zeros(size(directions));
0037 if isnumeric(directions) %cannot be an actual ratio
0038     if max(directions) == 1
0039         isratio = find(directions);
0040     else
0041         isratio = directions;
0042     end
0043 else % might be a ratio.  Gotta process strings
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; % total number of iterations.
0082 
0083 x0 = model.N\v0; % back substitute
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 % fit points if they are not currently feasible
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; % back substitute
0100 
0101 % safety check:
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 % compute scores for all points.
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 % pre-compute unnecesary directions.
0132 % if checkedbefore(i) ~= 0 then direction i is redundant
0133 % if checkedbefore(i) = j then direction i and j are identical and do not
0134 % need to be recomputed.
0135 % checkedbefore(i) = j < 0 means that direction j is the same as i except
0136 % for a sign switch.
0137 
0138 checkedbefore = zeros(length(isratio),1);
0139 for i = 2:length(isratio)
0140     if(isratio(i) < 0) % meaning it actually IS a ratio and no simplification possible
0141         continue
0142     end
0143     d = objective_coefficient(isratio(i), model);
0144     for j = 1:i-1
0145         if(isratio(j) < 0) % meaning it actually IS a ratio and no simplification possible
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 % initialize variables;
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); %initialize but fill in later.
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             % in case RXN is a ratio
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,... % set direction here too.
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; % fill in
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 %iterate through directions
0214 parfor itnum = 1:numiterations
0215     if exist('ttt.txt', 'file') % abort w/o crashing if file 'ttt.txt' found in current directory
0216         fprintf('quitting due to file found\n');
0217         continue;
0218     end
0219     [rxn, direction, point] = getValues(itnum, numpoints); % translate itnum to rxn, direction and point
0220     % direction == 1 means minimize.  direction == -1 means maximize
0221     % (opposite of what you might think.
0222     
0223     if checkedbefore(rxn) ~= 0 %if this reaction maps to a previous reaction, we can skip
0224         continue;
0225     end
0226     
0227     
0228     fLowBnd = fLowBnds(rxn, (direction+3)/2); %get the absolute bound in the space w/o regards to C13 constraints.
0229     fprintf('reaction %d of %d, direction %d, lowerbound %f point %d of %d\n', rxn ,length(isratio), direction, fLowBnd, point, numpoints);
0230      % short circuit if x0 already close to a bound.
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; % multiply by direction to correct sign.
0244         outputexitflag(itnum,1) = 111;
0245         outputfinalscore(itnum,1) = scores_valid(min_index);
0246 
0247     else % gotta actually do the computation.
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, ... % direction of optimization
0258                 'userParams', struct(...
0259                      'expdata', expdata,'model', model,'max_error', max_score,...
0260                      'diff_interval', diffInterval,'useparfor', true)...
0261                 ),...
0262                 'printLevel', printLevel, ...
0263                 ...%'intTol', 1e-7, ...
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                 ...%'intTol', 1e-7, ...
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; % multiply by direction to correct sign.
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 %short circuit if seen before.
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 % function [index] = getIndex(rxn, direction, point, numpoints)
0358 % % point goes from 1 .. NUMPOINTS
0359 % % rxn goes from 1 .. NUMRXNS
0360 % % direction is -1 or 1
0361 % % index goes from 1 to NUMPOINTS*NUMRXNS*2
0362 %
0363 % rxn = rxn - 1;
0364 % point = point -1;
0365 % direction = (direction + 1)/2;
0366 %
0367 %
0368 % index = rxn*numpoints*2 + direction*numpoints + point;
0369 % index = index+1;
0370 %
0371 % return;
0372 
0373 function [rxn, direction, point] = getValues(index, numpoints)
0374 % point goes from 1 .. NUMPOINTS
0375 % rxn goes from 1 .. NUMRXNS
0376 % direction is -1 or 1
0377 % index goes from 1 to NUMPOINTS*NUMRXNS*2
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; % remap to 1..NUMPOINTS
0388 direction = direction*2-1; % remap to -1,1
0389 rxn = rxn+1; % remap to 1 .. number of rxns;
0390 return;
0391 
0392 % function that returns the proper objective coefficient for each reaction
0393 % takes into account the reversibility of reactinos etc.
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)'; % transform to null space;
0401 return

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