solveQP Solves a quadratic fitting problem. model a model structure rxns either a cell array of reaction IDs, a logical vector with the same number of elements as reactions in the model, of a vector of indexes to fit to values the values to fit the fluxes to maxIter maximal number of iterations (opt, default 1000) restartIter run the fitting up to this many times in case it does not converge (opt, default 1) solution f Objective value x Primal stat Exit flag Usage: solution=solveQP(model,rxns,values,maxIter, restartIter) Rasmus Agren, 2013-02-13
0001 function solution=solveQP(model,rxns,values,maxIter, restartIter) 0002 % solveQP 0003 % Solves a quadratic fitting problem. 0004 % 0005 % model a model structure 0006 % rxns either a cell array of reaction IDs, a logical vector 0007 % with the same number of elements as reactions in the model, 0008 % of a vector of indexes to fit to 0009 % values the values to fit the fluxes to 0010 % maxIter maximal number of iterations (opt, default 1000) 0011 % restartIter run the fitting up to this many times in case it does 0012 % not converge (opt, default 1) 0013 % 0014 % solution 0015 % f Objective value 0016 % x Primal 0017 % stat Exit flag 0018 % 0019 % Usage: solution=solveQP(model,rxns,values,maxIter, restartIter) 0020 % 0021 % Rasmus Agren, 2013-02-13 0022 % 0023 0024 if nargin<4 0025 maxIter=1000; 0026 end 0027 0028 if nargin<5 0029 restartIter=1; 0030 end 0031 0032 %Check that it's feasible 0033 solution=solveLP(model); 0034 if isempty(solution.x) 0035 return; 0036 end 0037 0038 %Get the indexes of the fluxes to fit 0039 allIndexes=getIndexes(model,rxns,'rxns'); 0040 0041 f=zeros(numel(model.rxns),1); 0042 H=zeros(numel(model.rxns)); 0043 0044 %Get the fitIndexes for the experiment 0045 for j=1:numel(allIndexes) %Not that neat 0046 H(allIndexes(j),allIndexes(j))=2; 0047 end 0048 0049 f(allIndexes)=values.*-2; 0050 0051 %For a badly formulated problem it can occur that the solver stalls. 0052 %This can sometimes be fixed by trying to solve the problem again 0053 options=optimset('MaxIter',maxIter); 0054 for j=1:restartIter 0055 [x,fval,flag] = ... 0056 quadprog(H,f,[],[],model.S,model.b,model.lb,model.ub,[],options); 0057 if flag>0 0058 break; 0059 end 0060 end 0061 0062 solution.f=fval; 0063 solution.x=x; 0064 solution.stat=flag; 0065 end