solveCobraLP

PURPOSE ^

solveCobraLP Solve constraint-based LP problems

SYNOPSIS ^

function solution = solveCobraLP(LPproblem, varargin)

DESCRIPTION ^

solveCobraLP Solve constraint-based LP problems

 solution = solveCobraLP(LPproblem, parameters)

INPUT
 LPproblem Structure containing the following fields describing the LP
 problem to be solved
  A      LHS matrix
  b      RHS vector
  c      Objective coeff vector
  lb     Lower bound vector
  ub     Upper bound vector
  osense Objective sense (-1 max, +1 min)
  csense Constraint senses, a string containting the constraint sense for
         each row in A ('E', equality, 'G' greater than, 'L' less than).

OPTIONAL INPUTS
 Optional parameters can be entered using parameters structure or as
 parameter followed by parameter value: i.e. ,'printLevel',3)

 parameters    Structure containing optional parameters as fields.
               Setting parameters = 'default' uses default setting set in
               getCobraSolverParameters.
 printLevel    Printing level
               = 0    Silent (Default)
               = 1    Warnings and Errors
               = 2    Summary information 
               = 3    More detailed information
               > 10   Pause statements, and maximal printing (debug mode)
 saveInput     Saves LPproblem to filename specified in field. 
               i.e. parameters.saveInput = 'LPproblem.mat';
 minNorm       {(0), scalar , n x 1 vector}, where [m,n]=size(S); 
               If not zero then, minimise the Euclidean length 
               of the solution to the LP problem. minNorm ~1e-6 should be
               high enough for regularisation yet maintain the same value for 
               the linear part of the objective. However, this should be
               checked on a case by case basis, by optimization with and
               without regularisation.
 primalOnly    {(0),1} 1=only return the primal vector (lindo solvers)
               
 optional parameters can also be set through the
 solver can be set through changeCobraSolver('LP', value);
 changeCobraSolverParames('LP', 'parameter', value) function.  This
 includes the minNorm and the printLevel flags

OUTPUT
 solution Structure containing the following fields describing a LP
 solution
  full     Full LP solution vector
  obj      Objective value
  rcost    Reduced costs
  dual     Dual solution
  solver   Solver used to solve LP problem

  stat     Solver status in standardized form
            1   Optimal solution
            2   Unbounded solution
            0   Infeasible
           -1   No solution reported (timelimit, numerical problem etc)

  origStat Original status returned by the specific solver
  time     Solve time in seconds


 Markus Herrgard    08/29/06
 Ronan Fleming      11/12/08 'cplex_direct' allows for more refined control
                             of cplex than tomlab tomrun
 Ronan Fleming      04/25/09 Option to minimise the Euclidean Norm of internal
                             fluxes using either 'cplex_direct' solver or 'pdco'
 Jan Schellenberger 09/28/09 Changed header to be much simpler.  All parameters
                             now accessed through 
                             changeCobraSolverParams(LP, parameter,value)
 Richard Que        11/30/09 Changed handling of optional parameters to use
                             getCobraSolverParams().
 Ronan Fleming      12/07/09 Commenting of input/output
 Ronan Fleming      21/01/10 Not having second input, means use the parameters as specified in the
                             global paramerer variable, rather than 'default' parameters
 Steinn Gudmundsson 03/03/10 Added support for the Gurobi solver
 Tim Harrington     05/18/12 Added support for the Gurobi 5.0 solver

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function solution = solveCobraLP(LPproblem, varargin)
0002 %solveCobraLP Solve constraint-based LP problems
0003 %
0004 % solution = solveCobraLP(LPproblem, parameters)
0005 %
0006 %INPUT
0007 % LPproblem Structure containing the following fields describing the LP
0008 % problem to be solved
0009 %  A      LHS matrix
0010 %  b      RHS vector
0011 %  c      Objective coeff vector
0012 %  lb     Lower bound vector
0013 %  ub     Upper bound vector
0014 %  osense Objective sense (-1 max, +1 min)
0015 %  csense Constraint senses, a string containting the constraint sense for
0016 %         each row in A ('E', equality, 'G' greater than, 'L' less than).
0017 %
0018 %OPTIONAL INPUTS
0019 % Optional parameters can be entered using parameters structure or as
0020 % parameter followed by parameter value: i.e. ,'printLevel',3)
0021 %
0022 % parameters    Structure containing optional parameters as fields.
0023 %               Setting parameters = 'default' uses default setting set in
0024 %               getCobraSolverParameters.
0025 % printLevel    Printing level
0026 %               = 0    Silent (Default)
0027 %               = 1    Warnings and Errors
0028 %               = 2    Summary information
0029 %               = 3    More detailed information
0030 %               > 10   Pause statements, and maximal printing (debug mode)
0031 % saveInput     Saves LPproblem to filename specified in field.
0032 %               i.e. parameters.saveInput = 'LPproblem.mat';
0033 % minNorm       {(0), scalar , n x 1 vector}, where [m,n]=size(S);
0034 %               If not zero then, minimise the Euclidean length
0035 %               of the solution to the LP problem. minNorm ~1e-6 should be
0036 %               high enough for regularisation yet maintain the same value for
0037 %               the linear part of the objective. However, this should be
0038 %               checked on a case by case basis, by optimization with and
0039 %               without regularisation.
0040 % primalOnly    {(0),1} 1=only return the primal vector (lindo solvers)
0041 %
0042 % optional parameters can also be set through the
0043 % solver can be set through changeCobraSolver('LP', value);
0044 % changeCobraSolverParames('LP', 'parameter', value) function.  This
0045 % includes the minNorm and the printLevel flags
0046 %
0047 %OUTPUT
0048 % solution Structure containing the following fields describing a LP
0049 % solution
0050 %  full     Full LP solution vector
0051 %  obj      Objective value
0052 %  rcost    Reduced costs
0053 %  dual     Dual solution
0054 %  solver   Solver used to solve LP problem
0055 %
0056 %  stat     Solver status in standardized form
0057 %            1   Optimal solution
0058 %            2   Unbounded solution
0059 %            0   Infeasible
0060 %           -1   No solution reported (timelimit, numerical problem etc)
0061 %
0062 %  origStat Original status returned by the specific solver
0063 %  time     Solve time in seconds
0064 %
0065 %
0066 % Markus Herrgard    08/29/06
0067 % Ronan Fleming      11/12/08 'cplex_direct' allows for more refined control
0068 %                             of cplex than tomlab tomrun
0069 % Ronan Fleming      04/25/09 Option to minimise the Euclidean Norm of internal
0070 %                             fluxes using either 'cplex_direct' solver or 'pdco'
0071 % Jan Schellenberger 09/28/09 Changed header to be much simpler.  All parameters
0072 %                             now accessed through
0073 %                             changeCobraSolverParams(LP, parameter,value)
0074 % Richard Que        11/30/09 Changed handling of optional parameters to use
0075 %                             getCobraSolverParams().
0076 % Ronan Fleming      12/07/09 Commenting of input/output
0077 % Ronan Fleming      21/01/10 Not having second input, means use the parameters as specified in the
0078 %                             global paramerer variable, rather than 'default' parameters
0079 % Steinn Gudmundsson 03/03/10 Added support for the Gurobi solver
0080 % Tim Harrington     05/18/12 Added support for the Gurobi 5.0 solver
0081 
0082 
0083 %% Process arguments etc
0084 
0085 global CBTLPSOLVER
0086 if (~isempty(CBTLPSOLVER))
0087     solver = CBTLPSOLVER;
0088 else
0089     error('No solver found.  call changeCobraSolver(solverName)');
0090 end
0091 optParamNames = {'minNorm','printLevel','primalOnly','saveInput', ...
0092     'feasTol','optTol','EleNames','EqtNames','VarNames','EleNameFun', ...
0093     'EqtNameFun','VarNameFun','PbName','MPSfilename'};
0094 parameters = '';
0095 if nargin ~=1
0096     if mod(length(varargin),2)==0
0097         for i=1:2:length(varargin)-1
0098             if ismember(varargin{i},optParamNames)
0099                 parameters.(varargin{i}) = varargin{i+1};
0100             else
0101                 error([varargin{i} ' is not a valid optional parameter']);
0102             end
0103         end
0104     elseif strcmp(varargin{1},'default')
0105         parameters = 'default';
0106     elseif isstruct(varargin{1})
0107         parameters = varargin{1};
0108     else
0109         display('Warning: Invalid number of parameters/values')
0110         solution=[];
0111         return;
0112     end
0113 end
0114 [minNorm, printLevel, primalOnlyFlag, saveInput, feasTol, optTol] = ...
0115     getCobraSolverParams('LP',optParamNames(1:6),parameters);
0116 
0117 
0118 %Save Input if selected
0119 if ~isempty(saveInput)
0120     fileName = parameters.saveInput;
0121     if ~find(regexp(fileName,'.mat'))
0122         fileName = [fileName '.mat'];
0123     end
0124     display(['Saving LPproblem in ' fileName]);
0125     save(fileName,'LPproblem')
0126 end
0127 
0128 
0129 [A,b,c,lb,ub,csense,osense] = deal(LPproblem.A,LPproblem.b,LPproblem.c,LPproblem.lb,LPproblem.ub,LPproblem.csense,LPproblem.osense);
0130 
0131 % if any(any(~isfinite(A)))
0132 %     error('Cannot perform LP on a stoichiometric matrix with NaN of Inf coefficents.')
0133 % end
0134 
0135 % Defaults in case the solver does not return anything
0136 f = [];
0137 x = [];
0138 y = [];
0139 w = [];
0140 origStat = -99;
0141 stat = -99;
0142 
0143 t_start = clock;
0144 switch solver
0145     %% GLPK
0146     case 'glpk'
0147         params.msglev = printLevel; % level of verbosity
0148         params.tolbnd = feasTol; %tolerance
0149         params.toldj = optTol; %tolerance
0150         if (isempty(csense))
0151             clear csense
0152             csense(1:length(b),1) = 'S';
0153         else
0154             csense(csense == 'L') = 'U';
0155             csense(csense == 'G') = 'L';
0156             csense(csense == 'E') = 'S';
0157             csense = columnVector(csense);
0158         end
0159         %glpk needs b to be full, not sparse -Ronan
0160         b=full(b);
0161         [x,f,y,w,stat,origStat] = solveGlpk(c,A,b,lb,ub,csense,osense,params);
0162 
0163     case {'lindo_new','lindo_old'}
0164         %% LINDO
0165         if (strcmp(solver,'lindo_new'))
0166             % Use new API (>= 2.0)
0167             [f,x,y,w,s,origStat] = solveCobraLPLindo(A,b,c,csense,lb,ub,osense,primalOnlyFlag,false);
0168             % Note that status handling may change (see Lindo.h)
0169             if (origStat == 1 || origStat == 2)
0170                 stat = 1; % Optimal solution found
0171             elseif (origStat == 4)
0172                 stat = 2; % Unbounded
0173             elseif (origStat == 3 || origStat == 6)
0174                 stat = 0; % Infeasible
0175             else
0176                 stat = -1; % Solution not optimal or solver problem
0177             end
0178         else
0179             % Use old API
0180             [f,x,y,w,s,origStat] = solveCobraLPLindo(A,b,c,csense,lb,ub,osense,primalOnlyFlag,true);
0181             % Note that status handling may change (see Lindo.h)
0182             if (origStat == 2 || origStat == 3)
0183                 stat = 1; % Optimal solution found
0184             elseif (origStat == 5)
0185                 stat = 2; % Unbounded
0186             elseif (origStat == 4 || origStat == 6)
0187                 stat = 0; % Infeasible
0188             else
0189                 stat = -1; % Solution not optimal or solver problem
0190             end
0191         end
0192         %[f,x,y,s,w,stat] = LMSolveLPNew(A,b,c,csense,lb,ub,osense,0);
0193 
0194     case 'lp_solve'
0195         %% lp_solve
0196         if (isempty(csense))
0197             [f,x,y,origStat] = lp_solve(c*(-osense),A,b,zeros(size(A,1),1),lb,ub);
0198             f = f*(-osense);
0199         else
0200             e(csense == 'E') = 0;
0201             e(csense == 'G') = 1;
0202             e(csense == 'L') = -1;
0203             [f,x,y,origStat] = lp_solve(c*(-osense),A,b,e,lb,ub);
0204             f = f*(-osense);
0205         end
0206         % Note that status handling may change (see lp_lib.h)
0207         if (origStat == 0)
0208             stat = 1; % Optimal solution found
0209         elseif (origStat == 3)
0210             stat = 2; % Unbounded
0211         elseif (origStat == 2)
0212             stat = 0; % Infeasible
0213         else
0214             stat = -1; % Solution not optimal or solver problem
0215         end
0216         s = [];
0217         w = [];
0218     case 'mosek'
0219         %% mosek
0220         %if mosek is installed, and the paths are added ahead of matlab's
0221         %built in paths, then mosek linprog shaddows matlab linprog and
0222         %is used preferentially
0223         switch printLevel
0224             case 0
0225                options.Display='off';
0226             case 1
0227                 options.Display='final';
0228             case 2
0229                 options.Display='iter';
0230             otherwise
0231                 % Ask for default options for a function.
0232                 options  = optimset;
0233         end
0234                      
0235         if (isempty(csense))
0236             [x,f,origStat,output,lambda] = linprog(c*osense,[],[],A,b,lb,ub,[],options);
0237         else
0238             Aeq = A(csense == 'E',:);
0239             beq = b(csense == 'E');
0240             Ag = A(csense == 'G',:);
0241             bg = b(csense == 'G');
0242             Al = A(csense == 'L',:);
0243             bl = b(csense == 'L');
0244             clear A;
0245             A = [Al;-Ag];
0246             clear b;
0247             b = [bl;-bg];
0248             [x,f,origStat,output,lambda] = linprog(c*osense,A,b,Aeq,beq,lb,ub,[],options);
0249         end
0250         y = [];
0251         if (origStat > 0)
0252             stat = 1; % Optimal solution found
0253             f = f*osense;
0254             y = lambda.eqlin;
0255         elseif (origStat < 0)
0256             stat = 0; % Infeasible
0257         else
0258             stat = -1; % Solution did not converge
0259         end
0260         
0261     case 'gurobi'
0262         %% gurobi
0263         % Free academic licenses for the Gurobi solver can be obtained from
0264         % http://www.gurobi.com/html/academic.html
0265         %
0266         % The code below uses Gurobi Mex to interface with Gurobi. It can be downloaded from
0267         % http://www.convexoptimization.com/wikimization/index.php/Gurobi_Mex:_A_MATLAB_interface_for_Gurobi
0268 
0269         clear opts            % Use the default parameter settings
0270         if printLevel == 0
0271            % Version v1.10 of Gurobi Mex has a minor bug. For complete silence
0272            % Remove Line 736 of gurobi_mex.c: mexPrintf("\n");
0273            opts.Display = 0;
0274            opts.DisplayInterval = 0;
0275         else
0276            opts.Display = 1;
0277         end
0278 
0279         opts.FeasibilityTol = feasTol;
0280         opts.OptimalityTol = optTol;
0281         
0282         if (isempty(csense))
0283             clear csense
0284             csense(1:length(b),1) = '=';
0285         else
0286             csense(csense == 'L') = '<';
0287             csense(csense == 'G') = '>';
0288             csense(csense == 'E') = '=';
0289             csense = csense(:);
0290         end
0291     %gurobi_mex doesn't cast logicals to doubles automatically
0292     c = double(c);
0293         [x,f,origStat,output,y] = gurobi_mex(c,osense,sparse(A),b, ...
0294                                              csense,lb,ub,[],opts);
0295         if origStat==2
0296            stat = 1; % Optimal solutuion found
0297         elseif origStat==3
0298            stat = 0; % Infeasible
0299         elseif origStat==5
0300            stat = 2; % Unbounded
0301         elseif origStat==4
0302            stat = 0; % Gurobi reports infeasible *or* unbounded
0303         else
0304            stat = -1; % Solution not optimal or solver problem
0305         end
0306         
0307     case 'gurobi5'
0308         %% gurobi 5
0309         % Free academic licenses for the Gurobi solver can be obtained from
0310         % http://www.gurobi.com/html/academic.html
0311         resultgurobi = struct('x',[],'objval',[],'pi',[]);
0312         LPproblem.A = deal(sparse(LPproblem.A));
0313         clear params            % Use the default parameter settings
0314         
0315         if printLevel == 0 
0316            params.OutputFlag = 0;
0317            params.DisplayInterval = 1;
0318         else
0319            params.OutputFlag = 1;
0320            params.DisplayInterval = 5;
0321         end
0322 
0323         params.FeasibilityTol = feasTol;
0324         params.OptimalityTol = optTol;
0325         
0326         if (isempty(LPproblem.csense))
0327             clear LPproblem.csense
0328             LPproblem.csense(1:length(b),1) = '=';
0329         else
0330             LPproblem.csense(LPproblem.csense == 'L') = '<';
0331             LPproblem.csense(LPproblem.csense == 'G') = '>';
0332             LPproblem.csense(LPproblem.csense == 'E') = '=';
0333             LPproblem.csense = LPproblem.csense(:);
0334         end
0335     
0336         if LPproblem.osense == -1
0337             LPproblem.osense = 'max';
0338         else
0339             LPproblem.osense = 'min';
0340         end
0341         
0342         LPproblem.modelsense = LPproblem.osense;
0343         [LPproblem.rhs,LPproblem.obj,LPproblem.sense] = deal(LPproblem.b,double(LPproblem.c),LPproblem.csense);
0344         resultgurobi = gurobi(LPproblem,params);
0345         
0346         if strcmp(resultgurobi.status,'OPTIMAL')
0347            stat = 1; % Optimal solution found
0348            [x,f,y] = deal(resultgurobi.x,resultgurobi.objval,resultgurobi.pi);
0349         elseif strcmp(resultgurobi.status,'INFEASIBLE')
0350            stat = 0; % Infeasible
0351         elseif strcmp(resultgurobi.status,'UNBOUNDED')
0352            stat = 2; % Unbounded
0353         elseif strcmp(resultgurobi.status,'INF_OR_UNBD')
0354            stat = 0; % Gurobi reports infeasible *or* unbounded
0355         else
0356            stat = -1; % Solution not optimal or solver problem
0357         end
0358         
0359     case 'matlab'
0360         %matlab is not a reliable LP solver
0361         if (isempty(csense))
0362             [x,f,origStat,output,lambda] = linprog(c*osense,[],[],A,b,lb,ub);
0363         else
0364             Aeq = A(csense == 'E',:);
0365             beq = b(csense == 'E');
0366             Ag = A(csense == 'G',:);
0367             bg = b(csense == 'G');
0368             Al = A(csense == 'L',:);
0369             bl = b(csense == 'L');
0370             clear A;
0371             A = [Al;-Ag];
0372             clear b;
0373             b = [bl;-bg];
0374             [x,f,origStat,output,lambda] = linprog(c*osense,A,b,Aeq,beq,lb,ub);
0375         end
0376         y = [];
0377         if (origStat > 0)
0378             stat = 1; % Optimal solution found
0379             f = f*osense;
0380             y = lambda.eqlin;
0381         elseif (origStat < 0)
0382             stat = 0; % Infeasible
0383         else
0384             stat = -1; % Solution did not converge
0385         end
0386 
0387     case 'tomlab_cplex'
0388         %% Tomlab
0389         if (~isempty(csense))
0390             b_L(csense == 'E') = b(csense == 'E');
0391             b_U(csense == 'E') = b(csense == 'E');
0392             b_L(csense == 'G') = b(csense == 'G');
0393             b_U(csense == 'G') = 1e6;
0394             b_L(csense == 'L') = -1e6;
0395             b_U(csense == 'L') = b(csense == 'L');
0396         else
0397             b_L = b;
0398             b_U = b;
0399         end
0400         tomlabProblem = lpAssign(osense*c,A,b_L,b_U,lb,ub);
0401         %Result = tomRun('cplex', tomlabProblem, 0);
0402         % This is faster than using tomRun
0403         
0404         %set parameters
0405         tomlabProblem.optParam = optParamDef('cplex',tomlabProblem.probType);
0406         tomlabProblem.QP.F = [];
0407         tomlabProblem.PriLevOpt = printLevel;
0408         
0409         %if basis is availible use it
0410         if isfield(LPproblem,'basis') && ~isempty(LPproblem.basis)
0411             tomlabProblem.MIP.basis = LPproblem.basis;
0412         end
0413         
0414         %set tolerance
0415         tomlabProblem.MIP.cpxControl.EPRHS = feasTol;
0416         tomlabProblem.MIP.cpxControl.EPOPT = optTol;
0417         
0418         %solve
0419         Result = cplexTL(tomlabProblem);
0420 
0421         % Assign results
0422         x = Result.x_k;
0423         f = osense*sum(tomlabProblem.QP.c.*Result.x_k);
0424         %        [Result.f_k f]
0425 
0426         origStat = Result.Inform;
0427         w = Result.v_k(1:length(lb));
0428         y = Result.v_k((length(lb)+1):end);
0429         basis = Result.MIP.basis;
0430         if (origStat == 1)
0431             stat = 1;
0432         elseif (origStat == 3)
0433             stat = 0;
0434         elseif (origStat == 2 || origStat == 4)
0435             stat = 2;
0436         else
0437             stat = -1;
0438         end
0439     case 'cplex_direct'
0440         %% Tomlab cplex.m direct
0441         %Used with the current script, only some of the control affoarded with
0442         %this interface is provided. Primarily, this is to change the print
0443         %level and whether to minimise the Euclidean Norm of the internal
0444         %fluxes or not.
0445         %See solveCobraLPCPLEX.m for more refined control of cplex
0446         %Ronan Fleming 11/12/2008
0447         if isfield(LPproblem,'basis') && ~isempty(LPproblem.basis)
0448             LPproblem.LPBasis = LPproblem.basis;
0449         end
0450         [solution LPprob] = solveCobraLPCPLEX(LPproblem,printLevel,1,[],[],minNorm);
0451         solution.basis = LPprob.LPBasis;
0452         solution.solver = solver;
0453 
0454     case 'lindo'
0455         error('Solver type lindo is obsolete - use lindo_new or lindo_old instead');
0456     case 'pdco'
0457         %-----------------------------------------------------------------------
0458         % pdco.m: Primal-Dual Barrier Method for Convex Objectives (16 Dec 2008)
0459         %-----------------------------------------------------------------------
0460         % AUTHOR:
0461         %    Michael Saunders, Systems Optimization Laboratory (SOL),
0462         %    Stanford University, Stanford, California, USA.
0463         %Interfaced with Cobra toolbox by Ronan Fleming, 27 June 2009
0464         [nMet,nRxn]=size(LPproblem.A);
0465         x0 = ones(nRxn,1);
0466         y0 = ones(nMet,1);
0467         z0 = ones(nRxn,1);
0468  
0469         %setting d1 to zero is dangerous numerically, but is necessary to avoid
0470         %minimising the Euclidean norm of the optimal flux. A more
0471         %numerically stable way is to use pdco via solveCobraQP, which has
0472         %a more reasonable d1 and should be more numerically robust. -Ronan
0473         d1=0; 
0474         d2=1e-6;
0475         options = pdcoSet;
0476         options.FeaTol    = 1e-12;
0477         options.OptTol    = 1e-12;
0478         %pdco is a general purpose convex optization solver, not only a
0479         %linear optimization solver. As such, much control over the optimal
0480         %solution and the method for solution is available. However, this
0481         %also means you may have to tune the various parameters here,
0482         %especially xsize and zsize (see pdco.m) to get the real optimal
0483         %objective value
0484         xsize = 1000;
0485         zsize = 10000;
0486         
0487         options.Method=2; %QR
0488         options.MaxIter=100;
0489         options.Print=printLevel;
0490         [x,y,w,inform,PDitns,CGitns,time] = ...
0491             pdco(osense*c*10000,A,b,lb,ub,d1,d2,options,x0,y0,z0,xsize,zsize);
0492         f= c'*x;
0493         % inform = 0 if a solution is found;
0494 %        = 1 if too many iterations were required;
0495 %        = 2 if the linesearch failed too often;
0496 %        = 3 if the step lengths became too small;
0497 %        = 4 if Cholesky said ADDA was not positive definite.
0498         if (inform == 0)
0499             stat = 1;
0500         elseif (inform == 1 || inform == 2 || inform == 3)
0501             stat = 0;
0502         else
0503             stat = -1;
0504         end
0505         origStat=inform;
0506     case 'mps'
0507         %% BuildMPS
0508         % This calls buildMPS and generates a MPS format description of the
0509         % problem as the result
0510         % Build MPS Author: Bruno Luong
0511         % Interfaced with CobraToolbox by Richard Que (12/18/09)
0512         display('Solver set to MPS. This function will output an MPS matrix string for the LP problem');
0513         %Get optional parameters
0514         [EleNames,EqtNames,VarNames,EleNameFun,EqtNameFun,VarNameFun,PbName,MPSfilename] = ...
0515             getCobraSolverParams('LP',{'EleNames','EqtNames','VarNames','EleNameFun','EqtNameFun','VarNameFun','PbName','MPSfilename'},parameters);
0516         %split A matrix for L and E csense
0517         Ale = A(csense=='L',:);
0518         ble = b(csense=='L');
0519         Aeq = A(csense=='E',:);
0520         beq = b(csense=='E');
0521         
0522         %%%%Adapted from BuildMPS%%%%%
0523         [neq nvar]=size(Aeq);
0524         nle=size(Ale,1);
0525         if isempty(EleNames)
0526             EleNames=arrayfun(EleNameFun,(1:nle),'UniformOutput', false);
0527         end
0528         if isempty(EqtNames)
0529             EqtNames=arrayfun(EqtNameFun,(1:neq),'UniformOutput', false);
0530         end
0531         if isempty(VarNames)
0532             VarNames=arrayfun(VarNameFun,(1:nvar),'UniformOutput', false);
0533         end
0534         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0535         
0536         [solution] = BuildMPS(Ale, ble, Aeq, beq, c, lb, ub, PbName,'MPSfilename',MPSfilename,'EleNames',EleNames,'EqtNames',EqtNames,'VarNames',VarNames);
0537         
0538         
0539     otherwise
0540         error(['Unknown solver: ' solver]);
0541         
0542 end
0543 if ~strcmp(solver,'cplex_direct') && ~strcmp(solver,'mps')
0544     %% Assign solution
0545     t = etime(clock, t_start);
0546     if ~exist('basis','var'), basis=[]; end
0547     [solution.full,solution.obj,solution.rcost,solution.dual,solution.solver,solution.stat,solution.origStat,solution.time,solution.basis] = ...
0548         deal(x,f,w,y,solver,stat,origStat,t,basis);
0549 end
0550 
0551 %% solveGlpk Solve actual LP problem using glpk and return relevant results
0552 function [x,f,y,w,stat,origStat] = solveGlpk(c,A,b,lb,ub,csense,osense,params)
0553 
0554 % Old way of calling glpk
0555 %[x,f,stat,extra] = glpkmex(osense,c,A,b,csense,lb,ub,[],params);
0556 [x,f,origStat,extra] = glpk(c,A,b,lb,ub,csense,[],osense,params);
0557 y = extra.lambda;
0558 w = extra.redcosts;
0559 % Note that status handling may change (see glplpx.h)
0560 if (origStat == 180 || origStat == 5)
0561     stat = 1; % Optimal solution found
0562 elseif (origStat == 182 || origStat == 183 || origStat == 3 || origStat == 110)
0563     stat = 0; % Infeasible
0564 elseif (origStat == 184 || origStat == 6)
0565     stat = 2; % Unbounded
0566 else
0567     stat = -1; % Solution not optimal or solver problem
0568 end

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