solveCobraCPLEX

PURPOSE ^

[solution,LPProblem]=solveCobraCPLEX(model,printLevel,basisReuse,conflictResolve,contFunctName,minNorm)

SYNOPSIS ^

function [solution,model]=solveCobraCPLEX(model,printLevel,basisReuse,conflictResolve,contFunctName,minNorm)

DESCRIPTION ^

 [solution,LPProblem]=solveCobraCPLEX(model,printLevel,basisReuse,conflictResolve,contFunctName,minNorm)
 call CPLEX to solve an LP or QP problem using the matlab API to cplex written by ILOG

INPUT
 Model Structure containing the following fields describing the LP
 problem to be solved
  A or S       m x n LHS matrix
  b            m x 1 RHS vector
  c            n x 1 Objective coeff vector
  lb           n x 1 Lower bound vector
  ub           n x 1 Upper bound vector
  osense       scalar Objective sense (-1 max, +1 min)

OPTIONAL INPUT
 model.rxns    cell array of reaction abbreviations (necessary for
                   making a readable confilict resolution file).
 model.csense  Constraint senses, a string containting the constraint sense for
                   each row in A ('E', equality, 'G' greater than, 'L' less than).

 model.LPBasis Basis from previous solution of similar LP problem.
                   See basisReuse

 PrintLevel    Printing level in the CPLEX m-file and CPLEX C-interface.
               = 0    Silent
               = 1    Warnings and Errors
               = 2    Summary information (Default)
               = 3    More detailed information
               > 10   Pause statements, and maximal printing (debug mode)

 basisReuse = 0   Use this for one of soluion of an LP (Default)
            = 1   Returns a basis for reuse in the next LP
                  i.e. outputs model.LPBasis

 conflictResolve  = 0   (Default)
                  = 1   If LP problem is proven to be infeasible by CPLEX,
                        it will print out a 'conflict resolution file',
                        which indicates the irreducible infeasible set of
                        equaltiy & inequality constraints that together,
                        combine to make the problem infeasible. This is
                        useful for debugging an LP problem if you want to
                        try to resolve a constraint conflict

 contFunctName        = [] Use all default CLPEX control parameters, (Default)
                      = someString e.g. 'someFunctionName'
                        uses the user specified control parameters defined
                        in someFunctionName.m
                       (see template function CPLEXParamSet for details).
                      = cpxControl structure (output from a file like CPLEXParamSet.m)

 minNorm       {(0), 1 , n x 1 vector} If not zero then, minimise the Euclidean length
               of the solution to the LP problem. Gives the same objective,
               but minimises the square of flux. minNorm ~1e-6 should be
               high enough for regularisation yet keep the same objective

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [solution,model]=solveCobraCPLEX(model,printLevel,basisReuse,conflictResolve,contFunctName,minNorm)
0002 % [solution,LPProblem]=solveCobraCPLEX(model,printLevel,basisReuse,conflictResolve,contFunctName,minNorm)
0003 % call CPLEX to solve an LP or QP problem using the matlab API to cplex written by ILOG
0004 %
0005 %INPUT
0006 % Model Structure containing the following fields describing the LP
0007 % problem to be solved
0008 %  A or S       m x n LHS matrix
0009 %  b            m x 1 RHS vector
0010 %  c            n x 1 Objective coeff vector
0011 %  lb           n x 1 Lower bound vector
0012 %  ub           n x 1 Upper bound vector
0013 %  osense       scalar Objective sense (-1 max, +1 min)
0014 %
0015 %OPTIONAL INPUT
0016 % model.rxns    cell array of reaction abbreviations (necessary for
0017 %                   making a readable confilict resolution file).
0018 % model.csense  Constraint senses, a string containting the constraint sense for
0019 %                   each row in A ('E', equality, 'G' greater than, 'L' less than).
0020 %
0021 % model.LPBasis Basis from previous solution of similar LP problem.
0022 %                   See basisReuse
0023 %
0024 % PrintLevel    Printing level in the CPLEX m-file and CPLEX C-interface.
0025 %               = 0    Silent
0026 %               = 1    Warnings and Errors
0027 %               = 2    Summary information (Default)
0028 %               = 3    More detailed information
0029 %               > 10   Pause statements, and maximal printing (debug mode)
0030 %
0031 % basisReuse = 0   Use this for one of soluion of an LP (Default)
0032 %            = 1   Returns a basis for reuse in the next LP
0033 %                  i.e. outputs model.LPBasis
0034 %
0035 % conflictResolve  = 0   (Default)
0036 %                  = 1   If LP problem is proven to be infeasible by CPLEX,
0037 %                        it will print out a 'conflict resolution file',
0038 %                        which indicates the irreducible infeasible set of
0039 %                        equaltiy & inequality constraints that together,
0040 %                        combine to make the problem infeasible. This is
0041 %                        useful for debugging an LP problem if you want to
0042 %                        try to resolve a constraint conflict
0043 %
0044 % contFunctName        = [] Use all default CLPEX control parameters, (Default)
0045 %                      = someString e.g. 'someFunctionName'
0046 %                        uses the user specified control parameters defined
0047 %                        in someFunctionName.m
0048 %                       (see template function CPLEXParamSet for details).
0049 %                      = cpxControl structure (output from a file like CPLEXParamSet.m)
0050 %
0051 % minNorm       {(0), 1 , n x 1 vector} If not zero then, minimise the Euclidean length
0052 %               of the solution to the LP problem. Gives the same objective,
0053 %               but minimises the square of flux. minNorm ~1e-6 should be
0054 %               high enough for regularisation yet keep the same objective
0055 
0056 %OUTPUT
0057 % solution Structure containing the following fields describing a LP
0058 % solution
0059 %  full         Full LP solution vector
0060 %  obj          Objective value
0061 %  rcost        Lagrangian multipliers to the simple inequalties (Reduced costs)
0062 %  dual         Lagrangian multipliers to the equalities
0063 %  nInfeas      Number of infeasible constraints
0064 %  sumInfeas    Sum of constraint violation
0065 %  stat         COBRA Standardized solver status code:
0066 %               1   Optimal solution
0067 %               2   Unbounded solution
0068 %               0   Infeasible
0069 %               -1  No solution reported (timelimit, numerical problem etc)
0070 %  origStat     CPLEX status code. Use cplexStatus(solution.origStat) for
0071 %               more information from the CPLEX solver
0072 %  solver       solver used by cplex
0073 %  time         time taken to solve the optimization problem
0074 %
0075 %OPTIONAL OUTPUT
0076 % model.LPBasis When input basisReuse=1, we return a basis for reuse in
0077 %                   the next LP
0078 %
0079 % CPLEX consists of 4 different LP solvers which can be used to solve sysbio optimization problems
0080 % you can control which of the solvers, e.g. simplex vs interior point solver using the
0081 % CPLEX control parameter cpxControl.LPMETHOD. At the moment, the solver is
0082 % automatically chosen for you
0083 %
0084 % Ronan Fleming 23 Oct  09  ILOG-CPLEX 12.1 via  matlab API
0085 
0086 if ~exist('printLevel','var')
0087     printLevel=0;
0088 end
0089 if ~exist('basisReuse','var')
0090     basisReuse=0;
0091 end
0092 if ~exist('conflictResolve','var')
0093     conflictResolve=0;
0094 end
0095 
0096 if ~exist('minNorm','var')
0097     minNorm=0;
0098 end
0099 
0100 if basisReuse
0101     if isfield(model,'LPBasis')
0102         basis=model.LPBasis;
0103         %use advanced starting information when optimization is initiated.
0104         %cpxControl.ADVIND=1;
0105     else
0106         basis=[];
0107     end
0108 else
0109     basis=[];
0110     %do not use advanced starting information when optimization is initiated.
0111     cplex.Param.ADVIND=0;
0112 end
0113 
0114 if ~isfield(model,'A')
0115     if ~isfield(model,'S')
0116         error('Equality constraint matrix must either be a field denoted A or S.')
0117     end
0118     model.A=model.S;
0119 end
0120 
0121 if ~isfield(model,'csense')
0122     nMet=size(model.A);
0123     if printLevel>0
0124         fprintf('%s\n','Assuming equality constraints, i.e. S*v=b');
0125     end
0126     %assuming equality constraints
0127     model.csense(1:nMet,1)='E';
0128 end
0129 
0130 if ~isfield(model,'osense')
0131     %assuming maximisation
0132     model.osense=-1;
0133     if printLevel>0
0134         fprintf('%s\n','Assuming maximisation of objective');
0135     end
0136 end
0137 
0138 if size(model.A,2)~=length(model.c)
0139     error('dimensions of A & c are inconsistent');
0140 end
0141 
0142 if size(model.A,2)~=length(model.lb) || size(model.A,2)~=length(model.ub)
0143     error('dimensions of A & bounds are inconsistent');
0144 end
0145 
0146 %Conflict groups descriptor (cpxBuildConflict can be used to generate the input). Set this if
0147 %conflict refinement is desired in the case that infeasibility is detected
0148 %by CPLEX.
0149 if conflictResolve
0150     [m_lin,n]=size(model.A);
0151     m_quad=0;
0152     m_sos=0;
0153     m_log=0;
0154     %determines how elaborate the output is
0155     mode='full';%'minimal';
0156     fprintf('%s\n%s\n','Building Structure for Conflict Resolution...','...this slows CPLEX down so should not be used for repeated LP');
0157     confgrps = cpxBuildConflict(n,m_lin,m_quad,m_sos,m_log,mode);
0158     prefix=pwd;
0159     suffix='CPLEX_conflict_file.txt';
0160     conflictFile=[prefix '\' suffix];
0161 else
0162     confgrps=[]; conflictFile=[];
0163 end
0164 
0165 % Initialize the CPLEX object
0166 try
0167     cplex = Cplex('fba');
0168 catch ME
0169     error('CPLEX not installed or licence server not up')
0170 end
0171 
0172 % Now populate the problem with the data
0173 cplex.Model.sense = 'minimize';
0174 if model.osense==1
0175     %minimise linear objective
0176     cplex.Model.obj   = model.c;
0177 else
0178     %maximise linear objective by reversing sign
0179     cplex.Model.obj   = - model.c;
0180 end
0181 
0182 cplex.Model.lb    = model.lb;
0183 cplex.Model.ub    = model.ub;
0184 cplex.Model.A     = model.A;
0185 
0186 %cplex interface
0187 if isfield(model,'csense')
0188     %set up constant vectors for CPLEX
0189     b_L(model.csense == 'E',1) = model.b(model.csense == 'E');
0190     b_U(model.csense == 'E',1) = model.b(model.csense == 'E');
0191     b_L(model.csense == 'G',1) = model.b(model.csense == 'G');
0192     b_U(model.csense == 'G',1) =  Inf;
0193     b_L(model.csense == 'L',1) = -Inf;
0194     b_U(model.csense == 'L',1) = model.b(model.csense == 'L');
0195     cplex.Model.lhs   = b_L;
0196     cplex.Model.rhs   = b_U;
0197 else
0198     cplex.Model.lhs   = model.b;
0199     cplex.Model.rhs   = model.b;
0200 end
0201 
0202 %quadratic constraint matrix, size n x n
0203 if sum(minNorm)~=0
0204     if length(minNorm)==1
0205         % same weighting of min norm for all variables
0206         cplex.Model.Q=model.osense*speye(length(model.c))*minNorm;
0207     else
0208         if length(minNorm)~=length(model.c)
0209             error('Either minNorm is a scalar, or is an n x 1 vector')
0210         else
0211             % individual weighting of min norm for all variables
0212             cplex.Model.Q=model.osense*spdiags(minNorm,0,length(model.c),length(model.c));
0213         end
0214     end
0215 end
0216 
0217 %set the solver parameters
0218 if exist('contFunctName','var')
0219     if isstruct(contFunctName)
0220         cplex.Param=contFunctName;
0221     else
0222         if ~isempty(contFunctName)
0223             %calls a user specified function to create a CPLEX control structure
0224             %specific to the users problem. A TEMPLATE for one such function is
0225             %CPLEXParamSet
0226             %e.g. Param.lpmethod.Cur=0;
0227             cplex.Param=Param;
0228         end
0229     end
0230 end
0231 
0232 if printLevel==0
0233     cplex.DisplayFunc=[];
0234 else
0235     %print level
0236     cplex.Param.barrier.display.Cur = printLevel;
0237     cplex.Param.simplex.display.Cur = printLevel;
0238     cplex.Param.sifting.display.Cur = printLevel;
0239 end
0240 
0241 %limit the processing to 3 threads
0242 cplex.Param.threads.Cur = 3;
0243 
0244 % Optimize the problem
0245 cplex.solve();
0246 solution.origStat   = cplex.Solution.status;
0247 
0248 if printLevel>0 && solution.origStat~=1
0249     %use tomlab code to print out exit meassage
0250     [ExitText,ExitFlag] = cplexStatus(solution.origStat);
0251     solution.ExitText=ExitText;
0252     solution.ExitFlag=ExitFlag;
0253     if any(model.c~=0)
0254         fprintf('\n%s%g\n',[ExitText ', Objective '],  model.c'*cplex.Solution.x);
0255     end
0256 end
0257 
0258 if solution.origStat==1
0259     %extract the solution
0260     solution.obj        = model.osense*cplex.Solution.objval;
0261     solution.full       = cplex.Solution.x;
0262     solution.rcost      = cplex.Solution.reducedcost;
0263     solution.dual       = cplex.Solution.dual;
0264     solution.nInfeas    = NaN;
0265     solution.sumInfeas  = NaN;
0266     solution.solver     = cplex.Solution.method;
0267     solution.time       = cplex.Solution.time;
0268 else
0269     solution.time=NaN;
0270     %conflict resolution
0271     if conflictResolve ==1
0272         Cplex.refineConflict
0273         Cplex.writeConflict(suffix)
0274         if isfield(model,'mets') && isfield(model,'rxns')
0275             %this code reads the conflict resolution file and replaces the
0276             %arbitrary names with the abbreviations of metabolites and reactions
0277             [nMet,nRxn]=size(model.A);
0278             totAbbr=nMet+nRxn;
0279             conStrFind=cell(nMet+nRxn,1);
0280             conStrReplace=cell(nMet+nRxn,1);
0281             %only equality constraint rows
0282             for m=1:nMet
0283                 conStrFind{m,1}=['c' int2str(m) ':'];
0284                 conStrReplace{m,1}=[model.mets{m} ':  '];
0285             end
0286             %reactions
0287             for n=1:nRxn
0288                 conStrFind{nMet+n,1}=['x' int2str(n) ' '];
0289                 conStrReplace{nMet+n,1}=[model.rxns{n} ' '];
0290             end
0291             fid1 = fopen(suffix);
0292             fid2 = fopen(['COBRA_' suffix], 'w');
0293             while ~feof(fid1)
0294                 tline{1}=fgetl(fid1);
0295                 %replaces all occurrences of the string str2 within string str1
0296                 %with the string str3.
0297                 %str= strrep(str1, str2, str3)
0298                 for t=1:totAbbr
0299                     tline= strrep(tline, conStrFind{t}, conStrReplace{t});
0300                 end
0301                 fprintf(fid2,'%s\n', tline{1});
0302             end
0303             fclose(fid1);
0304             fclose(fid2);
0305             %delete other file without replacements
0306             %         delete(suffix)
0307         else
0308             warning('Need reaction and metabolite abbreviations in order to make a readable conflict resolution file');
0309         end
0310         fprintf('%s\n',['Conflict resolution file written to: ' prefix '\COBRA_' suffix]);
0311         fprintf('%s\n%s\n','The Conflict resolution file gives an irreducible infeasible subset ','of constraints which are making this LP Problem infeasible');
0312     else
0313         if printLevel>0
0314             fprintf('%s\n','No conflict resolution file. Perhaps set conflictResolve = 1 next time.');
0315         end
0316     end
0317 end
0318 
0319 % Try to give back COBRA Standardized solver status:
0320 %           1   Optimal solution
0321 %           2   Unbounded solution
0322 %           0   Infeasible
0323 %           -1  No solution reported (timelimit, numerical problem etc)
0324 if solution.origStat==1
0325     solution.stat = 1;
0326 else
0327     %use tomlab code to print out exit meassage
0328     [ExitText,ExitFlag] = cplexStatus(solution.origStat);
0329     solution.ExitText=ExitText;
0330     solution.ExitFlag=ExitFlag;
0331     if any(model.c~=0) && isfield(cplex.Solution,'x')
0332         fprintf('\n%s%g\n',[ExitText ', Objective '],  model.c'*cplex.Solution.x);
0333     end
0334     if solution.origStat==2
0335         solution.stat = 2;
0336     else
0337         if solution.origStat==3
0338             solution.stat = 0;
0339         else
0340             %this is a conservative view
0341             solution.stat = -1;
0342         end
0343     end
0344 end
0345 
0346 %return basis
0347 if basisReuse
0348     model.LPBasis=basis;
0349 end
0350 
0351 if sum(minNorm)~=0 && printLevel>0
0352     fprintf('%s\n','This objective corresponds to a flux with minimum Euclidean norm.');
0353     if length(minNorm)==1
0354         fprintf('%s%d%s\n','The weighting for minimising the norm was ',minNorm,'.');
0355     else
0356         fprintf('%s%d%s\n','The sum of the weighting for minimising the norm was ',sum(minNorm),'.');
0357     end
0358     fprintf('%s\n','Check that the objective is the same without minimising the norm.');
0359 end

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