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
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