Home > RAVEN > solveLP.m

solveLP

PURPOSE ^

solveLP

SYNOPSIS ^

function [solution hsSolOut]=solveLP(model,minFlux,params,hsSol)

DESCRIPTION ^

 solveLP
   Solves a linear programming problem

   model         a model structure
   minFlux       determines if a second optimization should be performed
                 in order to get rid of loops in the flux distribution
                 0: no such optimization is performed
                 1: the sum of abs(fluxes) is minimized. This is the
                 fastest way of getting rid of loops
                 2: the square of fluxes is minimized. This tends to
                 distribute fluxes across iso-enzymes, which results in a
                 larger number of reactions being used
                 3: the number of fluxes is minimized. This can result
                 in the flux distributions that are the easiest to
                 interpret. Note that this optimization can be very slow
                 (opt, default 0)
   params        parameter structure as used by getMILPParams. Only used when
                 minFlux==3 (opt)
   hsSol         hot-start solution for the LP solver. This can
                 significantly speed up the process if many similar
                 optimization problems are solved iteratively. Only used if
                 minFlux is 0 or 1 (opt)

   solution
         f       objective value
         x       primal (flux distribution)
         stat    exit flag 
                 1: the optimization terminated successfully
                -1: no feasible solution found by Mosek
                -2: the solution violated the bounds
   hsSolOut      solution to be used as hot-start solution (see the input
                 parameters). Only used if minFlux is 0 or 1

   Usage: [solution hsSolOut]=solveLP(model,minFlux,params,hsSol)

   Rasmus Agren, 2013-04-22

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [solution hsSolOut]=solveLP(model,minFlux,params,hsSol)
0002 % solveLP
0003 %   Solves a linear programming problem
0004 %
0005 %   model         a model structure
0006 %   minFlux       determines if a second optimization should be performed
0007 %                 in order to get rid of loops in the flux distribution
0008 %                 0: no such optimization is performed
0009 %                 1: the sum of abs(fluxes) is minimized. This is the
0010 %                 fastest way of getting rid of loops
0011 %                 2: the square of fluxes is minimized. This tends to
0012 %                 distribute fluxes across iso-enzymes, which results in a
0013 %                 larger number of reactions being used
0014 %                 3: the number of fluxes is minimized. This can result
0015 %                 in the flux distributions that are the easiest to
0016 %                 interpret. Note that this optimization can be very slow
0017 %                 (opt, default 0)
0018 %   params        parameter structure as used by getMILPParams. Only used when
0019 %                 minFlux==3 (opt)
0020 %   hsSol         hot-start solution for the LP solver. This can
0021 %                 significantly speed up the process if many similar
0022 %                 optimization problems are solved iteratively. Only used if
0023 %                 minFlux is 0 or 1 (opt)
0024 %
0025 %   solution
0026 %         f       objective value
0027 %         x       primal (flux distribution)
0028 %         stat    exit flag
0029 %                 1: the optimization terminated successfully
0030 %                -1: no feasible solution found by Mosek
0031 %                -2: the solution violated the bounds
0032 %   hsSolOut      solution to be used as hot-start solution (see the input
0033 %                 parameters). Only used if minFlux is 0 or 1
0034 %
0035 %   Usage: [solution hsSolOut]=solveLP(model,minFlux,params,hsSol)
0036 %
0037 %   Rasmus Agren, 2013-04-22
0038 %
0039 
0040 if nargin<2
0041     minFlux=0;
0042 end
0043 if nargin<3
0044     params.relGap=0.4;
0045 end
0046 if nargin<4
0047     hsSol=[];
0048 end
0049 hsSolOut=[];
0050 
0051 %Ignore the hot-start if the previous solution wasn't feasible
0052 if isfield(hsSol,'prosta')
0053    if strfind(hsSol.prosta,'INFEASIBLE')
0054         hsSol=[];
0055    end
0056 end
0057     
0058 % Setup the problem to feed to MOSEK.
0059 prob=[];
0060 prob.c=model.c*-1;
0061 prob.a=model.S;
0062 prob.blc=model.b(:,1);
0063 %If model.b has two column, then they are for lower/upper bound on the RHS
0064 prob.buc=model.b(:,min(size(model.b,2),2));
0065 prob.blx=model.lb;
0066 prob.bux=model.ub;
0067 
0068 %If hot-start should be used
0069 if ~isempty(hsSol)
0070    prob.sol.bas=hsSol;
0071    params.MSK_IPAR_SIM_HOTSTART=1;
0072 end
0073 
0074 %Use MSK_OPTIMIZER_FREE_SIMPLEX
0075 params.MSK_IPAR_OPTIMIZER=7;
0076 
0077 [crap,res] = mosekopt('minimize echo(0)',prob,getMILPParams(params));
0078 if isfield(res.sol,'bas') && ~isempty(hsSol)
0079     if isfield(res.sol.bas,'prosta')
0080        if strfind(res.sol.bas.prosta,'INFEASIBLE')
0081             %If the problem was infeasible using hot-start it is often
0082             %possible to re-solve it without hot-start and get a feasible
0083             %solution
0084             prob.sol=rmfield(prob.sol,'bas');
0085             [crap,res] = mosekopt('minimize echo(0)',prob,getMILPParams(params));
0086        end
0087     end
0088 end
0089 if res.rcode==1001
0090     throw(MException('','The Mosek licence has expired'));
0091 end
0092 if res.rcode==1010
0093     throw(MException('','The Mosek licence used only supports small problems (up to 300 variables). Have you requested the correct licence?'));
0094 end
0095 
0096 %Check that the problem was solved and construct the output structure
0097 if isfield(res,'sol')
0098     if isfield(res.sol,'bas')
0099         solution.x=res.sol.bas.xx;
0100         if minFlux<=1;
0101             hsSolOut=res.sol.bas;
0102         end
0103     else
0104         %Interior-point
0105         solution.x=res.sol.itr.xx;
0106     end
0107 else
0108     solution.x=[];
0109 end
0110 
0111 if length(prob.c)==length(solution.x)
0112     solution.f=prob.c'*solution.x; 
0113 else
0114     solution.f=[];
0115 end
0116 
0117 if ~isempty(res) && ~isempty(solution.x)
0118     if res.rcode~=0
0119         solution.stat=res.rcode;
0120     else
0121         solution.stat=1;
0122     end
0123 else
0124     solution.stat=-1; 
0125 end
0126 
0127 %If either the square, the number, or the sum of fluxes should be minimized
0128 %then the objective function value should be fixed before another
0129 %optimization. It is not correct to fix the reactions which participate in
0130 %the objective function to their values in solution.x, as there can be
0131 %multiple solutions with the same objective function value. In addition, this
0132 %approach could result in numerical issues when several fluxes are fixed.
0133 %Instead a new "fake metabolite" is added to the problem. This metabolite
0134 %is produced by each reaction with the stoichiometry that reaction has in
0135 %the objective function. The equality constraint of that "fake metabolite"
0136 %is then set to be at least as good as the objective function value.
0137 if minFlux~=0 && ~isempty(solution.x)
0138     model.S=[model.S;prob.c'];
0139     model.mets=[model.mets;'TEMP'];
0140     
0141     %If the constraint on the objective function value is exact there is a
0142     %larger risk of numerical errors. However for the quadratic fitting
0143     %intervals are not allowed
0144     if minFlux~=2
0145         if size(model.b,2)==1
0146             model.b=[model.b model.b];
0147         end
0148         if solution.f<0
0149             model.b=[model.b;[-inf solution.f*0.999999]];
0150         else
0151             model.b=[model.b;[-inf solution.f*1.000001]];
0152         end
0153     else
0154         model.b=[model.b;ones(1,size(model.b,2))*solution.f];
0155     end
0156     
0157     switch minFlux
0158         %The sum of fluxes should be minimized
0159         case 1
0160             %Convert the model to the irreversible format
0161             revRxns=find(model.rev);
0162             if ~isempty(revRxns)
0163                iModel=convertToIrrev(model);
0164             else
0165                 iModel=model;
0166             end
0167             
0168             %Minimize all fluxes
0169             iModel.c(:)=-1;
0170             sol=solveLP(iModel);
0171             
0172             %Map back to reversible fluxes
0173             if ~isempty(sol.x)
0174                 solution.x=sol.x(1:numel(model.c));
0175                 solution.x(revRxns)=solution.x(revRxns)-sol.x(numel(model.c)+1:end);
0176             else
0177                 solution.x=[];
0178                 solution.f=[];
0179                 solution.stat=-1;
0180             end
0181         %The square of fluxes should be minimized. This only works when
0182         %there is no interval on the mass balance constraints (model.b is a
0183         %vector)
0184         case 2
0185 %         if size(model.b,2)==1
0186 %             qsol=solveQP(model,model.rxns,zeros(numel(model.lb),1));
0187 %             %There is a problem that the solver seldom converges totally in this
0188 %             %kind of minimization. Print a warning but use the fluxes
0189 %             if any(qsol.x)
0190 %                 solution.x=qsol.x;
0191 %                 if qsol.stat==-1
0192 %                     fprintf('WARNING: The quadratic fitting did not converge\n');
0193 %                 end
0194 %             else
0195 %                 fprintf('WARNING: Could not solve the problem of minimizing the square of fluxes. Uses output from linear program\n');
0196 %             end
0197 %         else
0198 %             fprintf('WARNING: Cannot minimize square of fluxes when size(model.b,2)==2. Uses output from linear program\n');
0199 %         end
0200         fprintf('WARNING: Quadratic solver currently not working. Uses output from linear program\n');
0201         %The number of fluxes should be minimized
0202         case 3
0203             [qx,I]=getMinNrFluxes(model,model.rxns,params);
0204             qx(~I)=0;
0205             if any(I)
0206                 solution.x=qx;
0207             else
0208                fprintf('WARNING: Could not solve the problem of minimizing the number of fluxes. Uses output from linear program\n');
0209             end
0210     end
0211 end
0212 
0213 %Make a final check that no bounds were violated. This information could
0214 %probably be retrieved from the res-structure, but this is just to make
0215 %sure
0216 if ~isempty(solution.x)
0217     badSolution=false;
0218     if any((model.ub+10^-6)<solution.x | (model.lb-10^-6)>solution.x)
0219        badSolution=true;
0220     else
0221         %Also check that the mass constraints hold
0222         massBal=prob.a*solution.x;
0223         if size(model.b,2)==1
0224             if any((massBal+10^-6)<model.b(1:numel(massBal),1) | (massBal-10^-6)>model.b(1:numel(massBal),1))
0225                 badSolution=true;
0226             end
0227         else
0228            %This means that there are upper and lower constraints on the
0229            %mass balances
0230            %The reason for not always taking the full model.b array is
0231            %because a fake metabolite is added if minimizing for
0232            %square/number/sum of fluxes
0233            if any((massBal+10^-6)<model.b(1:numel(massBal),1) | (massBal-10^-6)>model.b(1:numel(massBal),2))
0234                badSolution=true;
0235            end
0236         end
0237     end
0238     if badSolution==true
0239         solution.f=[];
0240         solution.stat=-2;
0241         solution.x=[];
0242     end
0243 end
0244 end

Generated on Tue 23-Apr-2013 15:18:37 by m2html © 2005