solveQP

 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 (optional, default 1000)
   restartIter   run the fitting up to this many times in case it does
                 not converge (optional, default 1)

   solution
         f       Objective value
         x       Primal
         stat    Exit flag

 Usage: solution=solveQP(model,rxns,values,maxIter, restartIter)

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 (optional, default 1000)
0011 %   restartIter   run the fitting up to this many times in case it does
0012 %                 not converge (optional, 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 rxns=char(rxns);
0021 
0022 if nargin<4
0023     maxIter=1000;
0024 end
0025 
0026 if nargin<5
0027     restartIter=1;
0028 end
0029 
0030 %Check that it's feasible
0031 solution=solveLP(model);
0032 if isempty(solution.x)
0033     return;
0034 end
0035 
0036 %Get the indexes of the fluxes to fit
0037 allIndexes=getIndexes(model,rxns,'rxns');
0038 
0039 f=zeros(numel(model.rxns),1);
0040 H=zeros(numel(model.rxns));
0041 
0042 %Get the fitIndexes for the experiment
0043 for j=1:numel(allIndexes) %Not that neat
0044     H(allIndexes(j),allIndexes(j))=2;
0045 end
0046 
0047 f(allIndexes)=values.*-2;
0048 
0049 %For a badly formulated problem it can occur that the solver stalls. This
0050 %can sometimes be fixed by trying to solve the problem again
0051 options=optimset('MaxIter',maxIter);
0052 for j=1:restartIter
0053     [x,fval,flag] = ...
0054         quadprog(H,f,[],[],model.S,model.b,model.lb,model.ub,[],options);
0055     if flag>0
0056         break;
0057     end
0058 end
0059 
0060 solution.f=fval;
0061 solution.x=x;
0062 solution.stat=flag;
0063 end